SOLUCION NUMERICA DE ECUACIONES DIFERENCIALES ORDINARIAS. Roland England Marzo 2003.

1

Contents 1 Introducci´ on 1.1 Repaso de las propiedades de ecuaciones diferenciales ordinarias y de sus soluciones. . . . . . . . . . . . . . . . . . . . . 1.1.1 Soluciones anal´ıticas y num´ericas. . . . . . . . . . . . 1.1.2 Sistemas de ecuaciones, y ecuaciones de orden mayor que uno. . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Soluciones generales y condiciones iniciales. . . . . . . 1.1.4 Singularidades y ecuaciones impl´ıcitas. . . . . . . . . . 1.1.5 Teorema de existencia. . . . . . . . . . . . . . . . . . . 1.2 Ecuaciones lineales, linealizaci´ on de ecuaciones no-lineales, y estabilidad del problema. . . . . . . . . . . . . . . . . . . . . 1.2.1 Soluci´ on general de ecuaciones lineales. . . . . . . . . 1.2.2 Aplicaci´ on a la primera variaci´ on de ecuaciones nolineales, y estabilidad del problema. . . . . . . . . . .

4 4 4 5 6 7 9 10 10 12

2 CONCEPTOS BASICOS PARA PROBLEMAS CON CONDICIONES INICIALES. 14 2.1 M´etodo de Euler con an´ alisis del error. . . . . . . . . . . . . 14 2.1.1 M´etodo de Euler. . . . . . . . . . . . . . . . . . . . . 14 2.1.2 Pruebas num´ericas. . . . . . . . . . . . . . . . . . . . . 15 2.1.3 Una cota para el error. . . . . . . . . . . . . . . . . . . 16 2.1.4 Comprobaci´ on num´erica. . . . . . . . . . . . . . . . . 19 2.1.5 Estimaci´ on del error. . . . . . . . . . . . . . . . . . . . 19 3 ESTABILIDAD Y METODOS IMPLICITOS 21 3.1 Estabilidad del problema con condiciones iniciales . . . . . . 21 3.1.1 Problemas estables con cota grande para el error. . . . 21 3.1.2 Ecuaciones Diferenciales ”stiff”, y el m´etodo de Euler. 22

2

3.1.3

Un m´etodo con regi´on de estabilidad absoluta m´as grande. . . . . . . . . . . . . . . . . . . . . . . . . . .

4 METODOS CON UN SOLO PASO 4.1 Clasificaci´ on de m´etodos . . . . . . . . . . . . . . . . . . 4.1.1 M´etodos con un solo paso y con pasos m´ ultiples. 4.1.2 El m´etodo de la serie de Taylor. . . . . . . . . . 4.1.3 M´etodos de extrapolaci´ on. . . . . . . . . . . . . . 4.1.4 M´etodos del tipo Runge-Kutta. . . . . . . . . . .

. . . . .

. . . . .

. . . . .

5 METODOS DE RUNGE-KUTTA Y CONTROL DEL ERROR. 5.1 M´etodos del tipo Runge-Kutta. . . . . . . . . . . . . . . . . . 5.1.1 Forma aut´ onoma, y notaci´ on. . . . . . . . . . . . . . . 5.1.2 M´etodo del punto medio por extrapolaci´ on. . . . . . . 5.1.3 M´etodos del tipo Runge-Kutta de orden dos. . . . . . 5.1.4 La clase de m´etodos de orden cuatro. . . . . . . . . . 5.1.5 Suposiciones que permiten el ajuste pr´ actico del paso. 5.1.6 Estimaci´ on del error de truncamiento local. . . . . . .

25 28 28 28 29 29 34

35 35 35 36 36 37 41 42

6 METODOS CON PASOS MULTIPLES 45 6.1 Introducci´ on a los m´etodos con pasos m´ ultiples. . . . . . . . . 45 6.1.1 Los m´etodos de Adams. . . . . . . . . . . . . . . . . . 45 6.1.2 Estabilidad de los m´etodos de pasos m´ ultiples. . . . . 50 7 FORMULACION DE NORDSIECK Y METODOS PARA ECUACIONES ”STIFF” 7.1 Una t´ecnica para cambiar el paso, y las f´ormulas de Gear . . 7.1.1 La formulaci´on de Nordsieck. . . . . . . . . . . . . . . 7.1.2 Estabilidad y m´etodos para ecuaciones ”stiff”. . . . .

54 54 54 59

8 METODOS AUTOMATICOS. 64 8.1 Dise˜ no de rutinas con un m´ınimo de par´ ametros de control. . 64 8.1.1 Necesidades para un algoritmo autom´ atico. . . . . . . 64

3

Chapter 1

Introducci´ on 1.1 1.1.1

Repaso de las propiedades de ecuaciones diferenciales ordinarias y de sus soluciones. Soluciones anal´ıticas y num´ ericas.

El problema m´as sencillo en ecuaciones diferenciales ordinarias es el de encontrar una funci´ on y(t) cuya derivada dy = f (t, y) (1.1) dt donde f (t, y) es una funci´ on dada de t, y. En el caso particular donde f (t, y) no depende de y, y se escribe f (t), la soluci´on formal es .

y=




y=

f (t) dt.

(1.2)

Sin embargo, atr´ as de esta notaci´on quedan dos problemas inmediatos, uno f´ acil de resolver, y el otro causa de muchas dificultades pr´ acticas. Primero, la ecuacion (1.2) no se determina de manera u ´nica, pues si se le suma una constante arbitraria al lado derecho, sigue satisfaciendo la ecuaci´ on (1.1). Se resuelve este problema fijando el valor de y = y0 cuando t = t0 . Entonces, en lugar de la integral indefinida de (1.2), se tiene la integral definida
t

y(t) = y0 + que es la soluci´on u ´nica. 4

f (τ )dτ t0

(1.3)

El segundo problema es que la integral en (1.2) o (1.3) no puede generalmente evaluarse expl´ıcitamente como una expresi´on anal´ıtica para y(t). Es posible s´ olo para ciertas formas espec´ıficas de la funci´ on f (t), y a´ un as´ı, muchas veces es necesario adivinar el resultado, y verificarlo mediante la ecuaci´on ( 1.1). Ejemplo 1.1.1 Si f (t) = 1/(1+et sin2 t) no se conoce una integral expl´ıcita, y s´ olo puede resolverse mediante integraci´ on num´erica. Ejemplo 1.1.2 Sea f (t) = 1 + et sin2 t. Entonces f (t) = 1 + 12et (1 − cos 2t). 



f (t)dt = t + 12 et 1 − 15 (cos 2t + 2 sin 2t + constante arbitrario.

Cuando f depende expl´ıcitamente de y, como en ecuaci´on (1.1), puede pensarse que el primer problema desaparece, pues al sumar una constante . a y, se cambia f (t, y), sin cambiar y . Sin embargo, al analizar el caso, hay precisamente el mismo nivel de arbitrariedad, lo que puede resolverse fijando otra vez el calor de y(t0 ). Por otra parte, es de esperarse que el segundo problema vuelva a´ un m´as dif´ıcil cuando y aparece expl´ıcitamente en f . ¡A´ un no se conoce el integrando en ecuaci´ on (1.3) antes de encontrar la soluci´ on!. El hecho es que la soluci´ on queda bien determinada, excepto en la vecindad de ciertos puntos llamados ”singularidades” donde f (t, y) tiende hacia infinito o toma una forma indefinida tal como 0/0. Sin embargo, es cierto que el problema es m´as dif´ıcil, y en los casos sencillos que pueden reducirse por an´ alisis a unas integraciones simples tales como (1.3), en el pasado se ha considerado que vale la pena el trabajo de tal reducci´ on.

1.1.2

Sistemas de ecuaciones, y ecuaciones de orden mayor que uno.

Antes de considerar la soluci´ on del problema (1.1), vale la pena mencionar dos generalizaciones. Primero, y, f pueden ser vectores de dimensi´ on n, o sea (1.1) ser´ıa una notaci´ on reducida para el sistema de ecuaciones

5

.

y 1 = f1 (t, y1 , y2 , ..., yn ) . y 2 = f2 (t, y1 , y2 , ..., yn ) .. . . y n = fn (t, y1 , y2 , ..., yn )

(1.4)

donde f1 , f2 , ..., fn son funciones distintas de los n + 1 argumentos. No todos los argumentos deben de aparecer expl´ıcitamente en todas las funciones. En segundo lugar, puede tenerse una ecuaci´on de orden mayor, en la cual . dn y on no es y sino una derivada mayor y (n) = n que se especifica, como funci´ dt de t, y, y de las derivadas de o´rdenes mayores, de esta forma: . ..

y (n) = f (t, y, y , y , ..., y (n−1) ).

(1.5)

Afortunadamente, esta generalizaci´ on no es diferente, como puede verse al escribir y1 = y ,

yi = y (i−1) ,

i = 2, 3, ..., n

(1.6)

luego se puede reemplazar la ecuaci´on (1.5) por el sistema de ecuaciones de primer orden .

y 1 = y2 . y 2 = y3 .. . . y n−1 = yn . y n = f (t, y1 , y2 , ..., yn )

(1.7)

que es una forma particularmente sencilla de las ecuaciones (1.4). Normalmente, se supone que ecuaciones de orden mayor se reducen a sistemas de primer orden para la soluci´ on num´erica, aun que algunos m´etodos se presentar´ an para tratar directamente ecuaciones de orden mayor.

1.1.3

Soluciones generales y condiciones iniciales.

Volviendo al tema del nivel de arbitrariedad de la soluci´ on de (1.1) o de (1.4), se considerar´ a expl´ıcitamente el caso (1.1) en t´erminos de curvas en el espacio de dos dimensiones con coordenadas t, y, pero todos los argumentos ser´an v´ alidos para (1.4) en un espacio de n + 1 dimensiones con coordenadas t, y1 , y2 , ..., yn . 6

Consid´erese cualquier punto inicial especificado (t0 , y0 ) [o en caso de un on f en (1.1) [o cualquiera de las sistema (t0 , y10 , y20 , ..., yno )]. Si la funci´ funciones fi en (1.4)] no est´a definida en este punto, entonces este punto no puede encontrarse en una de las soluciones de (1.1). Por otra parte, si f [o . todas las. fi ] est´a definida en este punto, entonces fija el valor de y [o de toon de movimiento del punto cuando t se aumenta. das las y i ], y as´ı la direcci´ En un punto vecino en esta direcci´ on, una direcci´ on (un poco diferente) de movimiento se determina de manera similar, y se ve intuitivamente que la curva de la soluci´ on que empieza en este punto inicial se encuentra determinada de manera u ´nica, hasta cuando, posiblemente, llega a un punto en el cual f [o una de las fi ] no est´a definida, o sea, una singularidad. Si se empieza en un punto inicial distinto, se encuentra una curva de soluci´ on distinta, que no puede tener una intersecci´ on con la primera en un punto regular, como la direcci´ on de la curva en tal punto es u ´nica. Sigue que, en el caso general, pueden especificarse arbitrariamente las n + 1 coordenadas del punto inicial (al menos dentro de la regi´ on de definici´ on de la funci´ on f ), y que para cada punto se obtiene una curva u ´nica de soluci´ on. Desde luego, hay una simple infinidad de puntos en cada una de estas curvas, lo que deja n par´ ametros para especificar la soluci´on de manera u ´nica. Adem´ as, sin considerar las posibles dificultades en el caso de la existencia de singularidades, estos par´ametros pueden ser los n valores y1 (t0 ), y2 (t0 ), ..., yn (t0 ) para un t0 fijo, arbitrario, que es posible tomar como cero, cambiando, si es necesario, el origen de t.

1.1.4

Singularidades y ecuaciones impl´ıcitas.

Para ver un poco las dificultades que pueden presentarse, se dan a continuaci´ on unos ejemplos sencillos de problemas con singularidades. .

on u ´nica por cada punto en el Ejemplo 1.1.3 Si y = −t/y hay una soluci´ cual no se tiene y = 0. En la recta y = 0, f (t, y) es infinita excepto para t = 0, cuando no esta definida. Es f´ acil verificar que las soluciones son los c´ırculos conc´entricos t2 + y 2 = c, que claramente no tienen intersecciones. Adem´ as, es posible incluir las singularidades de la recta y = 0, en las cuales los c´ırculos tienen tangente vertical, por argumentos de continuidad. Sin embargo, no hay soluci´ on que pasa por el origen. .

on u ´nica por cada punto, excepto Ejemplo 1.1.4 Si y = y/t hay una soluci´ en la recta t = 0, en donde f (t, y) es infinita (con excepci´ on del punto (0, 0)

7

en donde no est´ a definida). Estas soluciones son las rectas y = ct, que pasan todas por la singularidad en (0, 0), el u ´nico punto en el cual tienen intersecci´ on. La recta t = 0 tiene tangente vertical, y puede considerarse tambi´en como una soluci´ on. .

Ejemplo 1.1.5 Si y = −y/t hay otra vez una soluci´ on u ´nica por cada punto en el cual t
= 0, pero ahora las soluciones son las hip´erbolas rectangulares ty = c. La recta particular t = 0 puede considerarse como una soluci´ on, parte de la hip´erbola degenerada con c = 0 (que constituyen las dos rectas t = 0, y = 0). Pasan as´ı dos soluciones por la singularidad (0, 0). .

on similar Ejemplo 1.1.6 Si y = t/y se tiene a primera vista una situaci´ a la del ejemplo 1.1.3. Esta vez las soluciones son hip´erbolas rectangulares t2 − y 2 = c, y para el caso particular c = 0, esto degenera para formar las dos rectas y = ±t con intersecci´ on en el punto singular (0, 0). .

on u ´nica por cada punto Ejemplo 1.1.7 Si y = (t + y)/(t − y) hay una soluci´ excepto en la recta y = t, en la cual f (t, y) es en general infinita, pero no est´ a definida en el origen (0, 0). Escribiendo t = r cos θ y = r sin θ se obtiene .

.

r sin θ + r θ cos θ cos θ + sin θ = . . cos θ − sin θ r cos θ − r θ sin θ .

.

lo cual se simplifica para dar r= r θ, o sea dr dθ = r, cuyas soluciones son las θ espirales equiangulares r = ce . Todas estas tienen tangentes verticales en o 5π la recta y = t, o sea cuando θ = π4 ´ 4 , estos puntos no siendo singularidades esenciales. Sin embargo, aunque todas las curvas convergen hacia la singularidad en (0, 0), no hay soluci´ on en este punto. Singularidades pueden tambi´en presentarse cuando las ecuaciones diferenciales no vienen en una de las formas expl´ıcitas consideradas. Tratando el caso m´as sencillo de una sola ecuaci´on de primer orden, esto puede presentarse en la forma impl´ıcita .

g(t, y, y) = 0

8

(1.8)

donde. g es una funci´ on dada de sus tres argumentos. La resoluci´ on de (1.8) y para , para ponerla en la forma (1.1) puede presentar problemas, y no siempre es factible. Existen m´etodos, similares al m´etodo de Gear, para resolver sistemas impl´ıcitos directamente, pero todav´ıa tienen unos problemas no-resueltos, y no se presentar´an aqu´ı. Siempre se supondr´ a que se tienen las ecuaciones en la forma (1.4) , o a veces (1.5). 2 Ejemplo 1.1.8 Consid´erese la ecuaci´ on y +y 2 − 1 = 0 que se satisface si . y y = sin(t + c)para cualquier valor del constante c. Si se resuelve para se . on con dos valores para | y |< 1, un valor tiene y= ± 1 − y 2 , una funci´ para | y |= 1, y ning´ un valor real para | y |> 1. La soluci´ on general dada necesita utilizar los dos valores de la ra´ız cuadrada. Las l´ıneas y = ±1 son lineas de singularidades, en las cuales la funci´ on cambia de signo. Son las envolventes de las curvas de soluci´ on, que s´ı tienen intersecciones, como la derivada tiene dos valores posibles en cualquier punto. Las envolventes constituyen ellas mismas unas soluciones que no se incluyen en la forma general de la soluci´ on. Es posible tambi´en construir otras soluciones, con derivadas continuas, por partes de la soluci´ on general, y partes de las envolventes. De todos modos, la soluci´ on no es u ´nica, una vez que se llega a las singularidades. Adem´ as, si, por errores tales como redondeo, la aproximaci´ on para y que se obtiene de la computadora es mayor que uno en valor absoluto, su derivada ya no tiene valor real. As´ı, aun una ecuaci´ on tan sencilla puede presentar problemas en la soluci´ on num´erica. .

1.1.5

Teorema de existencia.

Hay un teorema de existencia y unicidad de la soluci´ on del problema (1.1) con una condici´ on inicial y(t0 ) dada. Esencialmente define que es un punto regular, dejando que un punto inicial que no satisface las condiciones del teorema puede ser una singularidad. Sin dar la forma m´ as comprensiva, se tiene: .

Teorema 1.1.1 Dado la ecuaci´ on diferencial y = f (t, y), donde f (t, y) es una funci´ on cont´ınua para 0 ≤ t ≤ b (todo y), y existe un constante L tal que | f (t, y) − f (t, z) |≤ L | y − z |

9

para 0 ≤ t ≤ b (todo y, z) — Condici´ on de Lipschitz — entonces existe una . y funci´ on cont´ınua, diferenciable, y u ´nica y(t) tal que (t) = f (t, y(t)) , y(0) = y0 . La demostraci´on puede encontrarse en los libros sobre la teor´ıa de ecuaciones diferenciales ordinarias. La condici´ on de Lipschitz es m´as fuerte que una condici´ on de continuidad con respecto a y, pero no requiere que f (t, y) tenga derivada cont´ınua con respecto a y, el constante de Lipschitz L es una cota para la derivada parcial ∂f ∂y . Es posible plantear un teorema an´ alogo al Teorema 1.1.1, si las condiciones se satisfacen para y, z dentro de una regi´ on finita que contiene y0 . Adem´as, puede aplicarse a sistemas de ecuaciones, donde y, z, f son vectores de dimensi´ on n, si los valores absolutos se reemplazan por normas.

1.2 1.2.1

Ecuaciones lineales, linealizaci´ on de ecuaciones no-lineales, y estabilidad del problema. Soluci´ on general de ecuaciones lineales.

Para terminar esta introducci´ on, se mencionar´a la forma particular de la soluci´ on general para un sistema lineal, o sea un sistema en el cual los yi y sus derivadas aparecen linealmente, sin potencias o productos. Tal sistema puede escribirse .

y = A(t)y + b(t)

(1.9)

en lo cual A(t) es una matriz (nxn) y b(t) un vector de n elementos, las dos funciones de t, elemento por elemento. Si se tiene una soluci´on particular y = Y (t)

(1.10)

de (1.9), y n soluciones linealmente independientes y = Zr (t) ,

r = 1, 2, ..., n

(1.11)

de la ecuaci´ on homogenea .

y = A(t)y

10

(1.12)

que se obtiene de (1.9) reemplazando b(t) por un vector de n ceros, entonces la soluci´ on m´ as general posible de (1.9) es de la forma y = Y (t) +

n 

Cr Zr (t)

(1.13)

r=1

en la cual los Cr son constantes que pueden tomar cualquier valor fijo. El primer t´ermino en el lado derecho de (1.13) se llama Integral Particular, y la sumatoria se llama Funci´ on Complementaria. Se ve que (1.13) es efectivamente una soluci´ on de (1.9) por sustituci´ on directa: .

.

y −A(t)y − b(t) =Y −A(t)Y − b(t) +

n 

.



Cr Z r −A(t)Zr = 0

(1.14)

r=1

Para demostrar que la soluci´ on general es de la forma (1.13) se considera una soluci´ on y = W (t)

(1.15)

Z0 (t) = Y (t) − W (t).

(1.16)

de (1.9), y la funci´ on

on como Se ve que Z0 (t) es una soluci´on de (1.12), y por la misma sustituci´ antes y = b0 Z0 (t) +

n 

br Zr (t)

(1.17)

r=1

es una soluci´ on de (1.12) para valores cualesquiera de los n + 1 constantes as de n constantes b0 , b1 , ..., bn . Como la soluci´on general no puede tener m´ arbitrarios, una de las funciones Z0 (t), Z1 (t), ..., Zn (t) puede expresarse como combinaci´ on lineal de las dem´ as, para permitir la eliminaci´ on de uno de los constantes en (1.17). Concretamente, debe existir una relaci´on de la forma a0 Z0 (t) +

n 

ar Zr (t) = 0

(1.18)

r=1

en la cual no todos los coeficientes a0 , a1 , ..., an son ceros. No es posible que a0 = 0, como las n funciones Z1 , Z2 , ...., Zn se tomaron linealmente independientes. Entonces, se tiene 11

W (t) = Y (t) − Z0 (t)  = Y (t) + nr=1

(1.19)

ar a0 Zr (t)

utilizando (1.18). Como W (t) se tom´o como una soluci´ on cualquiera de (1.9), se tiene que la soluci´ on general es de la forma (1.13), con Cr = ar /a0 .

1.2.2

Aplicaci´ on a la primera variaci´ on de ecuaciones nolineales, y estabilidad del problema.

Si se tiene un sistema no-lineal (1.4), puede aplicarse este resultado a errores peque˜ nos en la soluci´ on. Sup´ ongase que la soluci´ on correcta es z(t), y que en lugar de resolver (1.4), se resuelve .

y = f (t, y) + δ (t)

(1.20)

donde δ (t) es un vector de errores, y que el resultado obtenido es y (t) = z (t) + ε (t) .

(1.21)

Se puede linealizar (1.20) para obtener .



.

z (t) + ε (t) = f (t, z) + A (t) ε + δ (t) + 0 ε2 .

(1.22)

En esta ecaci´  on A (t) es la matriz Jacobiana del sistema, que tiene elemento ∂fi , y 0(ε2 ) representa t´erminos del orden ε2 en la serie de Aij = ∂yj y=z(t)

Taylor para f (t, y). Se sabe que .

z= f (t, z)

(1.23)

como z (t) es la soluci´on correcta. Restando de (1.22), y despreciando los nos, se obtiene la ecuaci´on t´erminos de orden ε2 que se suponen peque˜ .

ε= A (t) ε + δ (t)

(1.24)

una ecuaci´ on lineal para ε. Ahora es posible aplicar el resultado anterior, e introducir el concepto de estabilidad que ser´ a muy importante en este curso. Si unas soluciones de (1.24) vuelven muy grandes en comparaci´ on con z(t), 2 ya no es posible despreciar t´erminos del orden ε , pero de todos modos los errores generan una desviaci´on grande de la soluci´ on correcta z(t). En este caso, el problema matem´atico es inestable. Si todas las soluciones de (1.24) quedan chicos, el problema matem´atico es estable. En tal caso (y a veces 12

unas de las funciones complementarias de (1.24) tienen un decrecimiento muy r´ apido), es todav´ıa posible que la soluci´on num´erica presente problemas, o requiera una cantidad excesiva de computaci´ on. Este fen´ omeno se llama ”stiffness” (tiesidad o rigidez), y se relaciona con ciertas restricciones de estabilidad de los m´etodos num´ericos. .

.

Ejemplo 1.2.1 Si y = λy+ s (t) − λs (t) donde λ es un constante, y s (t) . on general puede escribirse una funci´ on dada con derivada s, la soluci´ y (t) = [y (0) − s (0)] eλt + s (t) y si y(0) = s(0), la soluci´ on es y(t) = s(t), que puede ser una funci´ on sin variaciones r´ apidas. Si λ es grande positivo, errores pueden generar una desviaci´ on hasta una soluci´ on ”vecina” para la cual el primer t´ermino no es cero, y que desv´ıa r´ apidamente de la soluci´ on correcta s(t). Tal problema tiene inestabilidad matem´ atica. Si λ es grande negativo, toda soluci´ on vecina converge r´ apidamente hacia s(t) a medida que t aumenta, y el problema es muy estable. Sin embargo, es ”stiff ”, y los m´etodos num´ericos con restricciones de estabilidad pueden requerir grandes cantidades de computaci´ on. Cuando se tiene un sistema de ecuaciones, el comportamiento puede ser mucho m´as complejo. Vale la pena comentar que un problema que es estable para integraci´ on en el sentido de crecimiento de t, puede ser inestable en el sentido de decrecimiento de t.

13

Chapter 2

CONCEPTOS BASICOS PARA PROBLEMAS CON CONDICIONES INICIALES. 2.1

M´ etodo de Euler con an´ alisis del error.

2.1.1

M´ etodo de Euler.

El problema num´erico con condiciones iniciales (lo cual se supone es bien planteado, por satisfacer las condiciones de un teorema de existencia y unicidad) es el de encontrar la soluci´ on de la ecuaci´ on (1.1) para t ≥ t0 dado y (t0 ), sobre todo para los casos para los cuales no se conoce una soluci´on anal´ıtica, o si se conoce, no es. conveniente utilizarla. En t0 , ya se tienen y (t0 ), y (t0 ) = f (t0 , y (t0 )), y por serie de Taylor y (t0 + h) = y (t0 ) + hf (t0 , y (t0 )) + 0(h2 )

(2.1)

na. puede aproximarse por y1 = y (t0 ) + hf (t0 , y (t0 )) para una h peque˜ Con esta aproximaci´on, el proceso puede repetirse en t1 = t0 + h para dar una aproximaci´ on de y (t0 + 2h), o sea y2 = y1 + hf (t1 , y1 ). As´ı, se puede proceder paso por paso para obtener una secuencia yi de aproximaciones para y (ti ) donde

14

ti = t0 + ih, yi+1 = yi + hf (ti , yi ) ,

i = 0, 1, ...

(2.2)

con una magnitud constante de paso, h. No es necesario que los pasos sean todos iguales, y de hecho, en muchos casos es preferible variar h en cada paso. Para el momento, para simplificar el estudio, se considera h constante. El m´etodo (2.2) es el m´as sencillo de las t´ecnicas num´ericas para ecuaciones diferenciales, y se llama m´etodo de Euler (expl´ıcito). Si la soluci´ on de (1.1) se requiere en t = b, puede asegurarse que el m´etodo de Euler d´e una aproximaci´ on para y (b), con h constante, y sin interpolaci´ on, tomando b−t0 h = N , donde N es el n´ umero de pasos que se tomar´an para integrar de t0 hasta b. Entonces tN = b, yN es una aproximaci´on para y (tN ). La pregunta esencial es ”¿Qu´e tan buena es esta aproximaci´on?”. O puede preguntarse si tiene alguna relaci´ on con y (b).

2.1.2

Pruebas num´ ericas.

La manera m´as sencilla de investigar esas preguntas es de aplicar el m´etodo de Euler a un problema del cual se conoce la soluci´ on anal´ıtica. Ejemplo 2.1.1

.

y= y ,

y (0) = 1.

Aplicando el m´etodo de Euler, o bien ”resolviendo num´ericamente” la ecuaci´ on para y (1) (= e ≈ 2.71828) con varios valores de h, se encuentran los siguientes resultados: h 1 1/4 1/16 1/64

yN error (yN − y (1)) 2 −0.71828 2.44141 −0.27688 2.63793 −0.08035 2.69734 −0.02094

Se observa que el error se reduce con h. El error depende tambi´en de la ecuaci´ on diferencial. En general, y1 no coincide con la soluci´ on correcta en t1 , sino se encuentra sobre la curva de una soluci´ on vecina. Cuando i crece, los yi normalmente se encuentran en curvas de soluciones distintas. Si curvas de soluci´ on que se encuentran cerca unas a otras en ciertos valores de ti , i = 0, 1, ..., ellas se desv´ıan para tener una separaci´ on grande en tN ; as´ıes posible que yN sea lejos de y (tN ) . 15

2.1.3

Una cota para el error.

La soluci´ on y (t) de (1.1) satisface id´enticamente y (ti + h) = y (ti ) + hf (ti , y (ti )) − di

(2.3)

donde di = −y (ti + h) + y (ti ) + hf (ti , y (ti )) es el llamado error de truncamiento local, o cantidad por la cual y (t) no satisface la ecuaci´on de diferencias (2.2) localmente en ti . El error de truncamiento global εi = yi −y (ti ) es la diferencia entre la soluci´ on exacta de la ecuaci´ on de diferencias (2.2) y la soluci´ on correcta de la ecuaci´ on diferencial (1.1). Restando (2.3) de (2.2) se tiene εi+1 = εi + hf (ti , yi ) − hf (ti , y (ti )) + di .

(2.4)

Si el problema satisface las condiciones de un teorema de existencia y unicidad en una vecindad de yi , y (ti ), existe un constante de Lipschitz L , tal que f (ti , yi ) − f (ti , y (ti )) ≤ L yi − y (ti )

(2.5)

utilizando una norma vectorial si se tiene un sistema de ecuaciones, o un valor absoluto si se trata de una sola ecuaci´ on. Con esta condici´ on, es posible deducir de (2.4) las desigualdades εi+1  ≤ εi  + h f (ti , yi ) − f (ti , y (ti )) + di  ≤ εi  + hL yi − y (ti ) + di  ≤ εi  + hL εi  + di 

(2.6)

y si la condici´ on de Lipschitz se aplica uniformemente, con constante L, en una vecindad de la soluci´ on num´erica y de la soluci´on exacta para t en el rango t0 ≤ t ≤ tN , y D = max{di  : 0 ≤ i < N }, se tiene εi+1  ≤ (1 + hL) εi  + D ,

i = 0, 1, ..., N − 1.

(2.7)

on de diferencias Se demuestra f´ acilmente que, si ei es la soluci´on de la ecuaci´ ei+1 = (1 + hL) ei + D ,

i = 0, 1, ..., N − 1

(2.8)

con la condici´ on inicial e0 = ε0 , entonces εi  ≤ ei ,

i = 0, 1, ..., N

16

(2.9)

Efectivamente e0 ≥ ε0  ≥ 0, y se procede por inducci´ on para demostrar que ei ≥ εi  ≥ 0 para i > 0. Si es cierto para un i dado, entonces, restando (2.7) de (2.8). ei+1 − εi+1  ≥ (1 + hL) (ei − εi ) ≥ 0

(2.10)

y como, por definici´ on de norma εi+1  ≥ 0, se obtiene ei+1 ≥ εi+1  ≥ 0, lo que demuestra la proposici´on (2.9). Adem´as, la soluci´ on de (2.8) es muy f´ acil de obtener: ei = (1 + hL)i e0 +

 D

(1 + hL)i − 1 hL

(2.11)

de donde se obtiene la cota  D

(1 + hL)i − 1 . (2.12) hL El primer t´ermino en esta cota representa un error posible debido a ε0 , el error en los valores iniciales, que puede efectivamente existir, o por redondeo en su representaci´on, o por aproximaciones f´ısicas o num´eticas en el proceso que les genera. Para eliminar la dependencia de i, se nota que 1+ hL ≤ ehL , y se obtiene

εi  ≤ (1 + hL)i ε0  +

yi − y (ti ) ≤ eihL ε0  + =

D hL

e(ti −t0 )L ε0 



+



eihL − 1

D hL



e(ti −t0 )L



(2.13)

−1

y como esto es una funci´on creciente de ti  D (b−t0 )L e −1 (2.14) hL para t0 ≤ ti ≤ b. Para desplegar el comportamiento con h, se necesita todav´ıa una cota para D. Si f (t, y), adem´ as de satisfacer una condici´ on de Lipschitz en y, con constante L, satisface una condici´ on similar en t, con constante K, se obtiene

yi − y (ti ) ≤ e(b−t0 )L ε0  +

17

di  = y (ti + h) − y (ti ) − hf (ti , y (ti ))  ti +h

= ≤ ≤

ti

 ti+h ti

 ti +h ti



{f (t, y (t)) − f (ti , y (ti ))} dt

f (t, y (t)) − f (ti , y (ti )) dt {f (t, y (t)) − f (ti , y (t))

+ f (ti , y (t)) − f (ti , y (ti )) } dt ≤

 ti +h ti

(2.15)

{K (t − ti ) + L y (t) − y (ti )} dt  . t y

≤ h2 K/2 + hL maxti ≤t≤ti +h ≤ h2 K/2 + hL

ti



(τ ) dτ

 ti +h y. (t) dt ti .





≤ h2 K/2 + L maxti ≤t≤ti +h y (t) de donde .





D ≤ h2 K/2 + L maxt0 ≤t≤b y (t)

(2.16)

y sustituyendo en (2.14)



yi − y (ti ) ≤ e(b−t0 )L ε0  + h e(b−t0 )L − 1

.



K/2L + maxt0 ≤t≤b y (t) (2.17) on num´erica converge hacia la Esto demuestra que, si ε0 = 0, la soluci´ soluci´ on correcta cuando h tiende a cero. Es posible llegar a esta conclusi´ on sin suponer que f (t, y) satisface una condici´ on de Lipschitz en t, la continuidad siendo suficiente. Por otra parte, si se sabe que la soluci´ on y (t) tiene una segunda derivada cont´ınua (por ejemplo, si f (t, y)tiene derivadas cont´ınuas), es posible obtener una cota un poco mas fina. En efecto. 1 .. (2.18) di = − h2 y (ti + θh) para un θε [0, 1] 2 es una forma del resto de la serie de Taylor en ecuaci´on (2.3). En este caso

18

.. 1 2 h maxt0 ≤t≤b y (t) 2 que es una cota menor que la de ecuaci´on (2.16).

D = max0≤i 1/50 se tiene 1 − 100h < −1 y por consecuencia (1 − 100h)i crece en valor absoluta con i, y cambia de signo en cada paso. Por ejemplo, si h = 1/20. (−4)i − 1 2000 tomando ε0 = 0, y despu´es de 100 pesos |ε100 | >> 1050 . εi =

22

(3.6)

Por otra parte, para el problema .

y= 2t,

y (0) = 0

(3.7)

que tiene la misma soluci´ on i (3.8) 400 con el mismo paso h = 1/20 y ε0 = 0. En este caso ε100 = − 14 , un error de 1%. Se ve que la funci´ on complementaria con fuerte amortiguamiento y sin inter´es pr´actico est´a restringiendo el paso m´ aximo que puede utilizarse y requeriendo mucho m´ as trabajo que lo necesario para resolver con suficiente precisi´on el problema de inter´es. Sistemas de ecuaciones diferenciales con esta propiedad de ”stiffness” ocurren en el estudio de reacciones qu´ımicas, columnas de destilaci´ on, estabilidad de circuitos el´ectricos, y en la modelaci´ on de muchos otros sistemas f´ısicos, donde una de las cantidades involucrdas cambian mucho m´ as r´ apidamente que otras. Para el sistema . y = f (t, y) (3.9) εi = −ih2 = −

con soluci´ on y (t), es posible estudiar la variaci´ on local de la soluci´ on acerca de un valor fijo t = T mediante la sustituci´on 



y T + t = y(T ) + y  (t ) en cual caso

y  (t )

(3.10)

satisface

y  (t ) = f (T, y (T )) + .

= A (T ) y  (t ) + donde b (t ) = f (T, y (T )) +







∂f ∂t (T,y(T )) b (t )

· t + A (T ) y  (t )

(3.11)



∂f  ∂t (T,y(T )) · t , A (T ) es la matriz Jacobiana   ∂fi , y t queda peque˜ no para permitir elemento (i, j) igual a ∂y j (T,y(T )) 2 2 se desprecien t´erminos de orden t , y .

con que

Si A(T ) tiene n vectores propios linealmente independientes, estos forman una matriz U , no singular, tal que A (T ) = U ΛU −1 23

(3.12)

donde Λ es una matriz diagonal con elementos λi , i = 1, 2, ..., n, los valores propios de A (T ). Haciendo el cambio de variables  

 

Y t = U −1 y  t

(3.13)

la ecuaci´ on (3.11) se escribe .

 

 

 

y t = ΛY t + U −1 b t

(3.14)

un sistema de ecuaciones sin interconecciones, la i-´esima ecuaci´on conteon complementaria correspondiente niendo s´ olo el variable yi (t ). La funci´ λit
es e . Si unos de los λi tienen parte real grande negativo, el sistema (3.9) es ”stiff”, al menos localmente cerca de t = T . Como los λi son valores propios de la matriz real A(T ), unos o todos pueden tener valores complejos. La parte real representa crecimiento exponencial (si es positivo) o amortiguamiento (si es negativo) de la funci´ on complementaria, o modo, correspondiente. La parte imaginaria representa oscilaci´ on. Un concepto u ´til para m´etodos paso-por-paso, sobre todo si se aplican a ecuaciones ”stiff”, se define con respecto a la forma homog´enea de una de las ecuaciones del sistema (3.14), o sea .

y = λy

(3.15)

Para h real positivo, y λ en el plano complejo, hλ se encuentra en la regi´ on de estabilidad absoluta de un m´etodo si, cuando se aplica a la ecuaci´ on (3.15) con paso h, ninguna secuencia de valores iniciales yi , i = 0, 1, ..., K, genera una soluci´ on creciente en los pasos siguientes. Es claro que, para problemas ”stiff”, esta regi´on deber´ıa de incluir una parte grande de la mitad izquierda del plano complejo. Aplicando el m´etodo de Euler a la ecuaci´ on (3.15) se tiene yi+1 = yi + hλyi = (1 + hλ)yi

(3.16)

La soluci´ on no crece si |1 + hλ| < 1, o sea la regi´on de estabilidad absoluta para el m´etodo de Euler es el interior del c´ırculo con centro en −1 y radio 1.

24

.

3.1.3

Un m´ etodo con regi´ on de estabilidad absoluta m´ as grande.

El m´etodo de Euler por atr´ as se escribe yi+1 = yi + hf (ti+1 , yi+1 )

(3.17)

donde la derivada se evalua al final del paso. Aplic´ andolo a (3.15) se tiene yi+1 = yi + hλyi+1 =  

 

1 yi 1 − hλ

(3.18)

1 Si hλ tiene parte real neativo,  1−hλ  < 1 y toda la mitad izquierda del plano complejo se incluye en la regi´ on de estabilidad absoluta del m´etodo de Euler por atr´ as.

25

Aplic´ andolo al problema (3.1) se tiene 

yi+1 = yi + h 2ti+1 + 100t2i+1 − 100yi+1



(3.19)

y repitiendo la sustituci´ on yi = t2i + εi εi+1 =

h2 εi + 1 + 100h 1 + 100h

(3.20)

con soluci´ on ε0

h εi = i − 100 (1 + 100h)





1 (1 + 100h)i

−1 .

(3.21)

A´ un un error inicial ε0 bastante grande se amortigua despu´es de un rango normal de pasos. Con h = 1/20 y ε0 = 0 el error despu´es de 100 pasos es ε100

1 =− 2000

1 6100



−1 

1 2000

(3.22)

M´etodos que utilizan valores de la derivada al final del paso se llaman impl´ıcitos. El ejemplo dado, del m´etodo de Euler impl´ıcito, sugiere que el concepto es u ´til para el tratamiento de problemas ”stiff”. Para ecuaciones diferenciales lineales, la ecuaci´ on de diferencias que resulta puede resolverse directamente para yi+1 , pero en el caso no-lineal se necesita una t´ecnica iterativa. El m´etodo m´as sencillo para (3.17) es de (0) tomar una estimaci´on inicial yi+1 , por un predictor, como el m´etodo de Euler expl´ıcito, y de iterar por el proceso de sustituci´on sucesiva, con el corrector de Euler impl´ıcito. 

(j)

(j−1)

yi+1 = yi + hf ti+1 , yi+1



,

j = 1, 2, ...

(3.23)

(j)

Si se define δ(j) por yi+1 = yi+1 − δ(j) entonces 

yi+1 + δ(j) = yi + hf ti+1 , yi+1 + δ(j−1)



(3.24)

y despreciando t´erminos de orden δ2 δ(j) = hA (t) δ(j−1)

(3.25)

lo que sugiere que (3.23) s´ olo convergir´ a para problemas ”stiff” si h es muy peque˜ no, 1/h mayor que el valor absoluto m´ aximo de los valores propios

26

de A(t), una restricci´on m´ as fuerte que para Euler expl´ıcito. Esto se ve claramente con su aplicaci´on a (3.15) cuando 



(0)

δ(j) = hλδ(j−1) = yi+1 − yi+1 (hλ)j

(3.26)

lo que tiende a cero cuando j aumenta s´olo si |hλ| < 1. Otro m´etodo iterativo es el de Newton-Raphson cuya aplicaci´on a (3.17) puede derivarse de esta manera. Definiendo δ(j) de la misma manera 

(j)

(j)

yi+1 − δ(j) = yi + hf ti+1 , yi+1 − δ(j) 

(j)





(3.27)

 yi + hf ti+1 , yi+1 − hA(j) (t) δ(j)

despreciando erminos en δ2 , y con A(j) (t) una matriz con elemento (i, j)   t´ ∂fi (j) igual a ∂yj (j) . Resolviendo para δ y=yi+1





(j)



(j)

I − hA(j) (t) δ(j) = yi+1 − yi − hf ti+1 , yi+1



(3.28)

un sistema de n ecuaciones lineales en el caso de un sistema de ecuaciones diferenciales. Por la aproximaci´ on en (3.27) no se tiene el valor correcto (j+1) (j) (j) para δ , y se toma yi+1 = yi+1 − δ(j) . Normalmente, la convergencia de este m´etodo es buena con un predictor adecuado. En caso de problemas, puede reducirse el paso h. La aplicaci´ on del m´etodo de Newton-Raphson a una ecuaci´ on diferencial lineal, como (3.15), da convergencia en una sola iteraci´ on independiente del (0) valor de yi+1 . Aunque tiene mejor rapidez de convergencia, el m´etodo de NewtonRaphson puede parecer costoso, debido al c´ alculo de las derivadas parciales ∂fi 2 elementos), y la resoluci´ ( n o n del sistema lineal en cada iteraci´on. Sin ∂yj embargo, para un sistema ”stiff”, permite pasos mucho m´as grandes que otros m´etodos, lo que puede compensar completamente para el trabajo adicional por paso. Adem´ as, muchas veces las derivadas parciales no var´ıan r´ apidamente con y, y no es necesario reevaluarlos, y reinvertir la matriz on (si hay m´ as de una), ni tampoco en cada I − hA(j) (t), en cada iteraci´ paso en t.

27

Chapter 4

METODOS CON UN SOLO PASO 4.1 4.1.1

Clasificaci´ on de m´ etodos M´ etodos con un solo paso y con pasos m´ ultiples.

Los dos m´etodos — Euler expl´ıcito y Euler impl´ıcito — vistos hasta aqu´ı, son m´etodos con un solo paso, o sea, la ecuaci´ on de diferencias al i-´esimo paso involucra valores de la soluci´ on aproximada y de su derivada s´ olo en el intervalo de ti−1 a ti . Para una clase muy extensa de m´etodos paso-por-paso para ecuaciones diferenciales del primer orden, un error de truncamiento  r+1 genera un error propagado o global de orden 0 (hr ) local de orden 0 h y se dice que el m´etodo es de orden r. Para los dos m´etodos de Euler r = 1. En general, puede esperarse que mientras m´ as alto sea el orden del m´etodo, la soluci´ on que genera ser´ a m´as precisa para un tama˜ no dado del paso h. Para m´etodos de orden mayor, se necesita m´as informaci´ on para construir la soluci´on aproximada en un intervalo h, y esto puede conseguirse por una de dos maneras: (i) dentro del intervalo de ti−1 a ti (m´etodos con un solo paso, o mejor dicho, con pasos separados); (ii) de pasos anteriores (m´etodos con pasos m´ ultiples, o mejor dicho, con pasos ligados). En esta secci´on se tratar´ an m´etodos con un solo paso, m´as correctamente considerados como m´etodos que tratan cada paso como un nuevo problema con condiciones iniciales, sin tener memoria de valores anteriores

28

de la soluci´ on o de sus derivadas.

4.1.2

El m´ etodo de la serie de Taylor.

El m´etodo de Euler (expl´ıcito) utiliza los primeros dos t´erminos de la serie de Taylor en ti para aproximar y (ti+1 ). Si se tuvieran r derivadas, se podr´ıa poner h2 .. hr y i +... + yi(r) (4.1) 2 r! y el primer t´ermino del error de truncamiento local resultante estar´ıa dado (r+1) hr+1 yi , o sea, se tendr´ıa un m´etodo de orden r. Si f (t, y) tiene por − (r+1)! .

yi+1 = yi + h yi +

suficientes derivadas, abreviando f (t, y) por f , ∂f ∂t por ft , segunda y tercera derivadas de y pueden obtenerse por

∂f ∂y

por fy etc., las

. d (f (t, y)) = ft + fy y = ft + fy f dt

(4.2)

y = ftt + 2fty f + fyy f 2 + fy ft + fy2 f

(4.3)

..

y= ...

Es claro que la complejidad de estas expresiones crece r´apidamente, y que el m´etodo de la serie de Taylor resultar´ıa pr´ actico s´olo cuando f (t, y) tiene una forma bastante sencilla. Por ejemplo, no siempre se tiene f en forma anal´ıtica, sino como resultado de otro proceso num´erico, tal como la optimizaci´on. Sin embargo, cuando es posible, por ejemplo cuando f es un polinomio de grado reducido, el m´etodo tiene la ventaja que da directamente una estimaci´on del error de truncamiento local, lo que puede controlarse variando el n´ umero de t´erminos que se utilizan, adem´ as del paso h para lo cual no se necesita reevaluar las derivadas. Se consider´an despu´es las t´ecnicas para controlar el error.

4.1.3

M´ etodos de extrapolaci´ on.

Como el error global baja cuando el paso h se reduce, un m´etodo com´ un para averiguar la precisi´ on de una soluci´ on es de repetir la integraci´on, reduciendo el paso a la mitad, y de comparar los dos resultados. Si se utiliza el m´etodo de Euler, y se obtienen las aproximaciones yN1 , yN2 para y (b), entonces, utilizando los m´etodos de la secci´on 2.1.5 se tiene

29



yN1 = y (b) + hδ (b) + 0 h2 yN2 = y (b) +



(4.4)

 h δ (b) + 0 h2 2

(4.5)

Eliminando δ (b) 

y (b) = 2yN2 − yN1 + 0 h2 Ejemplo 4.1.1

.

y = −y ,



(4.6)

y (0) = 1

Aplicando el m´etodo de Euler con varios valores del paso h, para determinar y(1), y eliminando el error principal como en (4.6). Se obtienen los resultados siguientes: h 1 1/2 1/4 1/8 1/16 1/32 1/64

y1/h por Euler 0.0 0.25 0.316406 0.343609 0.356074 0.362055 0.364986

error −0.36788 −0.11788 −0.05147 −0.02427 −0.01182 −0.00582 −0.00289

2y2/h − y1/h 0.5 0.3828125 0.3708116 0.3685393 0.3680364 0.3679177

error 0.1321206 0.0149331 0.9929321 0.0006599 0.0001569 0.0000382

El error de la aproximaci´ on 2y1/2h − y1/h baja r´ apidamente cuando h tiende a cero, debido a su comportamiento como h2 . El error del valor obtenido de los casos h = 1, h = 1/2 es mayor que el error del m´ etodo de   Euler con h = 1/2. Esto se debe a que los t´erminos de orden 0 h2 son todav´ıa demasiado grandes para despreciarse como en (4.4). Esto es un ejemplo de la t´ecnica de extrapolaci´ on de Richardson hacia h = 0, y en este caso da un m´etodo de segundo orden. En la pr´actica, muchas veces la soluci´ on se necesita en puntos dentro del intervalo global de integraci´ on y pueden llevarse a cabo las dos integraciones juntas, y hacer la extrapolaci´ on (4.6) en cada paso de la integraci´ on con pasos m´as grandes. Sin embargo, parecer´ıa mejor, en lugar de guardar las dos secuencias independientes de yi , de reemplazar los dos valores por el valor mejorado despu´es de cada extrapolaci´on. Para el paso m´ as grande 30

y1/h (ti + h) = yi + hf (ti , yi )

(4.7)

y para los pasos chicos qi = y2/h (ti + h/2) = yi +

h f (ti , yi ) 2



(4.8)



h h (4.9) f ti + , q i 2 2 El valor extrapolado, que debe tomarse como yi+1 para los dos tama˜ nos de paso siguiente, es y2/h (ti + h) = qi +

yi+1 = 2y2/h (ti + h) − yh (ti + h) = 2yi + hf (ti , yi ) + hf ti + h2 , qi − yi − hf (ti , yi ) 

= yi + hf ti +

h 2 , qi



(4.10)

Los c´alculos que deben hacerse en la pr´ actica se dan en (4.8) y (4.10), y constituyen el m´etodo del punto medio. Para encontrar el error de truncamiento local de este m´etodo, se sustituye la soluci´on exacta en estas dos ecuaciones, dando lugar a qi (ti ) = y (ti ) + h2 f (ti , y (ti )) . = y (ti ) + h2 y (ti ) 

.

(4.11)

y (ti+1 ) −y (ti ) − hf ti + h2 , y (ti ) + h2 y (ti )   . 2 .. = h y + h2 y +0 h3   . 2 2 −hf (ti , y (ti )) − h2 ∂f − h2 A (t) y +O h3 ∂t   .. . y = O h3 pues y = ∂f ∂t + A (t)

(4.12)

El error de truncamiento local es de un orden mayor en h con una evaluaci´ on m´as de la funci´ on f . Otra forma del m´etodo del putno medio se obtiene reemplazando h por on en el punto medio. 2h, y notando que qi es una aproximaci´on para la soluci´ Se reescribe (4.10) como yi+1 = yi−1 + 2hf (ti , yi )

(4.13)

que es un m´etodo de dos pasos con el mismo nombre. El valor de yi+1 de (4.13) se utiliza para evaluar la derivada que se necesita junto con yi para determinar yi+2 , por medio de (4.13), reemplazando i por i + 1. 31

Este m´etodo de pasos m´ ultiples se ha introducido en este momento porque se utiliza como base de un m´etodo de exrapolaci´ on de Bulirsch y Stoer. Se trabaja sobre un paso b´ asico H, normalmente bastante grande. En cada refinamiento s del m´etodo, se divide H en un n´ umero par Ns de subintervalos, por ejemplo {Ns } = {2, 4, 6, 8, 12, 16, 24, ....} con tama˜ no hs = H/Ns , y se aplica el m´etodo del punto medio de dos pasos con h = hs . Esto necesita un m´etodo de un solo paso para arrancar, y se utiliza el m´etodo de Euler yi = y0 + hs f (t0 , y0 ) .

(4.14)

Puede demostrarse que el m´etodo (4.13), aunque convergente y estable cuando h tiende a cero, tiende una regi´ on de estabilidad absoluta constituido nada m´ as por el punto hλ = 0, y no es estable para hλ
= 0. Para amortiguar el crecimiento, se modifica el valor yNs obtenido, promediando con el resultado de un paso del tipo Euler por atr´ as para dar  1 yNs + yNs−1 + hs f (tNs , yNs ) (4.15) 2 El resultado es un m´etodo que tiene una regi´ on de estabilidad absoluta finita, incluyendo el intervalo del eje real negativo −3 ≤ Hλ < 0 para Ns = 2, y on creciendo para valores mayores de Ns . Gragg ha demostrado que la funci´ T (h) tiene una expansi´ on asint´ otica en potencias pares de h (s)

T0

= T (hs ) =

T (h) = y (t0 + H) + α1 h2 + α2 h4 + ...

(4.16)

si f (t, y) tiene suficientes derivadas. La idea general de extrapolaci´ on hacia ∧

h = 0 es de ajustar una funci´ on T (h) a los valores T (h) que se tienen, y ∧

de evaluar T (0) como una aproximaci´on mejor al valor buscado. Bulirsch y Stoer, adem´ as de estudiar extrapolaci´ on polinomial como la de Richardson, han utilizado funciones racionales, que en muchos casos es mejor. Definen ∧ (i)

T m (h) =

(i)

(i)

(i)

(i)

(i)

(i)

p0 + p1 h2 + ... + pµ h2µ q0 + q1 h2 + ... + qν h2ν

(4.17)

con µ = m/2 si m es par, o (m − 1) /2 si m es impar, ν = m − µ, y (i)

∧ (i)

(i)

con los coeficientes pj , qj determinados tales que T m (hs ) = T (hs ) para s = i, i + 1, i + 2, ..., i + m. Resulta que no es necesario determinar estos

32

(i)

∧ (i)

coeficientes, pues es posible determinar los valores extrapolados Tm =T m (0) por las f´ ormulas de recurrencia siguientes (i)

T−1 = 0 para arrancar (i) = T (hi ) la aproximaci´ on constante con h = hi T0 (i+1) (i) Tm−1 − Tm−1 (i) (i+1) para M ≥ 1 Tm = Tm−1 +    2 (i+1) (i) Tm−1 −Tm−1 hi 1 − (i+1) (i+1) − 1 hi+m Tm−1 −Tm−2

Esto es diferente de las f´ ormulas para extrapolaci´on polinomial por la introducci´ on del t´ermino en par´entesis cuadrados. Permite de formar la tabla siguiente (0)

T0 T−1 (1)

(1)

(0)

T1

T0 T−1 (2)

(2)

T0 .. .

.. .

(0)

T2

(1)

T1

(1)

T−1 (i)

(i)

(0)

T2 .. .

.. .

(i−1)

T0

(0)

T3

T4

(i−2)

(i−1)

T2

T1

T0 Puede demostrarse que



(i) − T (0) = 0 h2i h2i+1 ...h2i+m Tm



(4.18)

Las f´ ormulas de recurrencia que se obtienen en procesos de extrapolaci´on son muy sensitivos a posibles redondeos, y deben organizarse de manera cuidadosa, con un l´ımite sobre el valor de m, que depende de la computadora y de la precisi´ on requerida. Bulirsch y Stoer recomendaron m ≤ 6. Si  (0)   (0) T6 − Ts  no satisface el criterio de error, abandonan el primer giagonal, y siguen por la columna m = 6 hasta i = 10. Si todav´ıa no hay convergencia, repiten todo el c´ alculo con un valor reducido de H, que de todos mos escogen con base en (4.18), prediciendolo para satisfacer el criterio de error en (0) T6 . El m´etodo, que se considera un m´etodo de un solo paso (H grande), normalmente es muy bueno para problemas sin ”stiffness”, sobre todo si 33

la evaluaci´ on de f (t, y) no necesita m´as de 25 operaciones artim´eticas por componente. De otra manera, puede resultar costoso en evaluaciones de f .

4.1.4

M´ etodos del tipo Runge-Kutta.

N´ otese ahora que el m´etodo del punto medio (4.8) y (4.10) genera t´erminos que coinciden con la serie de Taylor truncada h2 .. yi (4.19) 2 como se vi´o en (4.12). Puede buscarse una clase de m´etodos m´as general que coinciden con la serie de Taylor truncada en un n´ umero determinado de t´erminos. Estos son los m´etodos de Runge-Kutta, que se derivan en el cap´ıtulo siguiente. .

yi+1 = yi + h y i +

34

Chapter 5

METODOS DE RUNGE-KUTTA Y CONTROL DEL ERROR. 5.1 5.1.1

M´ etodos del tipo Runge-Kutta. Forma aut´ onoma, y notaci´ on.

El sistema de ecuaciones (1.4) con condiciones iniciales y (0) = y0 puede, en muchas circunstancias, simplificarse por la adici´ on de una ecuaci´ on m´ as: y0 = 1, .

y 0 (0) = 0

(5.1)

nade al principio del vector y, ahora de n + 1 para el elemento y 0 , que se a˜ 0 nade de la misma manera f 0 = 1 al componentes. Como y (t) = t, si se a˜ principio del vector f , el problema puede escribirse en forma aut´ onoma, sin que t aparezca expl´ıcitamente: .

y = f (y) ,

(5.2)

y (0) = y0

Esto permite definir las siguientes notaciones: f i es el elemento i-´esimo de f , con f 0 = 1; fj es la derivada

∂f ∂y j



evaluada en ts

35

f0 =

∂f ∂t



.

5.1.2

M´ etodo del punto medio por extrapolaci´ on.

El m´etodo de Euler puede escribirse 

ys+1 − ys = hf (ys ) = y (ts+1 ) − (ts ) − ch2 + 0 h3



(5.3)

mientras sobre dos pasos de tama˜ no h/2 se escribe 

ys+ 1 2

2

 1 1 h − ys = h (ys ) = y ts+ 1 − y (ts ) − c 2 2 2



ys+1 − ys+ 1 2



+ 0 h3 2

 1  1 h = hf ys+ 1 = y (ts+1 ) − y ts+ 1 − c 2 2 2 2



(5.4) 

+ 0 h3



(5.5)

o sea 



ys+1 − ys = 12 h f (ys ) + f ys + 12 hf (ys ) = y (ts+1 ) − y (ts ) − 2c



2

1 2h

 



+ 0 h3 .

(5.6)

Aplicando la extrapolaci´ on de Richardson a este paso, como en ecuaci´on (4.10), o sea eliminando c entre (5.3) y (5.6), se tiene otra vez el m´etodo del punto medio 



 1 ys+1 − ys = hf ys + hf (ys ) = y (ts+1 ) − y (ts ) + 0 h3 2

(5.7)

un caso particular de un m´etodo del tipo Runge-Kutta de orden dos.

5.1.3

M´ etodos del tipo Runge-Kutta de orden dos.

Un m´etodo general del tipo Runge-Kutta, con dos evaluaciones de la funci´ on f , se escribir´ıa = hf (ys ) k0 = hf (ys + αk0 ) k1 ys+1 −ys = βk0 + γk1

(5.8)

donde α, β, γ deben ajustarse para que la expresi´ on para ys+1 − ys coincida hasta donde se puede (o se quiere) con la serie de Taylor acerca de ys , o sea 

3 y (ts+1 ) − y (ts ) = hf (ys ) + 12 h2 df dt + 0 h

= hf (ys ) + 36



1 2 n i i=0 fi f 2h

  + 0 h3 .

(5.9)

Para el m´etodo num´erico k1 = hf (ys ) + αh2



n 



fi f i + 0 h3 ,

(5.10)

i=0

ys+1 − ys = (β + γ) hf (ys ) + γαh2

n 





fi f i + 0 h3 ,

(5.11)

i=0

y el m´etodo es de orden dos si β = 1 − γ,

α = 1/2γ.

(5.12)

Si γ = 1, esto es el m´etodo del punto medio 

k0 = hf (ys ) ,



k1 = hf ys + 12 k0 ,

ys+1 − ys = k1 .

(5.13)

Si γ = 12 , se tiene un m´etodo modificado del trapecio k0 = hf (ys ) ,

ys+1 − ys =

k1 = hf (ys + k0 ) ,

1 2

(k0 + k1 ) . (5.14)

Si γ = 34 , se tiene uno de los m´etodos de Heun 

k0 = hf (ys ) ,



k1 = hf ys + 23 k0 ,

ys+1 − ys =

1 4

(k0 + 3k1 ) , (5.15) una f´ ormula de cuadratura del tipo Radau propuesta para minimizar una forma particular de la cota para el error de truncamiento local. La nomenclatura de estos m´etodos no es est´andar, y var´ıa mucho entre los diferentes textos.

5.1.4

La clase de m´ etodos de orden cuatro.

Los m´etodos del tipo Runge-Kutta, que se utilizan m´ as com´ unmente, son los del cuarto orden, que toman la forma general siguiente: k0 k1 k2 k3 ys+1

= hf (ys ) = hf (ys + α10 k0 ) = hf (ys + α20 k0 + α21 k1 ) = hf (ys + α30 k0 + α31 k1 + α32 k2 ) −ys = α40 k0 + α41 k1 + α42 k2 + α43 k3 37

(5.16)

mientras la serie de Taylor para y (ts+1 ) genera los t´erminos 

y (ts+1 ) − y (ts ) = hf (ys ) + 12 h2 ni=0 fi f i    + 16 h3 ni=0 nj=0 fij f i f j + fi fji f j 1 4 + 24 h

n

i=0

n

j=0

n

k=0



fijk f i f j f k + 3fij fki f j f k



i f j f k + f f i f j f k + 0 h5 +fi fjk i j k (5.17) Haciendo la expansi´ on de k1 , k2 , k3 en la misma forma, y comparando los coeficientes de cada t´ermino, el m´etodo es de orden cuatro si

α40 + α41 α41 β1 α41 β12 α41 β13

+ + + +

α42 α42 β2 α42 β22 α42 β23 α42 α21 β1 α42 β2 α21 β1 α42 α21 β12

+ + + + + + +

α43 α43 β3 α43 β32 α43 β33 α43 (α31 β1 + α32 β2 ) α43 β3 (α31 β1 + α32 β2 ) α43 α31 β12 + α32 β22 α43 α32 α21 β1

=1 = 1/2 = 1/3 = 1/4 = 1/6 = 1/8 = 1/12 = 1/24 (5.18)

donde β1 = α10 β2 = α20 + α21 β3 = α30 + α31 + α32

(5.19)

Estas ecuaciones tienen una infinidad de soluciones, de las cuales el m´etodo cl´asico de Runge-Kutta utiliza la f´ ormula de cuadratura de Simpson de la manera siguiente k0 k1 k2 k3 ys+1

= hf (ys ) = hf (ys + 12 k0 ) = hf (ys + 12 k1 ) = hf (ys + k2 ) −ys = 16 (k0 + 2k1 + 2k2 + k3 )

(5.20)

Para un m´etodo del tipo Runge-Kutta de orden ≤ 4, el n´ umero de evaluaciones de la funci´ on f es igual al orden, pero para orden 5 o´ 6, se necesita una evaluaci´ on m´ as, mientras para orden > 6 se necesitan al menos dos

38



evaluaciones suplementarias. Estos m´etodos difieren del m´etodo de la serie de Taylor en que no es f´ acil aumentar el orden durante una integraci´ on si los errores de truncamiento son demasiado grandes, como para esto se necesita un conjunto de multiplicadores totalmente nuevo. Hasta unas propuestas muy recientes, ha sido normal utilizar un m´etodo de cuarto orden todo el tiempo, y de ajustar el error de truncamiento mediante el paso h. Sin embargo, se necesita un criterio para estos ajustes. Sea el m´etodo (5.21) ys+1 − ys = hψ (ys , h) con error de truncameinto local ds (h) = hψ (y (ts ) , h) − (y (ts+1 ) − y (ts ))  6 5 = h φ (y) + 0 h

(5.22)

donde se ha especializado al caso de un m´etodo de cuarto orden. Sea el error de truncamiento global 

s = ys − y (ts ) = h4 δ (ts ) + 0 h5



(5.23)

Para el problema (5.2) la funci´ on δ (t) satisface .

δ = A (t) δ + φ (y)

(5.24)

ecuaci´ on an´ aloga al (2.27), con A (t) la matriz Jacobiana de f . Cuando el paso no es constante, es posible utilizar una funci´ on θ (t) para describir su variaci´ on, mediante la definici´ on hs = hθ (ts ) ,

ts+1 = ts + hs

(5.25)

=1

(5.26)

donde se requiere
b dt 0

θ (t)

para que el n´ umero total de pasos en el intervalo de 0 a b sea aproximadamente 1/h, y el trabajo hecho dependa solamente de h, y no de θ. Entonces, la ecuaci´ on (5.24) se modifica por la introducci´ on de θ (ts ) en (5.22), y queda de la forma .

4 δ = A(t)δ + θ (t) φ (y)

de donde δ (t) puede escribirse como 39

(5.27)


t

δ (t) =

0

G (v, t) φ (y (v)) θ 4 (v) dv

(5.28)

donde G (v, t) es la matriz de funciones complementarias, que satisface la ecuaci´ on diferencial homog´enea, y la condici´ on inicial .

G= A (t) G,

G (v, v) = I,

(5.29)

y que representa el efecto propagado en t de un error unitario cometido en v. Deben notarse las siguientes propiedades de G (v, t): G (v, t) G (u, v) tambi´en satisface la ecuaci´on diferencial homog´enea, con condici´ on inicial G (v, v) G (u, v) = G (u, v) ; pero la soluci´ on u ´nica de este problema es G (u, t), de donde G (v, t) G (u, v) = G (u, t) La mejor funci´ on θ (t) es la que minimiza δ (b). El m´ınimo de
b 0

G (v, b) φ (y (v)) θ 4 (v) dv

(5.30)

que no es exactamente lo mismo, sino es lo mejor posible, ocurre cuando θ (t) = G (t, b) φ (y (t))

−1/5


b 0

G (v, b) φ (y (v))1/5 dv.

(5.31)

La contribuci´ on del paso de ts a ts+1 al error de truncamiento global en b es 5
t 
b s+1 1/5 4 4 G (t, b) φ (y (t)) θ (t) h dt G (v, b) φ (y (v)) dv  h 0

ts

(5.32) para la funci´ on θ (t) de la ecuaci´ on (5.31). Resulta que la mejor distribuci´ on de tama˜ nos de paso es tal que cada paso hace la misma contribuci´ on al error de truncamiento global en b, independientemente de su tama˜ no.

40

5.1.5

Suposiciones que permiten el ajuste pr´ actico del paso.

Este resultado te´orico es interesante, pero normalmente no es pr´actico aplicarlo. Generalmente, no se sabe el efecto propagado de un error de truncamiento local, representado por la funci´ on matricial G(t, b), que puede obtenerse s´olo por la integraci´ on de la ecuaci´ on (5.29) o algo similar, de antemano. En la ausencia de un conocimiento de G(t, b), una de las dos suposiciones siguientes parece razonable: (a) G(t, b nunca es mucho mayor que uno. Esto es el caso cuando la soluci´ on general (funci´ on complementaria) de la ecuaci´on diferencial (5.29) tiene solamente componentes oscilantes o con decaimiento, o sea cuando el problema original tiene estabilidad matem´atica. Una estrategia pr´ actica ser´ıa de ajustar el tama˜ no de cada paso para mantener la magnitud del error de truncamiento local aproximadamente constante. Unas implementaciones cambian el tama˜ no del paso solamente por factores de 2, pero se reconoce generalmente ahora que puede ajustarse en cada paso para que h5 φ (y) sea un poco menos que un cierto valor constante. (b) G (v, t) tiene componentes con crecimiento exponencial de una forma similar a la de y(t). Debe notarse que G(v, t) tiene una propiedad multiplicativa G (t, b) G (v, t) = G (v, b)

(5.33)

que es independiente de t, as´ı que para guardar G (t, b) φ (y) h5

(5.34)

aproximadamente constante, se debe guardar h5 φ (y) aproximadamente proporcional a G (v, t), o sea por la suposici´on (b), proporcional a y (t). En este caso, es el error de truncamiento local relativo a la soluci´on que debe mantenerse constante. Estas estrategias tiene la desventaja que no se tiene idea de antemano del n´ umero de pasos que van a utilizarse, cuando se especifica el error permisible por paso. Si quiere obtenerse una cierta precisi´on en la soluci´on en b, otra estrategia posible es de mantener la contribuci´ on al error por paso proporcional al tama˜ no h del paso. Con la suposici´ on (a) se ajusta h para que bh4 φ (y) sea un poco menos que la cota impuesta para el error en b. Con la suposici´ on (b) se debe mantener bh4 φ (y) / y un poco menos que la cota para el error relativo en b.

41

5.1.6

Estimaci´ on del error de truncamiento local.

En todos los casos, se necesita una estimaci´on para el error de truncamiento  6 5 local ds (h) = h φ (y) + 0 h , que no siempre se tiene f´acilmente. En el caso del m´etodo de la serie de Taylor, cuando pueden calcularse las derivadas elevadas de y (t), el t´ermino siguiente de la serie es la estimaci´on apropiada. Sin embargo, para los m´etodos del tipo Runge-Kutta, no hay una estimaci´ on tan evidente. El australiano, Merson, propuso un m´etodo con cinco evaluaciones de la funci´ on f = hf  (ys ) = hf ys + 13 k0

k0 k1









= hf ys + 16 k0 + 16 k1

k2

= hf ys + 18 k0 + 38 k2

k3

(5.35)



k4 = hf ys + 12 k0 − 32 k2 + 2k3 ys+1 −ys = 16 (k0 + 4k3 + k4 )



1 (2k0 − 9k2 + 8k3 − k4 ) era la estimaci´on propuesta para el error. donde 30 El error de truncamiento local de este m´etodo puede escribirse

n 1 5 n n n i=0 j=0 k=0 l=0 120 h



1 i j k l 24 fijkl f f f f i f jf kf l + 16 fij fkl i f j f kf l − 12 fi fjk l

+ 14 fijk fli f j f k f l + 18 fij fki flj f k f l i f j f kf l + 16 fij fki f j flk f l − 16 fi fjkl  j k l − 16 fi fji fkl f f

mientras la expansi´on de la estimaci´on de Merson es 1 4 h − 120

1 5 − 120 h

n

i=0

n

n

n

n

n

n

i=0 j=0

j=0 k=0

k=0 l=0





2 i f jf k f f i f j f k + 23 fij fki f j fk + 13 fi fjk  9 ijk 11 11 11 i j k l i j k l i j k l 108 fijkl f f f f + 18 fijk fl f f f + 36 fij fk fl f f i f j f kf l + 1 f f i f j f kf l + 5 f f i f j f kf l + 49 fij fkl l  3 ij k 54 i jkl   i f j f k f l + 1 f f i f j f k f l + 0 h6 + 16 fi fjk l 6 i j k l

(5.37) que claramente no es una buena estimaci´on, siendo en general de un orden menor en h que el error estimado. Coincide exactamente cuando f (y) es lineal en todos sus argumentos, cuando todas las derivadas de f de mayor orden que el primero desaparecen, y dejan el u ´nico t´ermino con el mismo 42



− 16 fi fji fkj flk f l + 0 h6 (5.36)



coeficiente en las dos expresiones. Ser´ıa injusto hacia Merson de decir que propuso su estimaci´on como v´alido para ecuaciones no-lineales, pero se ha empleado mucho como un m´etodo de aplicaci´ on general. Normalmente su estimaci´on es mayor que el error de truncamiento local, y algoritmos que ajustan los pasos en base en su estimaci´ on toman pasos m´as peque˜ nos que lo necesario. Aunque todav´ıa no se conocen muy ampliamente, existen otros m´etodos del mismo tipo general que estiman correctamente los t´erminos de orden h5 . Se da aqu´ıun ejemplo que, al costo de una evaluaci´ on adicional de la funci´ on f , da una estimaci´on del error sobre dos pasos consecutivos iguales. k0 k1 k2 k3 ys+1 k8

= hf (ys ) = hf (ys + 12 k0 ) = hf (ys + 14 k0 + 14 k1 ) = hf (ys − k1 + 2k2 ) −ys = 16 (k0 + 4k2 + k3 ) = hf ys − 16 k0 − 16k1 +

k4 = hf (ys+1 ) k5 = hf (ys+1 + 12 k4 ) k6 = hf (ys+1 + 14 k4 + 14 k5 )

− 121 6 k3 + 24k4 + k5 − 2k6 (5.38) donde el error de truncamiento local en ts+2 sobre los dos pasos es 1 90

46 3 k2





(k0 − 4k2 − 17k3 + 23k4 + 4k6 + k8 ) + 0 h6 .

Otra t´ecnica, de aplicaci´on m´ as general, para la estimaci´on del error sobre dos pasos para un m´etodo de Runge-Kutta, se basa en la extrapolaci´ on, o sea se repite la integraci´ on con un paso dos veces m´as grande, al costo de tres evaluaciones adicionales, como f (ys ) ya se ha evaluado. (Puede decirse que el m´etodo de Merson cuesta dos evaluaciones adicionales para dos pasos). Sobre dos pasos de tama˜ no h, el error de truncamiento local es G(ts+1 , ts+2 )ds (h) + ds+1 (h) = 2h5 φ(y) + 0(h6 )

(5.39)

mientras sobre un paso de tama˜ no 2h es ds (2h) = (2h)5 φ(y) + 0(h6 ).

(5.40)

La diferencia entre las soluciones es tambi´en la diferencia entre estes errores de truncamiento, y es 30h5 φ(y) + 0(h6 ) = 15 [G(ts+1 , ts+2 )ds (h) + ds+1 (h)] , 43

(5.41)

15 veces el error de truncamiento local en ts+2 sobre los dos pasos. Un aspecto de los m´etodos de Runge-Kutta que no se ha considerado es su estabilidad para tama˜ nos finitos del paso, o sea, sus regiones de estabilidad absoluta. Se considera entonces su aplicaci´ on a la ecuaci´ on .

y = λy.

(5.42)

Para cualquier m´etodo de cuatro evaluaciones (5.16) que satisface las condiciones (5.18 ) y (5.19), k0 k1 k2 k3 ys+1

= hλys = hλ (ys + α10 hλys ) = hλ (1 + hλβ1 ) ys = hλ (ys + α20 hλys + α21 hλ (1 + hλβ1 ) ys ) = hλ (1 + hλ (β2 + hλα21 β1 )) ys = hλ (ys + α30 hλys + α31 hλ (1 + hλβ1 ) ys + α32 hλ (1 + hλ (β2 + hλα21 β1 ))) ys = hλ (1 + hλ (β3 + hλ (α31 β1 + α32 β2 + hλα32 α21 β1 ))) ys = ys + α40 hλys + α41 hλ (1 + hλβ1 ) ys + α42 hλ (1 + hλ (β2 + hλα21 β1 )) ys +α43 hλ α32 β2 + hλα32 α21 β1 ))) ys  (1 + hλ  (β3 + hλ  (α31 β1 + 1 1 1 = ys + hλ 1 + hλ 2 + hλ 6 + 24 hλ ys 

= ys 1 + hλ + 12 h2 λ2 + 16 h3 λ3 +

1 4 4 24 h λ



que hab´ıa que esperarse como es una aproximaci´on del cuarto orden a la on de estabilidad es la parte del plano complejo para soluci´ on ys ehλ . La regi´ la cual     1 + hλ + 1 h2 λ2 + 1 h3 λ3 + 1 h4 λ4  < 1   2 6 24

(5.43)

y se encuentra casi completamente entre e (hλ) = −2.7 y e (hλ) = 0 y entre Im (hλ) = ±2.5, as´ı que m´etodos del tipo Runge-Kutta siempre convergen cuando h tiende a cero, pero son muy poco eficaces para problemas ”stiff”. La regi´on de estabilidad absoluta es un poco m´as grande para m´etodos de orden tres y cuatro que para m´etodos de ´orden uno y dos, pero como utilizan aproximaciones polinomiales para ehλ , no puede esperarse que se mejoren mucho para o´rdenes mayores.

44

Chapter 6

METODOS CON PASOS MULTIPLES 6.1 6.1.1

Introducci´ on a los m´ etodos con pasos m´ ultiples. Los m´ etodos de Adams.

Los m´etodos que se describen en seguida utilizan informaci´on acumulada sobre la soluci´ on, que no puede obtenerse durante un solo paso de integraci´ on. Por esto, es m´as f´ acil entenderlos si, para el momento, no se consideran los problemas particulares de arranque, sino que se supone que un n´ umero suficiente de pasos ya se han generado de alguna manera. Se concentra la atenci´on en el m´etodo para avanzar la soluci´ on por otro paso. Para fijar las ideas, sup´ ongase que ya se ha alcanzado el punto ts = sh, y que se tiene informaci´ on suficiente para definir un vector de polinomios zs (t − ts ), de lo cual cada componente es un polinomio de grado K, o sea zsi (τ ) =

K 

ζsik τ k

(6.1)

k=0

y que aproxima la soluci´on mente. Hay varias maneras de la soluci´ on que se tienen. yr , r = s − p, s − p + 1, ..., s

del sistema de ecuaciones diferenciales, localde relacionar estos polonomios con los valores Por ejemplo, pueden haberse guardado valores de la soluci´on calculada, y/o de sus derivadas .

y r = f (tr , yr )

45

(6.2)

Si se tienen K + 1 valores, estos definen de manera u ´nica el vector de polion com´ un es la adoptada en los nomios zs (τ ) que los interpola. Una selecci´ . m´etodos del tipo Adams, para los cuales la informaci´ on guardada es ys , y s , . . . . . y s−1 , ..., y s−k+1 (o de manera equivalente ys , ys , ∇ y s , ..., ∇k−1 y s donde ∇ . . . es el operador de diferencias por atr´ as, ∇ y s =y s − y s−1 ). Otra selecci´on es los coeficientes ζsik en el polinomio (6.1), que corresponden a los coeficientes . .. (k) de una serie de Taylor, o a las derivadas zs (0) , z s (0) , z s (0) , ..., zs (0) . Cualquiera que sea la relaci´ on entre estos polinomios y la soluci´on, se utilizan para predecir el valor de la soluci´ on en ts+1 = ts + h, por medio de (p)

ys+1 = zs (h) .

(6.3)

En seguida se utiliza la ecuaci´ on diferencial en el tiempo ts+1 para evaluar la  . (p) (p) (p) y derivada s+1 = f ts+1 , ts+1 y se determina otro polinomio zs+1 (t − ts+1 ) . (p)

con la informaci´ on anterior y el valor de ys+1 , que se utiliza para corregir el valor predicho, obteniendo valores corregidos (c)

(p)

ys+1 = zs+1 (0) = ys+1 = zs+1 (0) .



(c) y s+1 =z s+1 (0) = f ts+1 , ys+1 .

(6.4)



(6.5)

Finalmente, se completa el paso desechando la informaci´ on guardada que no se necesita para definir zs+1 (τ ), y definiendo esto en t´erminos de lo dem´as de la informaci´ on guardada y de la nueva informaci´ on de (6.4),(6.5). Ejemplo 6.1.1 El m´etodo Adams del segundo orden utiliza la f´ ormula corrector .  1 . (c) ys+1 = ys + h y s+1 + y s 2 que . es .equivalente a la regla del trapecio. Inicialmente se tienen valores de ys , y s , y s−1 que definen un polinomio predictor .

zs (τ ) = ys + τ y s +

 . τ2  . y s − y s−1 . 2h

De esto se obtiene  . 1  . (p) ys+1 = zs (h) = ys + h 3 y s − y s−1 . 2

46

Hay varias maneras de aplicar la f´ ormula corrector. Es posible simplemente evaluar 

(p) y (p) s+1 = f ts+1 , ys+1 .



y definir ys+1 = .

(c) ys+1



= ys + (c)

1 2h





. y (p) s+1

+

. ys



ys+1 = f ts+1 , ys+1 Esto es una manera expl´ıcita, clasificado como una aplicaci´ on de modo PECE; lo que quiere decir que las operaciones que se llevan a cabo son (p)

Predicci´ on de

ys+1 . y (p) s+1

Evaluaci´ on de (c) Correcci´ on para obtener ys+1  . (c) (c) y s+1 = f ts+1 , ys+1 Evaluaci´ on de Otra manera de considerar la f´ ormula corrector es como un conjunto de ecuaciones impl´ıcitas (y generalmente no-lineales) .  1   (c) (c) ys+1 = ys + h f ts+1 , ys+1 + y s 2 que deben resolverse por iteraci´ on. Se toman los valores expl´ıcitos anteriores (c) (1) on ys+1 , y se generan aproximaciones de ys+1 como una primera aproximaci´ sucesivas .  1   (j) (j−1) ys+1 = ys + h f ts+1 , ys+1 + y s , j = 2, 3, ... 2 continuando hasta satisfacer un criterio de convergencia, como en el proceso de sustituci´ on sucesiva con el corrector de Euler impl´ıcito (3.23) que es el corrector de Adams de primer orden. Esto puede clasificarse como una aplicaci´ on de modo P (EC)m , donde se toma m tan grande como sea necesario para convergencia, si la iteraci´ on es convergente.

En la pr´actica, cuando se utiliza un tal m´etodo se restringe h para conseguir una convergencia r´ apida, por ejemplo en no m´ as de tres iteraciones. El m´etodo general de Adams es similar. La f´ormula predictor es

47

(p)

ys+1 = ys + h

k 

.

βk∗ y s+1−k

(6.6)

k+1

ormula de que puede derivarse, para obtener los constantes βk∗ , como una f´  ts +h . . . . y (t) dt en t´erminos de y s , y s−1 , ..., y s+1−k . La f´ ormula integraci´ on para ts corrector es (c)

ys+1 = ys + h

k 

.

βk ys+1−k

(6.7)

k=0

que puede escribirse .

(c) (p) ys+1 = ys+1 + hβ0 y s+1 +h

K 

.

(βk − βk∗ ) y s+1−k

(6.8)

k=1

y puede aplicarse en modo expl´ıcito PECE (o un modo P (EC)M E con M fijo > 1), o en modo impl´ıcito P (EC)m con iteraci´on hasta convergencia. Puede ser interesante de escribir en forma expl´ıcita el algor´ıtmo para el modo PECE:

ys+1 = ys + h

K 

 .

βk ys+1−k +hβ0 f ts+1 , ys + h

k=1

K 

 . βk∗ y s+1−k

(6.9)

k=1 .

y s+1 = f (ts+1 , ys+1 )

(6.10)

Se ha supuesto desde un principio que la f´ ormula predictor es de orden K, o sea que ser´ıa exacta si los componentes del vector soluci´on fueran polinomios en t de grado no mayor que K. Es posible que la f´ ormula corrector sea de orden K o de orden K + 1. Una manera sencilla de ver esto es de considerar el corrector de Adams (6.7) como una f´ ormula de integraci´ on . . . ts +h y. y y y (t) dt en t´erminos de s+1 , s , ..., s+1−K . Si se utiliza toda para dar ts esta informaci´ on, pueden escogerse los K + 1 coeficientes βk , k = 0, 1, ..., K . y para dar resultados exactos si (t) es un polinomio de grado K, o sea si la soluci´ on y(t) es un polinomio de .grado K + 1. Por otra parte, si se toma βK = 0, lo que permite desechar y s+1−K antes de aplicar el corrector, y no despu´es, y as´ı reduce la memoria necesaria, se reduce el orden de la f´ormula corector hasta K. Otros m´etodos, que no son del tipo Adams, utilizan un corrector de orden K con βK
= 0, utilizando el grado de libertad para mejorar la regi´ on de estabilidad absoluta. 48

Por otra parte, si el corrector es del tipo Adams y de orden K, el predictor no puede ser de orden mayor que K sin aumentar la memoria utilizada. Ciertos expertos consideran que hay una ventaja utilizar un predictor y un corrector del mismo orden, como esto permite utilizar la llamada ”t´ecnica de Milne” para estimar el error de truncamiento local. Otros prefieren utilizar esta estimaci´on para mejorar la soluci´ on por extrapolaci´ on, lo que es equivalente al utilizar un corrector de orden K + 1, y tener una estimaci´on del u ´ltimo t´ermino de la serie de Taylor incluido. La t´ecnica de Milne para estimar el error es precisamente un proceso de on exacta de la ecuaci´ on extrapolaci´ on. Si Yi+1 es el valor en ts+1 de la soluci´ diferencial que pasa por el punto (ts , ys ), se tiene (p)

ys+1 = Ys+1 + PK+1 hK+1 y (K+1) (ts + θ ∗ h) ,

(6.11)

(c)

ys+1 = Ys+1 + CK+1 hK+1 y (K+1) (ts + θh) ,

(6.12)

ormulas donde PK+1 , CK+1 son constantes conocidos, que dependen de las f´ umeros desconocidos entre 1 − K y 1. particulares utilizadas, y θ ∗ , θ son n´ La diferencia, (c)



(p)



ys+1 − ys+1 = (CK+1 − PK+1 ) hK+1 y (K+1) + 0 hK+2 y (K+2) ,

(6.13)

de donde el error de truncamiento local del corrector de orden K es (c)

ys+1 − Ys+1 =

  CK+1 (c) (p) ys+1 − ys+1 + 0 hK+2 y (K+2) , CK+1 − PK+1

(6.14)

lo que permite obtener una estimaci´ on desde la diferencia entre los valores predichos y corregisos. Es normal ahora, implementar una estrategia de ajuste del tama˜ no del paso para guardar este valor estimado del error local dentro de ciertos l´ımites, de acuerdo con las consideraciones de la secci´on anterior. Si se especifica un vector de tolerancias ε por paso, sus componentes se comparan con los de la estimaci´on del error Ts =

 CK+1 (c) (p) ys+1 − ys+1 . CK+1 − PK+1

(6.15)

Si us = max |Tsi /εi | ≤ 1 se acepta el paso, si no se rechaza. Se calcula un nuevo tama˜ no de paso igual a ξ multiplicado por el valor anterior, donde etodos con pasos m´ ultiples, esto requiere la creaci´on de ξ K+1 < u−1 s . Para m´ 49

un pasado artificial para los pasos siguientes, que no es trivial para ciertas representaciones del polinomio z (τ ). Por esto se evita cambiar el tama˜ no de paso con demasiada frecuencia, y ciertos programas hacen un cambio solamente cuando us > 1 o es posible tomar un ξ = 2, aunque estas condiciones se consideran demasiado fuertes ahora. Ejemplo 6.1.2 Considerese el predictor de Euler, de primer orden .

(p)

ys+1 = ys + h y s y el corrector de Euler impl´ıcito iterado hasta convergencia 

(c)

(c)



ys+1 = ys + hf ts+1 , ys+1 . Se tiene ..



...

..



...

(p) ys+1 = Ys − 12 h2 y s +0 h3 y (c) ys+1 = Ys + 12 h2 y s +0 h3 y

de donde ..



...

(c) (p) ys+1 − ys+1 = h2 y s +0 h3 y

y la estimaci´ on del error de truncamiento local es 1  (c) (p) ys+1 − ys+1 . 2

Ts =

6.1.2

Estabilidad de los m´ etodos de pasos m´ ultiples.

Para investigar la estabilidad de los m´etodos de pasos m´ ultiples, se considera la ecuaci´ on .

y = λy + g (t)

(6.16)

que puede representar un componente de un sistema, con λ igual a un valor propio de la matriz Jacobiana A(t). Se podr´ıa considerar la ecuaci´ on .

y = λy

(6.17)

como antes, y la generalidad de (6.16) no introduce nada nuevo, como se ve en adelante. Sup´ ongase que la soluci´ on correcta de (6.16) es 50

y = S(t)

(6.18)

as´ı que .

g(t) =S (t) − λS (t)

(6.19)

y que la soluci´ on calculada es ys = S(ts ) + εs

(6.20)

y s = λys + g (ts ) . = λS (ts ) + λεs + S (ts ) − λS (ts ) . = S (ts ) + λεs

(6.21)

as´ı que .

y el error global εs satisface la ecuaci´on (6.17). Para el m´etodo impl´ıcito de tipo Adams de segundo orden (la regla del trapecio del ejemplo 6.1.1) iterado hasta convergencia, se tiene 1 (6.22) εs+1 = εs + h [λεs+1 + λεs ] + Ts+1 2 un si S(t) es un polinomio donde Ts+1 es el error de truncamiento local. A´ on (6.22) tiene soluciones cuadr´ atico y Ts+1 = 0, la ecuaci´ ε0 ξ s εs =   ξ = 1 + 12 hλ / 1 − 12 hλ

(6.23)

que pueden ser diferentes de cero. Si se empieza con valores iniciales exactos an un entonces ε0 = 0, pero de todos modos, errores de redondeo introducir´ a error intermedio diferente de cero, y en seguida un Ts+1
= 0 constituir´ otro fuerte de errores, los cuales se propagar´ an de acuerdo con (6.23). Si hλ > 0, hay una funci´ on complementaria creciente eλt , que crece por un hλ no, ξ es una buena aproximaci´ on factor e en cada paso. Para hλ peque˜ hλ a e , pero cuando hλ tiende a 2, ξ tiende a infinito. Se ve que, con hλ positivo (o con parte real positiva, si se trata de un sistema) el m´etodo es menos estable que el problema, pero controlando el tama˜ no del paso, es posible evitar que ξ tome valores grandes. Si la parte real de hλ es negativa, entonces |ξ| < 1 y se tiene estabilidad absoluta para toda la mitad izquierda del plano complejo, una propiedad que se llama A-estabilidad. Sin embargo, cuando hλ → −∞, ξ → −1, y hay una oscilaci´on finita para problemas muy 51

”stiff”. (Para conseguir convergencia, como en (6.22), para |hλ| grande, es esencial utilizar el proceso iterativo de Newton-Raphson, como en el caso del m´etodo impl´ıcito de Euler). Considerese ahora un m´etodo impl´ıcito de Adams de orden mayor, o sea el corrector (6.7). En lugar de las ecuaciones (6.23), se tiene ξ K − ξ K−1 − hλ

K 

βk ξ K−k = 0

(6.24)

k=0

una ecuaci´ on de grado K en ξ. Una de sus ra´ıces (la ra´ız principal) es una no. Las otras (ra´ıces ficticias) son buena aproximaci´ on a ehλ para |hλ| peque˜ todas cero para hλ = 0, y entonces para hλ peque˜ no y diferente de cero satisfacen 

ξ K−1 + hλβK = 0 (hλξ) + 0 ξ K



(6.25)

on de estabilidad de donde son del orden |hλβK |1/K−1 en magnitud. La regi´ absoluta es la regi´on dentro de la cual todas las ra´ıces de (6.24) tienen valor absoluto menor que uno. T´ıpicamente, para orden mayor que dos, esta regi´ on es de la forma siguiente

Es necesario que h sea tal que hλ caiga dentro de esta regi´ on, y para problemas ”stiff” esto quiere decir valores muy chiquitos de h a´ un cuando la soluci´ on var´ıa lentamente. 52

Para un m´etodo Adams en modo PECE (6.9) se tiene

ξ

K

−ξ

K−1

− hλ

K 



βk ξ

K−k

− hλβ0 ξ

K−1

+ hλ

k+1

K 



βk∗ ξ K−k

=0

(6.26)

k=1

una ecuaci´ on con propiedades similares a los de (6.24). Sin embargo, tienen regiones de estabilidad absoluta un poco m´as grande para valores mayores de K, y por esto varios usuarios prefieren el modo PECE de los m´etodos Adams, y no la iteraci´ on hasta convergencia, cuando no se trata de problemas ”stiff”. Para un m´etodo de pasos m´ ultiples m´ as general que los de Adams, con corrector de la forma

(c)

α0 ys+1 = −

K 



(c)



αk ys+1−k + hβ0 f ts+1 , ys+1 + h

k=1

K 

.

βk ys+1−k

(6.27)

k=1

la ecuaci´ on correspondiente a (6.24) es K  k=0

αk ξ K−k − hλ

K 

βk ξ K−k = 0

(6.28)

k=0

o sea α (ξ) − hλβ (ξ) = 0

(6.29)

definiendo los polinomios α (ξ) , β (ξ) con coeficientes αk , βk , k = 0, 1, ..., K. Si el m´etodo es de orden K, se sabe que la ra´ız principal de (6.29) satisface 

ξ = ehλ + 0 (hλ)K+1



(6.30)

de donde con hλ = 0, α (1) = 0. Resulta que puede escogerse el polinomio β (ξ) con ra´ıces en donde sea, y la condici´on (6.30) determina el polinomio α (ξ) de manera u ´nica. Los m´etodos de Adams resultan si se toma β (ξ) de tal manera que las ra´ıces ficticias de (6.29) tienden a cero cuando hλ → 0, o sea si las ra´ıces ficticias de α (ξ) son cero. Para problemas ”stiff” se quiere mantener estabilidad cuando hλ tiende a −∞ lo que resulta posible para orden K hasta 11.

53

Chapter 7

FORMULACION DE NORDSIECK Y METODOS PARA ECUACIONES ”STIFF” 7.1 7.1.1

Una t´ ecnica para cambiar el paso, y las f´ ormulas de Gear La formulaci´ on de Nordsieck.

Para simplificar la notaci´ on en esta secci´on, se considera una sola ecuaci´ on diferencial .

y = f (t, y)

(7.1)

on guardada en el y se reserva la notaci´on ys para el vector de informaci´ m´etodo predictor-corrector. Para el m´etodo de Adams de orden K, los componentes de ys ser´ıan         

ys0 ys1 − − − ysK





    =   

       

ys . h ys − − − . h y s+1−K

54

        

(7.2)

donde se han incluido multiplicadores h, para las derivadas, para simplificar f´ ormulas tales como (6.7). Si en realidad se tiene un sistema de ecuaciones, es posible sustituir un arreglo de dos dimensiones para este vector, con una columna para cada componente de la soluci´ on, pero la notaci´ on resulta muy compleja Dado ys , es f´acil construir un vector predicho, multiplicando por una matriz apropiada, por ejemplo       (p)  ys+1 =      



(p) ys+1 . (p) h y s+1 .

h ys − − − . h y s+2−K



            =           

1 0 0 0 − − − 0

β1∗ γ1 1 0 − − − 0

β2∗ γ2 0 1 − − − 0

−−− −−− −−− −−− −−− −−− −−− −−−

∗ βK−1 ∗ γK−1 0 0 − − − 1

∗ βK ∗ γK 0 0 − − − 0

        ys = Cys.      

(7.3) Aqu´ı, los γk son coeficientes nuevos, que todav´ıa pueden escogerse para . (p) consistencia, y ys+1 claramente no tiene el mismo significado como antes. . (p) En efecto, y s+1 debe ser la derivada de zs (τ ) en τ = h. Si se tiene otro (p)

vector de informaci´ on guardada, la interpretaci´ on de los vectores ys , ys+1 es diferente, y la matriz C tiene que ser diferente. (p) (0) on iterativa Se tomar´ a ys+1 como ys+1 , el valor inicial para la aplicaci´ del corrector. Una aplicaci´ on del corrector (6.7), o de su forma m´ as general (6.27), puede escribirse Fm . (m+1)

h y s+1 (m+1) ys+1

. (m)



(m)

= h y s+1 −hf ts+1 , ys+1



. (m)

(7.4)

= h y s+1 −Fm (m) = ys+1 − Fm β0 /α0 . (m+1)

donde las primeras dos ecuaciones corresponden a la evaluaci´ on, ys+1 es lo . (m) . (p) (o sea y s+1 para la primera aplicaci´ on con m = 0). que antes se llamaba y s+1 . A cualquier cambio a h ys+1 corresponde el mismo cambio, multiplicado por β0 /α0 para seguir satisfaciendo (6.27). El resultado de aplicar este corrector una o varias veces es (c)

(p)

ys+1 = ys+1 − F q 55

(7.5)

donde, en este caso, q es el vector       q=     

β0 /α0 1 0 − − − 0

           

(7.6)

(donde α0 = 1 para un m´etodo Adams), y F es un escalar que depende del modo en lo cual se aplica el corrector. Finalmente, tomando (c)

ys+1 = ys+1

(7.7)

se ha avanzado un paso. Por la manera de incluir la evaluaci´ on en la correcci´on siguiente, el modo de aplicaci´on es necesariamente P (EC)m , con m fijo ≥ 1 o con iteraci´on hasta convergencia. Consid´erese ahora un cambio de la informaci´ on guardada, sin cambiar de manera esencial el m´etodo. Si el nuevo vector de informaci´ on guardada ´ nica el es zs , cualquier de los dos vectores ys , zs especifica de manera u mismo polinomio de grado K, y debe haber una matriz P no-singular que les relaciona, o sea = P zs ys ys+1 = P zs+1

(7.8)

En t´erminos de los nuevos vectores, (7.3) se escribe (p)

zs+1 = P −1 CP zs

(7.9)

y (7.5) se convierte en (c)

(p)

zs+1 = zs+1 − P −1 qF con el mismo valor de F . En el caso particular

56

(7.10)

      zs =      

ys . h ys .. 1 2 y s 2! h − − − 1 K (K) K! h ys

           

(7.11)

el polinomio predictor es zs (t − ts ) =

K 



zsj

j=0

t − ts h

j

(7.12)

y se tienen sus derivadas escaladas K  1 k (k) h zs (τ ) = k! j=k



j k



 j−k

zsj

τ h

.

(7.13)

Con τ = h se obtiene la matriz del predictor        −1 P CP = A =       

1 0 0 0 − − − 0

1 1 0 0 − − − 0

1 2 1 0 − − − 0

1 3 3 1 − − − 0

−−− −−− −−− −−− −−− −−− −−− −−−

1 K K(K−1) 2 K(K−1)(K−2) 2·3

− − − 1

             

(7.14)

donde 

Akj =

j k



= Ak,j−1 + Ak−1,j−1 ,

j ≥ k.

(7.15)

Escribiendo l = P −1 q se tiene la formulaci´on de Nordsieck (p)

zs+1 = Azs (p) zs+1 = zs+1 − lF en donde se normaliza l por la condici´ on l1 = 1. 57

(7.16)

Con esta notaci´ on, no hay problema ninguno en cambiar el tama˜ no del paso h cuando es necesario, como esto implica solamente un nuevo escalamiento de los componenetes de zs por potencias de h. En la representaci´ on (7.3) del m´etodo de Adams, no se definieron los on de C, pero es f´ acil de encontrar sus valores, elementos γk del segundo rengl´ considerando la consistencia del primer corrector, (7.4) con m = 0:  (1) ys+1

=

(p) ys+1



+ β0 hf

(p) ts+1, ys+1



−

k 

 .

γk h y s+1−k

(7.17)

k=1

de donde, comparando con (6.8) γk = (βk∗ − βk ) /β0 .

(7.18)

Ejemplo 7.1.1 Para el m´etodo Adams se segundo orden β2∗ = − 12 β1∗ = 32 , 1 β0 = β1 = 2 , β2 = 0 α0 = 1, α1 = −1, α2 = 0 1 2,

de donde γ1 = 2, γ2 = −1 

1 32  C= 0 2 0 1  

1 2



− 12  −1  0  

q= 1  0 Para polinomios de segundo grado 









ys ys 1 0 0 . .       ys=  h y s  =  0 1 0   h y s  = P zs .. . 1 2 ys 0 1 −2 h y s−1 2h de donde

58



1  P −1 CP =  0 0

0 1 1 2



0 1  0  0 − 12 0





3 2



− 12 1 0 0   2 −1   0 1 0  1 0 0 1 −2



1 1 1   = 0 1 2 =A 0 0 1 como debe de ser, y 

1  l = P −1 q =  0 0

0 1

1 2









1 1 0  2   2  0  1  =  1 . 1 − 12 0 2

Este vector caracteriza el m´etodo impl´ıcito de Adams de segundo. orden en la formulaci´ on de Nordsieck. Cuando se hace una correcci´ on a h y s+1 (o sea una evaluaci´ on) se hace la mitad de la misma correcci´ on a ys+1 , y tambi´en .. 1 2 y on polinomial de segundo a 2 h s+1 . Esto quiere decir que la aproximaci´ ultiplo de grado en ts+1 se corrige por un m´ 1 t − ts+1 1 + + 2 h 2



t − ts+1 h

2

1 = 2



t − ts h

2

y es obvio que esto es consistente con el proceso de guardar los valores de . ys , h y s sin cambio para el paso siguiente. Adem´as, se ve que F es la suma de varias correcciones a la derivada, una para cada iteraci´ on del corrector, y de (7.16), para un m´etodo dado (o sea un vector l dado, como A es fijo), F es proporcional a la diferencia entre los valores corregidos y predichos y entonces proporcional al error de truncamiento local.

7.1.2

Estabilidad y m´ etodos para ecuaciones ”stiff”.

Es claro que para ecuaciones ”stiff” se necesita un m´etodo impl´ıcito, como m´etodos expl´ıcitos tienen β0 = 0 y al menos una de las ra´ıces de (6.28) crece sin l´ımite cuando hλ tiende a −∞. Se necesita un m´etodo impl´ıcito para lo cual todas las ra´ıces de α (ξ) − hλβ (ξ) quedan con valor absoluto menor que uno para hλ real negativo, cualquiera que sea su magnitud, lo que implica que las ra´ıces de β (ξ) no deben ser mayores que uno. La regi´ on

59

de estabilidad absoluta de tal m´etodo incluir´ıa todo el eje real negativo, y ser´ıa de la forma general del exterior de una curva as´ı

.

Ser´ıa ideal que el m´etodo des A-estable como los m´etodos de Euler por atr´ as, y del trapecio: esto implica que la curva se encontrar´ıa completamente del lado derecho del eje imaginario, pero no es posible para K > 2. La regi´ on de inestabilidad del lado izquierdo del eje imaginario corresponde a funciones complementarias oscilantes con amortiguamiento d´ebil, para las cuales la restricci´on sobre el paso h es normalmente m´as fuerte debido a precisi´on, que la debido a la estabilidad. Gear propuso los m´etodos de orden K con β (ξ) = ξ K , o sea con todas sus ra´ıces iguales a cero, y demostr´o que ten´ıan regiones de estabilidad absoluta de la forma descrita para K ≤ 6. Otros investigadores han encontrado m´etodos con las mismas propiedades hasta K = 11. Para todos estos m´etodos, queda un problema ya mencionado para los m´etodos de Euler por atr´ as, y del trapecio. Consid´erese la iteraci´on del corrector (7.4) para la ecuaci´on lineal 60

.

y = λy + g (t) .

(7.19)

Se tiene . (m+1)

Fm+1 = h y s+1 . (m)



(m+1)

−hf ts+1 , ys+1 

(m+1)

= h y s+1 −Fm − h λys+1 

(m)



 



+ g (ts+1 )



(m)



= hf ts+1 , ys+1 − h λ ys+1 − Fm β0 /α0 + g (ts+1 ) = hλFm β0 /α0

(7.20)

y si |hλ| < |α0 /β0 | las iteraciones del corrector divergen. Es necesario reconocer que el objeto de las iteraciones es de encontrar una F en (7.16) tal que zs+1 satisface la ecuaci´on diferencial, lo que quiere decir que zs+1,1 − hf (ts+1 , zs+1,0 ) = 0

(7.21)

y de resolver esta ecuaci´on por el m´etodo iterativo de Newton-Raphson. Se tiene (m+1)



(m+1)

zs+1,1 − hf ts+1 , zs+1,0





(m)

(m)

= zs+1,1 − l1 Fm − hf ts+1 , zs+1,0 − l0 Fm 

(m)

(m)

 zs+1,1 − hf ts+1 , zs+1,0 







− l1 I − l0 hA(m) (t) Fm (7.22) 2 , y con A(m) (t) definido como en (3.27). Por despreciando t´erminos en Fm (7.21), esta expresi´ on debe de ser cero aproximadamente, dando 



. (m)



(m)

l1 I − l0 hA(m) (t) Fm = h y s+1 −hf ts+1 , ys+1



(7.23)

un sistema de n ecuaciones lineales para determinar Fn (en el caso de un sistema de n ecuaciones diferenciales), en lugar de la definici´ on en (7.4). Si el sistema es lineal, y se tiene la matriz Jacobina exacta, una iteraci´on con (7.23) da la soluci´ on de (7.21). En general, no se tiene la matriz Jacobiana exacta, como no se quiere repetir el c´alculo e inversi´ on (factorizaci´ on) de   (m) (t) en cada paso (a´ un menos en cada iteraci´ on la matriz l1 I − l0 hA para un sistema no-lineal), si no es relativamente sencillo. Por lo tanto, es necesario normalmente, iterar el corrector (7.23) m´as de una vez. Gear sugeri´ o la siguiente estrategia: 61

(a) Si hay convergencia del corrector en 3 iteraciones o menos,   se sigue (m) (t) . utlitizando la misma factorizaci´ on de la matriz l1 I − l0 hA (b) Si en 3 iteraciones no hay convergencia, y si no se hab´ıa calculado la matriz inmediatamente antes, se calcula y se factoriza nuevamente, y se repiten las iteraciones. (c) Si no hay convergencia con una nueva matriz, se reduce h en un factor 16, se calcula la nueva matriz, y se repiten las iteraciones. Con un poco m´ as memoria y programaci´ on, ser´ıa posible mejorar la  (m) estrategia, guardando A (t) para facilitar el calculo de l1 I − l0 hA(m) (t) cuando h o K (que determina l0 , l1 ) cambian. En caso (b) un nuevo c´ alculo (m) de A (t) ser´ıa necesario de todos modos. Ejemplo 7.1.2 Para la f´ ormula de diferenciaci´ on por atr´ as (de Gear) de segundo orden, se pone β0 = 1 ,

β1 = 0 ,

β2 = 0

y resulta 3 1 , α1 = −2 , α2 = 2 2 para que el m´etodo sea exacto para cualquiera soluci´ on polinomial de segundo grado. Entonces la informaci´ on guardada es α0 =











ys ys 1 0 0 . .      1 0   h y s  = P zs ys =  h y s  =  0 .. 1 2 y 1 −1 1 ys−1 s 2h y los predictores son .

(p) ys+1 = 2h y s +ys−1 . (p) . h y s+1 = −2ys + 3h y s +2ys−1

de donde la matriz 



0 2 1   C =  −2 3 2  . 1 0 0 Se verifica que P −1 CP = A, y como

62









2 β0 /α0    3  q= 1 = 1  0 0

se tiene 









2 2 1 0 0   3   3  −1 l = P q =  0 1 0  1  =  1 . 1 0 −1 1 1 3

ultiplo de Cualquiera correcci´ on a zs+1 es un m´ 2 t − ts+1 1 + + 3 h 3



t − ts+1 h

2

=

1 3



t − ts h





t − ts−1 . h .

Entonces deja sin cambio los valores de ys , ys−1 , pero en este caso no h y s que ya no se utiliza en el paso siguiente. Al comparar los m´etodos de Gear y de Adams, se encuentra que los de Gear necesitan mas tiempo de c´alculo para problemas sin ”stiffness”, y mucho menos tiempo para problemas ”stiff”. Algunos investigadores consideran que las ventajas de los m´etodos de Adams en los casos sin ”stiffness” no son suficientemente grandes para justificar la utilizaci´ on de los dos, y sugieren los m´etodos de Gear para fines generales. Sin emabargo, los m´etodos de Adams son todav´ıa m´as conocidos y generalmente disponibles y otros trabajadores consideren que es necesario tener rutinas para los dos m´etodos, dejando al usuario decidir si su problema es ”stiff”. Puede concluirse que en caso de duda es mejor utilizar los m´etodos de Gear, pero ciertas implementaciones de los m´etodos de Adams salen con un diagn´ ostico si consideran que el sistema de ecuaciones es ”stiff”.

63

Chapter 8

METODOS AUTOMATICOS. 8.1 8.1.1

Dise˜ no de rutinas con un m´ınimo de par´ ametros de control. Necesidades para un algoritmo autom´ atico.

Ya se ha mencionado varias veces la idea de m´etodos autom´aticos de integraci´ on, en particular la selecci´on autom´ atica del paso h en base en un error de truncamiento local estimado. Se estudiar´ an ahora con m´ as detalle las necesidades para un algoritmo autom´ atico. Deben incluirse los siguientes puntos: (A) un m´etodo para avanzar la soluci´ on de un tiempo dado t por un paso hasta tiempo t + h, y estimar el error de truncamiento local correspondiente; (B) una t´ecnica para cambiar el tama˜ no del paso h (y eventualmente el orden K), y para actualizar la informaci´ on guardada que depende de h (o de K); (C) una estrategia para decidir el valor de h (y eventualmente de K) que debe utilizarse, como resultado de la utilizaci´ on de (A) para un paso de prueba, y un examen del error estimado; (D) una manera de entregar al usuario los valores y(Tr ), r = 0, 1, 2, ..., de la soluci´on en tiempos Tr , especificados de antemano, y frecuentemente con intervalos iguales. Se ha considerado (A) en las secciones anteriores. Hay que acordarse que normalmente existe una constante K (el orden del m´etodo) tal que el

64

error de truncamiento local, para un paso de tama˜ no h, es 

dr = CK+1 hK+1 y (K+1) (tr ) + O hK+2



(8.1)

ormula puede donde CK+1 es un constante caracter´ıstico del m´etodo. Tal f´ tener validez, por cierto, s´ olo si las ecuaciones diferenciales son tales que y(t) tiene una derivada continua de orden K + 1, o sea que f (t, y) tiene derivadas continuas de orden K con respecto a sus argumentos. A´ un si esto es cierto para la ecuaci´on diferencial exacta que quiere resolverse, debe cuidarse que cualquiera aproximaci´ on a las funciones involucradas es suficientemente suave para no violar esta condici´on, si se quiere que un m´etodo autom´ atico funcione sin problemas. Sup´ ongase por ejemplo que se quiere resolver la ecuaci´on sencilla .

y +αy = g(t)

(8.2)

donde g(t) es una funci´ on f´ısica compleja para la cual se tienen aproximaciones diferentes en los rangos 0 ≤ t ≤ 1, 1 ≤ t ≤ 4, 4 ≤ t. A´ un si se cuida que los errores de las aproximaciones sean peque˜ nas, y que sean continuas entre s´ı en los puntos t = 1, t = 4, de todos modos no tendr´ an una suavidad perfecta en estos puntos; una derivada, en muchos casos la primera, tendr´ a una discontinuidad. Es claro que la derivada de la soluci´ on de orden uno mayor ser´a tambi´en discontinua, y se puede esperar que un m´etodo autom´atico de soluci´ on de alto orden no podr´ıa atravesar esas discontinuidades, o al menos generar´ıa errores inesperadamente grandes, o tomar´ıa pasos muy peque˜ nos as´ı desgastando el tiempo de c´alculo. El punto (B) es trivial para un m´etodo de un solo paso. Para un m´etodo de pasos m´ ultiples, puede ser relativamente sencillo para la formulaci´on de Nordsieck, o bastante dif´ıcil para otras formulaciones. Por esta raz´ on, muchos programas autom´aticos en el pasado han utilizado pasos h relacionados por potencias de 2, o sea de la forma 2p h0 para un entero p (positivo o negativo) y un paso de prueba h0 , posiblemente dado por el usuario, o calculado en base a un criterio del programa. El mejor valor de h en cualquier punto de la integraci´ on puede ser de la forma 2p+θ h0 para 0 ≤ θ < 1, y si se hace la suposici´ on razonable que θ puede tomar cualquier valor en este intervalo con la misma probabilidad se concluye que tal m´etodo hace aproximadamente 1.4 veces m´as trabajo que lo necesario. M´etodos m´as modernos no restringen los cambios de h a factores de dos, y reducen el trabajo a aproximadamente 1.1 veces m´as que lo ideal, o sea 30% menos que antes.

65

El punto (C) es lo m´ as importante para un m´etodo autom´atico, y se considera con m´ as detalle despu´es. Antes se va a ver el punto (D) que puede parecer trivial a primera vista. La soluci´ on obvia es de probar, antes de hacer un paso de la integraci´ on, si se va a rebasar el punto siguiente Tr en lo cual se necesita la soluci´on, y en el caso que s´ı, reducir h hasta Tr − t. Esto funciona bien si el espaciamiento entre los puntos Tr es grande comparado con el paso h´ optimo, en cual caso no cuesta mucho c´alculo suplementario, y si no hay complicaciones al cambiar h por un factor arbitrario. Muchas veces esto no es el caso. Para m´etodos de pasos m´ ultiples, a´ un con la formulaci´ on de Nordsieck, hay problemas si se cambia el paso h con frecuencia, o por factores grandes. Adem´ as muchas veces el usuario quiere valores de la soluci´on en nos que los pasos h que puntos Tr , con intervalos iguales pero m´as peque˜ pueden utilizarse por razones de precisi´ on y estabilidad. Esto puede ser el caso para problemas ”stiff”, si el usuario quiere resultados con intervalos t´ıpicos de las funciones complimentarias r´ apidas, a´ un cuando la soluci´ on particular se ha suavizada. Consid´erese los distintos tipos de m´etodos por separado. Los m´etodos de un solo paso no tienen, desde luego, problemas para cambiar el valor del paso h que utilizan. Para problemas ”stiff” tienen que tomar pasos chiquitos (al menos los m´etodos que se han visto en este curso), y es probable que necesitar´ıan muchos pasos entre valores consecutivos de Tr . Los m´etodos predictor-corrector pueden utilizar la formulaci´ on de Nordsieck, que facilita el cambio del paso, pero estas f´ormulas tienen un m´etodo muy econ´omico de interpolar dentro de un paso, como en cada paso se tiene un polinomio interpolante. Pueden, entonces, evitar el problema tomando un paso normal que rebasa uno (o varios) puntos Tr , y despues interpolar la soluci´ on en estos puntos. Por ejemplo, con la formulaci´ on de Nordsieck, puede utilizarse el polinomio predictor y(t) =

  K  1 k (k) t − ts k h ys k! h

(8.3)

k=0

para interpolaci´ on en el rango ts < t < ts + h, o igualmente utilizar de manera atrasada el predictor siguiente   K  1 k (k) t − ts+1 k h ys+1 y(t) = k=0

k!

h

(8.4)

despu´es de completar el paso. Esto tiene ventajas bien claras sobre (8.3), como da una soluci´on continua; (8.4) evaluado en t = ts da ys , mientras 66

(p)

(8.3) evaluado en t = ts+1 da ys+1
= ys+1 . Hay otra posibilidad, que es de interpolar todos los valores guardados, tanto en ts como en ts+1 ; por un polinomio de grado m´ as elevado, que da una soluci´ on interpolada m´ as suave. Para los m´etodos de un solo paso, no parece que hay una manera de utilizar los valores intermedios dentro del paso, que generan los m´etodos de Runge-Kutta y de extrapolaci´ on, como esos valores no tienen el mismo orden de precisi´on como los valores del principio y al final del paso. En una rutina autom´ atica publicada, utilizando el m´etodo de extrapolaci´ on de Bulirsch y Stoer, cada vez que se ha pasado por un punto Tr , se guardan los valores actuales de la soluci´ on, y se ejecuta un paso m´ as corto por atr´ as, hasta Tr . Con este m´etodo, el paso adicional puede ser bastante costoso. Consid´erese ahora la estrategia para fijar h (y eventualmente K, si se deja la libertad de variarlo). Es importante realizar que casi todos los m´etodos autom´ aticos utilizan estrategias relativamente sencillas, para evitar c´alculos adicionales costosos. Escogen h (dado el orden K) lo m´as grande posible para guardar, en cada paso, una norma del error de truncamiento local estimado dentro de ciertos l´ımites. Tomando en cuenta las suposiciones (a) y (b) de la secci´ on 5.1.5, una norma t´ıpica es    di   Y 

d = max1≤i