UNIVERSIDAD TECNOLÓGICA DE LOS ANDES
FACULTAD DE INGENIERÍA
Escuela Profesional de Ingeniería Ambiental y Recursos Naturales
Ejercicios Rk simple, RK segundo orden y Heun
CURSO: simulación y modelamiento ambiental
DOCENTE: Yhon fuentes Huamán
PRESENTADO POR:
AGRADA VALLE JUSTINO
Abancay - Apurímac – Perú 2023
Método de RK simple, Euler simple
Con el método clásico de Rk simple con h = 0.5 resuelva el siguiente problema el
valor inicial en el intervalo de x=0 a 2:
𝐷𝑦
= y𝑥 2 – 1.2y
𝐷𝑥
Donde y (0) = 1
Sol:
𝐷𝑦
= F (x, y) = y𝑥 2 – 1.2y
𝐷𝑥
1
𝑌2 = 𝑌1 + (𝐾1 + 2𝐾2 + 2𝐾3 + 𝐾4 ) h
6
𝐾1 = hʄ (x, y) = y𝑥 2 – 1.2y = (1) (0)2 – 1.2(1) = -1.2
1 1
𝐾2 = ʄ (𝑥1 + 2 h, 𝑦1 + 2 𝐾1 h
Y = (0 + 0.5(0.5))2 – 1.2y = y (0.0625) – 1.2y
= (𝑦1 + 0.5𝐾1 h) (0.0625) – 1.2 (𝑦1 + 0.5𝑘1 h)
(1 + 0.5(-1.2) (0.5)) (0.0625) – 1.2(1 + 0.5 (- 1.2) (0.5)) = - 0.79625
1 1
𝑘3 = ʄ (𝑥1 + 2 h, 𝑦1 + 2 𝑘2 h)
(1 + 0.5(- 0.79625) (0.5) (0 + 0.5))2 – 1.2(1 + 0.5(- 0.79625) (0.5)) = -0.91106641
𝑘4 = ʄ (𝑥1 + h, 𝑦1 + 𝑘3 h)
(1 + (-0.91106641) (0.5)) (0 + 0.5)2 – 1.2(1 + (-0.91106641) (0.5)) = -0.51724346
1
𝑦2 = 𝑦1 + 6 (𝑘1 + 2𝑘2 + 2𝑘3 + 𝑘4 ) h
1
= 1 + 6 (-1.2 + 2(-0.79625) + 2(-0.91106641) -0.51724346) (0.5) = 0.57234364
Ejercicio 1: 𝟐𝒙𝟐 − 𝟓𝒙 + 𝟏𝟎, 𝒙 = 𝟎
𝑦=0 , ℎ=5
Sol:
dy
= 2𝑥 2 − 5𝑥 + 10
𝑑𝑥
∫ dy = ∫ 2x 2 dx − ∫ 5xdx + ∫ 10dx
2𝑥 3 5𝑥 2
𝑦=0 − + 10𝑥 + 𝑐
3 2
𝑐=4
2 3 5 2
𝑦= 𝑥 − 𝑥 + 10 + 4
3 4
𝑥=0 𝑦=4
2 5
𝑥 = 0.5 𝑦= (0.5)3 − (0.5)2 + 10 ∗ 0.5 + 4
3 2
𝑦 = 8.458
2 5
𝑥=1 , 𝑦 = (1)3 − (1)2 + 10(1) + 4
3 2
𝑦 = 12.166
Sol numérico:
yi+1 = yi + f(x, t) ∗ h
𝑦0 = 4
𝑖 = 0 → 𝑦1 = 𝑦0 + (2𝑥 2 − 5𝑥 + 10) ∗ ℎ
𝑥=0
𝑦1 = 4 + (10) ∗ 0.5 = 9.00
𝑥 = 0.5
𝑦2 = 𝑦1 + (2𝑥 2 − 5𝑥 + 10) ∗ 0.5
𝑦2 = 9 + (2 ∗ 0.52 − 5 ∗ 0.5 + 10) ∗ 0.5
𝑦2 = 13
Ejercicio 2: Resolver el problema de valor inicial
𝑑𝑦
= (y + 1) (x + 1) cos (𝑥 2 + 2x), y (0) =4
𝑑𝑥
Sol:
Primer paso
1 1
𝑘1 = ʄ (0,4) = 5, 𝑘2 = ʄ (0 + 2 * 0.5, 4 + 2 * 5 * 0.5) = 6.609
1
𝑘3 = ʄ (0 + 2 * 0.5, 4 + 7.034 * 0.5) = 7.034
𝑘4 = ʄ (0 + 0.5, 4 + 7.034 * 0.5) = 4.028
𝑥1 = 0 + 0.5 = 0.5
1
𝑦1 = 4 + * (5 + 2 * 6.609 + 2 * 7.034 + 4.028) * 0.5 = 7.026
6
7.036 − 7.026
𝐸𝑇 = = 0.1%
7.036
Segundo paso
1 1
𝑘1 = ʄ (.5, 7.026) = 3.796, 𝑘2 = ʄ (.5 + 2 * 0.5, 7.026 + 2 * 3.796 * .5) = -7.415
1 1
𝑘3 = ʄ (0.5 + 2 * 0.5, 7.026 + 2 * (-7.415) * 0.5) = -5.1
𝑘4 = ʄ (0.5 + 0.5, 7.026 + (-5.1) * 0.5) = -10.84
𝑥1 = 0.5 + 0.5 = 1.0
1
𝑦1 = 7.026 + 6 * (3.796 + 2 * (-7.415) + 2 * (-5.1) + (-10.84)) * 0.5 = 4.353
4.366 − 4.353
𝐸𝑇 = = 0.3%
4.366
Método RK segundo orden:
yi+1 = yi + (a1 k1 + a2 k 2 ) ∗ h
Donde:
k1 = f(xi , yi )
𝑘2 = 𝑓(𝑥1 + 𝑝1 𝑘1 , 𝑦𝑖 + 𝑞11 𝑘1 ∗ ℎ)
𝑒𝑛𝑡𝑜𝑛𝑐𝑒𝑠: 𝑎1 + 𝑎2 = 1
1
𝑎2 𝑝1 =
2
1
𝑎2 𝑞11 =
2
Ejercicio 1:
𝑑𝑦
= 𝑓(𝑦 + 1) (𝑥 + 1) cos(𝑥 2 + 2𝑥) , 𝑦(0) = 4
𝑑𝑥
𝒔𝒐𝒍:
Primer paso
k1 = f(0,4) = 5
3 3
𝑘2 = 𝑓 (0 + ∗ 0.5,4 + ∗ 5 ∗ 0.5) = 𝑓(0.375,5.875) = 5.945
4 4
𝑥1 = 0 + 0.5 = 0.5
1 2
𝑦1 = 4 + ( ∗ 5 + ∗ 5.945) ∗ 0.5 = 6.815
3 3
7.036 − 6.815
𝜖𝜏 = = 3.1%
7.036
Segundo paso
𝒌𝟏 = 𝒇(𝟎. 𝟓, 𝟔. 𝟖𝟏𝟓) = 𝟑. 𝟔𝟗𝟔
3 3
𝑘2 = 𝑓 (0.5 + ∗ 0.5, 6.815 + ∗ 3.696 ∗ 0.5) = 𝑓(0.875, 8.201) = −13.98
4 4
𝑥1 = 0.5 + 0.5 = 1.0
1 2
𝑦1 = 6.815 + ( ∗ 3.696 + ∗ (−13.98)) ∗ 0.5 = 2.771
3 3
4.366 − 2.771
𝜖𝜏 = = 57%
2.771
Sol analítico:
1
𝑦(𝑥) = 5 exp ( 𝑠𝑒𝑛(𝑥 2 + 2𝑥)) − 1
2
Ejercicio 2:
Dada la ecuación diferencial:
𝑑𝑦 𝑥 2
𝑒 , 𝑦(0) = 1
𝑑𝑥
𝑎𝑝𝑟𝑜𝑥𝑖𝑚𝑎𝑟 𝑌(0), ℎ = 0.1
Sol:
2
𝑓(𝑥, 𝑦) = 𝑒 𝑥 , 𝑥0 = 0; 𝑦0 = 1; ℎ = 0.1
𝑥𝑖+1 = 𝑥𝑖 + ℎ
1 2
𝑦𝑖+1 = 𝑦𝑖 + ( 𝑘1 + 𝑘2 ) ∗ ℎ
3 3
𝒅𝒐𝒏𝒅𝒆: 𝑘1 = 𝑓(𝑥𝑖 ; 𝑦𝑖 )
3 3
𝑘2 = 𝑓 (𝑥𝑖 + ℎ; 𝑦𝑖 + 𝑘1 ∗ ℎ)
4 4
Interacción:
i = 0:
2
k1 = f(x0 ; y0 ) = ex0 = 1
3 3 3 3
𝑘2 = 𝑓 (𝑥0 + ℎ; 𝑦0 + 𝑘1 ℎ) = 𝑓 (0 + (0.1); 1 + (1)(0.1)) = 𝑓(0.075; 1.075)
4 4 4 4
= 1.00564
1 2 1 2
𝑦1 = 𝑦0 + ( 𝑘1 + 𝑘2 ) ℎ = 1 + ( (1) + (1.00564)) (0.1) = 1.100376
3 3 3 3
𝑝𝑜𝑟 𝑙𝑜 𝑡𝑎𝑛𝑡𝑜 𝑦(𝑥1 ) = 𝑦(0.1) = 𝑦1 = 1.100376
Método de Heun
Ejercicio 1:
𝑑𝑦
Integrar 𝑑𝑥 = 0.05𝑦 + 4𝑒 0.8𝑥 ; 𝑥 = 0 ℎ𝑎𝑠𝑡𝑎 𝑥 = 4, con un tamaño de paso igual a
1. X=0; y=2
Sol:
4
𝑦= (𝑒 0.8𝑥 − 𝑒 −0.5𝑥 ) + 2𝑒 −0.5𝑥
1.3
4
𝑦= (𝑒 0.8(0) − 𝑒 −0.5(0) ) + 2𝑒 −0.5(0) = 2
1.3
4
𝑦= (𝑒 0.8(1) − 𝑒 −0.5(1) ) + 2𝑒 −0.5(1) = 6.1946314
1.3
4
𝑦= (𝑒 0.8(2) − 𝑒 −0.5(2) ) + 2𝑒 −0.5(2) = 14.8439219
1.3
4
𝑦= (𝑒 0.8(3) − 𝑒 −0.5(3) ) + 2𝑒 −0.5(3) = 33.6771718
1.3
4
𝑦= (𝑒 0.8(4) − 𝑒 −0.5(4) ) + 2𝑒 −0.5(4) = 75.3389626
1.3
0
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 )ℎ
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦𝑖+1 = 𝑦𝑖 + ℎ
2
𝒊=𝟎
𝑥0 = 0 ; 𝑥1 = 1 ; 𝑦0 = 2 ; ℎ = 1
0
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 ) ℎ
𝑦10 = 2 + (4𝑒 0,8(0) − 0.5(2))1
𝑦10 = 5
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦𝑖+1 = 𝑦𝑖 + ℎ
2
(4𝑒 0,8(0) − 0.5(2)) + (4𝑒 0,8(1) − 0.5(5))
𝑦1 = 2 + [ ]1
2
𝑦1 = 6.7010
𝒊=𝟏
𝑥1 = 1 ; 𝑥2 = 2 ; 𝑦1 = 6.7010 ; ℎ = 1
0
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 ) ℎ
𝑦20 = 6.7010 + (4𝑒 0,8(1) − 0.5(6.7010))1
𝑦20 = 12.2527
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦𝑖+1 = 𝑦𝑖 + ℎ
2
(4𝑒 0,8(1) − 0.5(6.7010)) + (4𝑒 0,8(2) − 0.5(12.2527))
𝑦2 = 6.7010 + [ ]1
2
𝑦2 = 16.3197
𝒊=𝟐
𝑥2 = 2 ; 𝑥3 = 3 ; 𝑦2 = 16.3197 ; ℎ = 1
0
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 ) ℎ
𝑦30 = 16.3197 + (4𝑒 0,8(1) − 0.5(16.3197))1
𝑦30 = 27.9720
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦𝑖+1 = 𝑦𝑖 + ℎ
2
(4𝑒 0,8(2) − 0.5(16.3197)) + (4𝑒 0,8(3) − 0.5(27.9720))
𝑦3 = 16.3197 + [ ]1
2
𝑦3 = 37.1992
𝒊=𝟑
𝑥3 = 3 ; 𝑥4 = 4 ; 𝑦3 = 37.1992 ; ℎ = 1
0
𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑥𝑖 , 𝑦𝑖 ) ℎ
𝑦40 = 37.1992 + (4𝑒 0,8(1) − 0.5(37.1992))1
𝑦40 = 62.6923
0
𝑓(𝑥𝑖 , 𝑦𝑖 ) + 𝑓(𝑥𝑖+1 , 𝑦𝑖+1 )
𝑦𝑖+1 = 𝑦𝑖 + ℎ
2
(4𝑒 0,8(3) − 0.5(37.1992)) + (4𝑒 0,8(4) − 0.5(62.6923))
𝑦4 = 37.1992 + [ ]1
2
𝑦4 = 83.337
𝑖 𝑌𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑖𝑣𝑜 𝑌𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑖𝑣𝑜 𝐸𝑟𝑟𝑜𝑟 (%)
0 2.0000 2.0000 0.00
1 6.1946 6.7010 8.18
2 14.8439 16.3197 9.94
3 33.6771 37.1992 10.46
4 75.3389 83.337 10.62
Ejercicio 2: 𝒚𝟏𝒊 = 𝒇(𝒙𝒊 − 𝒚𝒊 )
𝒚𝟎(𝒊+𝟏) = 𝒚𝒊 + 𝒇(𝒙𝒊 − 𝒚𝒊 ) ∗ 𝒉 𝒍𝒂 𝒆𝒄𝒖𝒂𝒄𝒊𝒐𝒏 𝒅𝒆𝒍 𝒑𝒓𝒆𝒅𝒊𝒄𝒕𝒐𝒓
𝒚𝟏(𝒊+𝟏) = 𝒇(𝒙𝒊 + 𝟏, 𝒚𝟎𝒊+𝟏 )
𝒇(𝒙𝒊 − 𝒚𝒊 ) + 𝒇(𝒙𝒊 + 𝟏, 𝒚𝟎𝒊+𝟏 )
𝒚𝒊+𝟏 = 𝒚𝒊 + ∗𝒉
𝟐
Ejemplo: con el método heun integre 𝒚′ = 𝟒𝒆𝟎.𝟖𝒙 − 𝟎. 𝟓 y desde 4, con un tamaño de paso
igual a 1. La condición inicial es en x=0, y =2.
SOL:
4
𝑦= (𝑒 0.8𝑥 − 𝑒 0.5𝑥 ) + 2𝑒 −0.5𝑥
1.3
Para calcular 𝒚𝟎 , 𝒓𝒆𝒆𝒎𝒑𝒍𝒂𝒛𝒂𝒎𝒐𝒔 𝒍𝒐𝒔 𝒗𝒂𝒍𝒐𝒓𝒆𝒔 𝒊𝒏𝒊𝒄𝒊𝒂𝒍𝒆𝒔
i=0
𝑥0 = 0, 𝑦0 = 2
𝑦 ′ 0 = 𝑓(𝑥0 . 𝑦0 )
4
𝑦′0 = (𝑒 0.8(0) − 0.5(2)) = 3
1.3
Calculando en función al predictor
0
𝑦(1) = 𝑦0 + 𝑓(𝑥0 . 𝑦0 ) ∗ (1)
0
𝑦(1) = 2 + 3 ∗ (1) = 5
Ahora reemplazamos en la siguiente:
0
𝑦 ′ 𝑖+1 = 𝑓(𝑥(𝑖+1) . 𝑦(𝑖+1) )
𝑥=1
𝑌 ′1 = 𝑓(𝑥1 . 𝑦10 ) = 𝑓(1,5)
𝑓(1,5) = 4𝑒 0.8(1) − 0.5(5) = 6.40021164
0
𝑓(𝑥𝑖 . 𝑦𝑖 ) + 𝑓(𝑥(𝑖+1) . 𝑦(𝑖+1) )
𝑦(𝑖+1) = 𝑦𝑖 + ∗ℎ
2
ℎ ℎ
0
)
𝑓(𝑥𝑖 . 𝑦𝑖 + 𝑓(𝑥(𝑖+1) . 𝑦(𝑖+1) ) 𝑓(𝑥0 . 𝑦0 ) + 𝑓(𝑥1 . 𝑦10 )
=
2 2
3 + 6.40021164
= = 4.701082
2
𝑦1 = 𝑦0 + 4.701082 ∗ (1) = 2 + 4.701082
= 6.701082
Cuando I=1
𝑥1 = 1, 𝑦1 = 6.701082
𝑦 ′1 = 𝑓(𝑥1 . 𝑦1 )
𝑦 ′1 = 4𝑒 0.8(1) − 0.5(6.701082) = 5.521
𝑪𝑨𝑳𝑪𝑼𝑳𝑨𝑵𝑫𝑶 𝑬𝑵 𝑭𝑼𝑵𝑪𝑰𝑶𝑵 𝑫𝑬 𝑷𝑹𝑬𝑫𝑰𝑪𝑻𝑶𝑹:
0
𝑦(2) = 𝑦1 + 𝑓(𝑥1 . 𝑦1 ) ∗ (1)
0
𝑦(2) = 6.701082 + 5.521 ∗ (1) = 12.25233
Ahora reemplazamos en la siguiente:
0
𝑦 ′ 𝑖+1 = 𝑓(𝑥(𝑖+1) . 𝑦(𝑖+1) )
𝑥=2
𝑌 ′ 2 = 𝑓(𝑥2 . 𝑦20 ) = 𝑓(2,12.25233)
𝑓(2,12.25233) = 4𝑒 0.8(2) − 0.5(12.25233) = 13.68596
ℎ
𝑓(𝑥1 . 𝑦1 ) + 𝑓(𝑥2 . 𝑦20 ) 5.521 + 13.68596
=
2 2
= 9.603482
𝑦2 = 𝑦1 + 9.603482 ∗ (1) = 6.701082 + 9.603482 = 𝟏𝟔. 𝟑𝟎𝟒𝟓𝟔
Eje
mpl
ific
and
o en MATLAB
%método heun
%datos de entrada
clc; clear all;
x0=0;%x0 inicial
y0=2;
h=1; %tamaño de paso
y1=y0;
xf=8;
hi=0;
fprintf ("x \t y(verdadero) \t y(Heun)\t Error\n");
i=1;
for xi=x0: xf
dy=f(xi, yi);
%para extrapolar linealmente yi+1:
xi_1=xi+h;
y0i_1=yi+f(xi,yi)*h;
dyi_1=f(xi_1,y0i_1);
yi_1=yi+((dy+dyi_1)/2)*h;
Ea=abs((yi_1-yi)/yi_1); %cálculo del error
fprintf("%.of \t %.4f \t %.4f \t %.4f \n",h,fv(xi),yi,Ea);
pointgraf(i,1)=xi;%valores de x
pointgraf(i,2)=fv(xi); %y verdadero
pointgraf(1,3)=yi; %y Heun
pointgraf(1,4)=Ea; %Error
hi=hi+h;
yi=yi_1;
end
plot(pointgraf(:,1),pointgraf(:,2)'r') % grafica de solución real
hold on
plot(pointgraf(:1),pointgraf(:,3),'b') % grafixa solucion Eun
title('APLICACION DEL METODO EUN')
xlabel('iteracion con paso h')
ylabel('funciones')
legend(('y=solucion analitica','y=solucion
Heun'),'location','southweat')
