Método De Heun Con Ejemplo En Matlab 1s2o6q

Capítulo III: Métodos Numéricos para Resolver Ecuaciones Diferenciales Ordinarias

MÉTODO DE HEUN [9, 10, 12] El Método de Euler puede tener serios problemas de eficiencia en aquellos puntos donde las variables están cambiando rápidamente, debido a que el método asume que la pendiente es constante sobre todo el intervalo de tiempo (𝑡𝑖 , 𝑡𝑖+1 ). Un método para mejorar la estimación de la pendiente, emplea la determinación de dos derivadas en el intervalo (una en el punto inicial y otra en el final). Las dos derivadas se promedian después con la finalidad de obtener una mejor estimación de la pendiente en todo el intervalo. Este procedimiento, es también conocido como Método PredictorCorrector, en la Figura 3.3.1 se presenta gráficamente el método.

Figura 3.3.1.

Representación gráfica del Método de Heun. (a) Predictor y (b) Corrector. Fuente: [10].

En el Método de Euler, la pendiente al inicio de un intervalo es la ecuación 3.2.1: 𝑑𝑦 = 𝑓(𝑡, 𝑦) 𝑑𝑡

3.2.1

La solución numérica de Euler es la ecuación 3.2.6: 𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑡𝑖 , 𝑦𝑖 )ℎ

3.32.6

El método estándar de Euler se completa en la ecuación 3.2.6. Sin embargo, en el Método de Heun la 𝑦𝑖+1 calculada en la ecuación 3.2.6 no es la solución final, sino una 54

Capítulo III: Métodos Numéricos para Resolver Ecuaciones Diferenciales Ordinarias

predicción intermedia. Por lo tanto, para distinguirla se va a utilizar una notación con 0 superíndice 𝑦𝑖+1 . En este método la ecuación 3.2.6 recibe el nombre de ecuación

predictora o simplemente predictor. Se va a suponer que en lugar del lado derecho de la ecuación 3.2.1 se utiliza el promedio de 𝑓(𝑡, 𝑦) sobre el intervalo (𝑡𝑖 , 𝑡𝑖+1 ), entonces la aproximación del Método de Heun seria: 𝑦𝑖+1 = 𝑦𝑖 + (𝑓𝑖 + 𝑓𝑖+1 )

ℎ 2

3.3.1

donde, 𝑓𝑖 = 𝑓(𝑡𝑖 , 𝑦𝑖 )

3.3.2

con una definición similar para 𝑓𝑖+1 . La problema con la ecuación 3.3.1 es que 𝑓𝑖+1 no puede ser evaluada hasta que 𝑦𝑖+1 sea conocida, pero esta es precisamente la cantidad que se busca. Una salida de esta dificultad es emplear la ecuación de la solución de Euler 3.2.6 para obtener una 0 estimación preliminar de 𝑦𝑖+1 , es decir, utilizar a 𝑦𝑖+1 . Y con esta estimación se puede

calcular a 𝑓𝑖+1 como se muestra en la ecuación 3.3.3: 0 ) 𝑓𝑖+1 = 𝑓(𝑡𝑖+1 , 𝑦𝑖+1

3.3.3

y luego empleando la ecuación 3.3.1 obtener el valor de 𝑦𝑖+1 . Y así mediante la ecuación 3.2.6 se pudo dar un valor de 𝑓𝑖+1 que permite el cálculo de una estimación de la pendiente al final del intervalo. La ecuación 3.3.1 es conocida como ecuación correctora o simplemente corrector y es equivalente a integrar la ecuación 3.2.1 con la regla trapezoidal. Entonces, se obtiene la siguiente descripción del proceso Predictor-Corrector: 0 𝐸𝑢𝑙𝑒𝑟 𝑝𝑟𝑒𝑑𝑖𝑐𝑡𝑜𝑟: 𝑦𝑖+1 = 𝑦𝑖 + 𝑓(𝑡𝑖 , 𝑦𝑖 )ℎ

0 )] 𝑇𝑟𝑎𝑝𝑒𝑧𝑜𝑖𝑑𝑎𝑙 𝑐𝑜𝑟𝑟𝑒𝑐𝑡𝑜𝑟: 𝑦𝑖+1 = 𝑦𝑖 + [𝑓(𝑡𝑖 , 𝑦𝑖 ) + 𝑓(𝑡𝑖+1 , 𝑦𝑖+1

3.3.4

ℎ 2

3.3.5 55

Capítulo III: Métodos Numéricos para Resolver Ecuaciones Diferenciales Ordinarias

Este algoritmo también es veces llamado el Método de Euler Modificado. Por otro lado, se observa que ningún algoritmo puede ser probado como solo predictor o corrector. Así que muchos otros métodos pueden ser clasificados como predictor corrector. Por propósitos de comparación con los Métodos Runge-Kutta, se va a expresar el Método de Heun como se muestra a continuación: 𝑦𝑖+1 = 𝑦𝑖 + (𝑘1 + 𝑘2 )

ℎ 2

3.3.6

𝑘1 = 𝑓(𝑡𝑖 , 𝑦𝑖 )

3.3.7

𝑘2 = 𝑓(𝑡𝑖 + ℎ, 𝑦𝑖 + 𝑘1 ℎ)

3.3.8

Como ejemplo, en el siguiente M-archivo tipo script se resuelve la ecuación diferencial 𝑦̇ = 𝑟𝑦 y se grafica la solución sobre el intervalo de 0 ≤ 𝑡 ≤ 0.5 para el caso donde 𝑟 = −10 y la condición inicial es 𝑦(0) = 2. La constante de tiempo es 𝜏 = − 1⁄𝑟 = 0.1, y la solución analítica es 𝑦(𝑡) = 2𝑒 −10𝑡 . Para ilustrar el efecto del tamaño de paso sobre la exactitud de la solución, se emplea un tamaño de paso de ∆𝑡 = 0.02, el cual es 20 por ciento de la constante de tiempo. r=-10; delta=0.02; y(1)=2; k=0; for time=[delta:delta:0.5] k=k+1; x(k+1)=y(k)+delta*r*y(k); y(k+1)=y(k)+(delta/2)*(r*y(k)+r*x(k+1)); end t=[0:delta:0.5]; y_true=2*exp(-10*t); plot(t,y,'o',t,y_true), xlabel('t'), ylabel('y') 56

Capítulo III: Métodos Numéricos para Resolver Ecuaciones Diferenciales Ordinarias

La Figura 3.3.2 muestra el resultado. La solución numérica es mostrada por medio de los pequeños círculos y la solución verdadera (solución analítica) por la línea sólida . Hay menos error que con el Método de Euler utilizando el mismo tamaño de paso.

Figura 3.3.2.

Solución con Método de Heun de 𝑦̇ = −10𝑦 con la condición inicial 𝑦(0) = 2. Fuente: elaboración propia

Para observar como el Método de Heun funciona con una solución oscilatoria, se considera la ecuación 𝑦̇ = sin 𝑡 con la condición inicial 𝑦(0) = 0 sobre el intervalo 0 ≤ 𝑡 ≤ 4𝜋. Para poder comparar los resultados con los obtenidos en el Método de Euler, se puede emplear un tamaño de paso de ∆𝑡 = 2𝜋⁄13. Se debe notar que el lado derecho de la E.D.O. no está en función de y, así que por lo tanto no se requiere del predictor Euler de la ecuación 3.3.4. El M-archivo tipo script programado es el siguiente:

57

Capítulo III: Métodos Numéricos para Resolver Ecuaciones Diferenciales Ordinarias

delta=2*pi/13; y(1)=0; k=0; for time=[delta:delta:4*pi] k=k+1; y(k+1)=y(k)+(delta/2)*(sin(time -delta)+sin(time)); end t_true=[0:delta/10:4*pi]; y_true=1-cos(t_true); t=[0:delta:4*pi]; plot(t,y,'o',t_true,y_true), xlabel('t'), ylabel('y'), axis([0 14 0 2]) Los resultados se muestran en la Figura 3.3.3. El error es menor que con el Método de Euler, pero aún se notan errores cerca de los picos, donde la función cambia más rápido. Este error puede ser disminuido por medio de la disminución del tamaño de paso.

Figura 3.3.3. Solución con el Método de Euler Modificado de 𝑦̇ = sin 𝑡, con la condición inicial 𝑦(0) = 0. Fuente: Elaboración propia.

58