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 qu e la pendien pen diente te es e s consta co nstante nte sobre sob re todo to do el inter i ntervalo valo de tiempo tiem po ( , + ). Un método para mejorar la estimación de la p endiente, 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éto como Método do Predicto Pred ictorrCorrector , 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:
+ = ( , )ℎ
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 + 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 inte rmedia. Por lo tanto, para distinguirla se va a utilizar una nota ción con superíndice + . 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 ( , + ), entonces la aproximación del Método de Heun seria:
+ = ( + ) donde,
ℎ
3.3.1
2
= ( , )
3.3.2
con una definición similar para + . La problema con la ecuación 3.3.1 es que + no puede ser evaluada hasta que + 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 estimación preliminar de + , es decir, utilizar a + . Y con esta estimación se puede
calcular a + como se muestra en la ecuación 3.3.3: ) + = (+ , +
3.3.3
y luego empleando la ecuación 3.3.1 obtener el valor de + . Y así mediante la ecuación 3.2.6 se pudo dar un valor de + 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: : + = ( , )ℎ
)] : + = [ ( , ) (+ , +
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 predictorcorrector. 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:
+ = ( )
ℎ 2
3.3.6
= ( , )
3.3.7
= ( ℎ, ℎ)
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 − . 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);
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 ve rdadera (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]; 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
