El método de Heun es un método implícito $$ \begin{equation*} y_{n+1} = y_{n} + \frac{h}{2} [f(t_{n}, y_{n}) + f(t_{n+1}, y_{n+1})] \end{equation*} $$
Puede hacerse una aproximación en forma iterativa con el corrector $$ \begin{align*} \tilde{y}_{n+1} &= y_{n} + h f(t_{n}, y_{n}) & \mbox{predictor} \\ y_{n+1} &= y_{n} + \frac{h}{2} [f(t_{n}, y_{n}) + f(t_{n+1}, \tilde{y}_{n+1})] & \mbox{corrector} \end{align*} $$
Resolver $$ \begin{align*} y' &= 4 e^{0.8 x} - 0.5 y \\ y(0) &= 2 \\ y(4) &= ? \end{align*} $$
La solución analítica es $$ \begin{align*} y &= \frac{4}{1.3} (e^{0.8 x} - e^{-0.5 x}) + 2 e^{-0.5 x} \\ y(4) &= \frac{4}{1.3} [e^{0.8 (4)} - e^{-0.5 (4)}] + 2 e^{-0.5 (4)} = 75.33 \end{align*} $$
Para la solución numérica elegimos \( n = 4 \) $$ \begin{equation*} h = \frac{x_{f} - x_{i}}{n} = \frac{4 - 0}{4} = 1 \end{equation*} $$
Iteracion 0 $$ \begin{align*} y_{0} &= 2 \\ x_{0} &= 0 \end{align*} $$
Iteración 1
Se calculará una iteración con el corrector $$ \begin{align*} \tilde{y}_{1} &= y_{0} + h f(x_{0}, y_{0}) = 2 + 1[4 e^{0.8 (0)} - 0.5 (2)] = 5 \\ x_{1} &= x_{0} + h = 0 + 1 = 1 \\ y_{1} &= y_{0} + \frac{h}{2} [f(x_{0}, y_{0}) + f(x_{1}, \tilde{y}_{1})] = 2 + \frac{1}{2} [4 e^{0.8 (0)} - 0.5 (2) + 4 e^{0.8 (1)} - 0.5 (5)] = 6.701 \end{align*} $$
Iteración 2 $$ \begin{align*} \tilde{y}_{2} &= y_{1} + h f(x_{1}, y_{1}) = 6.701 + 1[4 e^{0.8 (1)} - 0.5 (6.701)] = 12.253 \\ x_{2} &= x_{1} + h = 1 + 1 = 2 \\ y_{2} &= y_{1} + \frac{h}{2} [f(x_{1}, y_{1}) + f(x_{2}, \tilde{y}_{2})] = 6.701 + \frac{1}{2} [4 e^{0.8 (1)} - 0.5 (6.701) + 4 e^{0.8 (2)} - 0.5 (12.253)] = 16.320 \end{align*} $$
Iteración 3 $$ \begin{align*} \tilde{y}_{3} &= y_{2} + h f(x_{2}, y_{2}) = 16.320 + 1[4 e^{0.8 (2)} - 0.5 (16.320)] = 27.972 \\ x_{3} &= x_{2} + h = 2 + 1 = 3 \\ y_{3} &= y_{2} + \frac{h}{2} [f(x_{2}, y_{2}) + f(x_{3}, \tilde{y}_{3})] = 16.320 + \frac{1}{2} [4 e^{0.8 (2)} - 0.5 (16.320) + 4 e^{0.8 (3)} - 0.5 (27.972)] = 37.199 \end{align*} $$
Iteración 4 $$ \begin{align*} \tilde{y}_{4} &= y_{3} + h f(x_{3}, y_{3}) = 37.199 + 1[4 e^{0.8 (3)} - 0.5 (37.199)] = 62.692 \\ x_{4} &= x_{3} + h = 3 + 1 = 4 \\ y_{4} &= y_{3} + \frac{h}{2} [f(x_{3}, y_{3}) + f(x_{4}, \tilde{y}_{4})] = 37.199 + \frac{1}{2} [4 e^{0.8 (3)} - 0.5 (37.199) + 4 e^{0.8 (4)} - 0.5 (62.692)] = 83.338 \end{align*} $$
xi = 0
yi = 2
xf = 4
n = 4
h = (xf - xi)/n
y = yi
x = xi
for i=1 to n do
ya = y + h*(4*exp(0.8*x) - 0.5*y)
xa = x + h
y = y + (h/2)*(4*exp(0.8*x) - 0.5*y + 4*exp(0.8*xa) - 0.5*ya)
x = xa
end for
function g(x, y)
4*exp(0.8*x) - 0.5*y
end
g (generic function with 1 method)
function MétodoHeun(xi, yi, h)
yf = yi + h*g(xi, yi)
x = xi + h
y = yi + (h/2)*(g(xi, yi) + g(x, yf))
return x, y
end
MétodoHeun (generic function with 1 method)
xi = 0
yi = 2
xf = 4
n = 50
h = (xf - xi) / n
y = yi
x = xi
@printf("%s \t %s \t %17s\n", "i", "x", "y")
@printf("%s \t %s \t %17s\n", 0, x, y)
for i= 1:n
x, y = MétodoHeun(x, y, h)
@printf("%d \t %0.15f \t %0.15f\n",i , x, y)
end
i x y
0 0 2
1 0.080000000000000 2.245774783801841
2 0.160000000000000 2.503340685158105
3 0.240000000000000 2.773651554228560
4 0.320000000000000 3.057717441178962
5 0.400000000000000 3.356608578790039
6 0.480000000000000 3.671459617772590
7 0.560000000000000 4.003474131902822
8 0.640000000000000 4.353929411206959
9 0.720000000000000 4.724181562613496
10 0.800000000000000 5.115670938759900
11 0.880000000000000 5.529927916993519
12 0.960000000000000 5.968579052049326
13 1.040000000000000 6.433353627425923
14 1.120000000000000 6.926090632122263
15 1.200000000000000 7.448746191147498
16 1.280000000000000 8.003401480082484
17 1.360000000000000 8.592271155961400
18 1.440000000000000 9.217712338863757
19 1.520000000000000 9.882234180869657
20 1.600000000000000 10.588508061443616
21 1.680000000000000 11.339378450884560
22 1.760000000000000 12.137874486222342
23 1.840000000000001 12.987222306865421
24 1.920000000000001 13.890858200422342
25 2.000000000000000 14.852442612443983
26 2.080000000000001 15.875875077377975
27 2.160000000000001 16.965310131805623
28 2.240000000000001 18.125174275060736
29 2.320000000000001 19.360184046625317
30 2.400000000000001 20.675365294276887
31 2.480000000000001 22.076073711844899
32 2.560000000000001 23.568016730639457
33 2.640000000000001 25.157276854165467
34 2.720000000000001 26.850336531652445
35 2.800000000000001 28.654104672238510
36 2.880000000000001 30.575944908372595
37 2.960000000000001 32.623705724169042
38 3.040000000000001 34.805752572093120
39 3.120000000000001 37.131002109505609
40 3.200000000000002 39.608958695283206
41 3.280000000000002 42.249753295994395
42 3.360000000000002 45.064184960985713
43 3.440000000000002 48.063765036261614
44 3.520000000000002 51.260764298265208
45 3.600000000000002 54.668263200633277
46 3.680000000000002 58.300205439756176
47 3.760000000000002 62.171455058573692
48 3.840000000000002 66.297857322537212
49 3.920000000000002 70.696303617126731
50 4.000000000000002 75.384800632790643
Implementación iterativa del corrector:
function MétodoHeunIterativo(xi, yi, h)
es = 0.01
y = yi + h*g(xi, yi)
x = xi + h
while true
y_anterior = y
y = yi + (h/2)*(g(xi, yi) + g(x, y_anterior))
ea = ((y - y_anterior) / y) * 100
if ea <= es
break
end
end
return x, y
end
MétodoHeunIterativo (generic function with 1 method)
xi = 0
yi = 2
xf = 4
n = 50
h = (xf - xi) / n
y = yi
x = xi
@printf("%s \t %s \t %17s\n", "i", "x", "y")
@printf("%s \t %s \t %17s\n", 0, x, y)
for i= 1:n
x, y = MétodoHeunIterativo(x, y, h)
@printf("%d \t %0.15f \t %0.15f\n", i, x, y)
end
i x y
0 0 2
1 0.080000000000000 2.245659288125804
2 0.160000000000000 2.503104772258407
3 0.240000000000000 2.773289942176862
4 0.320000000000000 3.057224462431668
5 0.400000000000000 3.355978153450947
6 0.480000000000000 3.670685225242664
7 0.560000000000000 4.002548780800291
8 0.640000000000000 4.352845607433842
9 0.720000000000000 4.722931275436943
10 0.800000000000000 5.114245564768579
11 0.880000000000000 5.528318241780513
12 0.960000000000000 5.966775209463633
13 1.040000000000000 6.431345056224662
14 1.120000000000000 6.923866029844989
15 1.200000000000000 7.446293465022624
16 1.280000000000000 8.000707694763614
17 1.360000000000000 8.589322477878358
18 1.440000000000000 9.214493976959218
19 1.520000000000000 9.878730323477448
20 1.600000000000000 10.584701809048932
21 1.680000000000000 11.335251744489456
22 1.760000000000000 12.133408031021849
23 1.840000000000001 12.982395490920409
24 1.920000000000001 13.885649007994747
25 2.000000000000000 14.846827531638148
26 2.080000000000001 15.869829001708514
27 2.160000000000001 16.958806255287346
28 2.240000000000001 18.118183980389627
29 2.320000000000001 19.352676785991239
30 2.400000000000001 20.667308462318534
31 2.480000000000001 22.067432510225302
32 2.560000000000001 23.558754023686031
33 2.640000000000001 25.147353014981977
34 2.720000000000001 26.839709278071265
35 2.800000000000001 28.642728891939893
36 2.880000000000001 30.563772472453330
37 2.960000000000001 32.610685288395544
38 3.040000000000001 34.791829365023553
39 3.120000000000001 37.116117706611860
40 3.200000000000002 39.593050778146157
41 3.280000000000002 42.232755395584817
42 3.360000000000002 45.046026183977929
43 3.440000000000002 48.044369773257522
44 3.520000000000002 51.240051912732163
45 3.600000000000002 54.646147697280405
46 3.680000000000002 58.276595110989348
47 3.760000000000002 62.146252107579599
48 3.840000000000002 66.270957461451360
49 3.920000000000002 70.667595638637934
50 4.000000000000002 75.354165953425863