• Author / Uploaded
• xabihdez

Citation preview

Algoritmos básicos para la resolución de ecuaciones diferenciales ordinarias A la hora de establecer algoritmos para la resolución de ecuaciones diferenciales ordinarias es preciso tener en cuenta la casuística de cada problema. Así se establece la división entre los problemas con condiciones iniciales y los problemas con condiciones de contorno. En los primeros los valores de la función solución se van construyendo sucesivamente a partir de valores previamente calculados, comenzando en las condiciones iniciales. Los problemas con condiciones de contorno pueden resolverse, en su mayoría, seleccionando la solución correcta de entre todas las obtenidas con el apoyo de los algoritmos de condiciones iniciales. a) Problemas de condiciones iniciales. En lo que sigue se describen los fundamentos de los algoritmos más sencillos para resolver una ecuación diferencial ordinaria de primer orden sujeta a cierta condición inicial. La generalización a sistemas de ecuaciones de primer orden y/o ecuaciones diferenciales de orden superior es inmediata. Consideremos una ecuación diferencial del tipo:

y'(t) =

dy(t) = f[t, y(t)] dt

ta≤t≤tb

y(ta)=ya

Conocido el valor de la función en un punto yi≡y(ti), podemos calcular el valor yi+1 a partir del desarrollo en serie de Taylor alrededor de yi. Llamando h≡ti+1-ti (en lo que sigue asumiremos que h es constante) se tiene:

h2 hn y(t ) = y(t ) + hy'(t ) + y''(t ) + L + y(n) (t ) + O(hn+1 ) i i i i i+1 2 n! La aproximación más simple consiste en cortar el desarrollo en primer orden:

y(t ) = y(t ) + hy'(t ) + O(h2 ) i i i+1 de este modo se obtienen los valores de y(t) por propagación a partir de y(ta). Este es el Método de Euler que, utilizando la notación más habitual, puede reescribirse como:

wo ≡ ya

wi + 1 = wi + hf(ti, wi )

t −t h≡ b a N

i = 0,L, N - 1

El método de Euler requiere muchos valores {es decir muchas evaluaciones de f[t,y(t)]} de h pequeños para obtener una precisión aceptable. Desgraciadamente el error aumenta linealmente con N y, por tanto, salvo para tanteos iniciales el método de Euler es de poca utilidad. La solución más obvia consiste en retener términos de orden superior en el desarrollo de Taylor mencionado.

wo ≡ ya

wi + 1 = wi + hT (n) (ti, wi )

i = 0,L, N - 1

donde:

h h2 hn - 1 (n - 1) T (n) (t , w ) ≡ f(t , w ) + f'(t , w ) + f''(t , w ) + L + f (t , w ) + O(hn ) i i i i 2 i i i i i i 2 n! Sin embargo, el algoritmo directo requiere evaluaciones de la función f[t,y(t)] y de sus derivadas. Los algoritmos conocidos como Runge-Kutta simplifican en gran medida el cálculo de T(n)(ti,wi), sin perder precisión, encontrando la respuesta a la siguiente pregunta: ¿Es posible aproximar el valor de T(n)(ti,wi) a partir del cálculo de la función derivada, f[tRK,yRK(tRK)], en varios puntos dentro del intervalo [ti,ti+1]?. Para explicar la idea consideremos T(2). Por definición: ⎫ h ⎧ ∂f(t, y) ∂f(t, y) T (2) (t, y) ≡ f(t, y) + ⎪⎨ + f(t, y)⎪⎬ 2 ⎪⎩ ∂t ∂y ⎪⎭

Por otra parte es fácil ver que dadas las constantes desconocidas a1, α1 y β1 se tiene: ⎡

a1f(t + α1, y + β1 ) ≈ a1 ⎢f(t, y) + ⎢⎣

Y, por comparación, si a1=1, α1=

∂f(t, y) ∂f(t, y) ⎤ β⎥ α + 1 1 ⎥⎦ ∂t ∂y

h h y β1= f(t, y): 2 2

h h T (2) (t, y) ≈ f[t + , y + f(t, y)] 2 2 Como consecuencia se obtiene el algoritmo de Runge-Kutta de orden 2 ó del punto medio:

wo ≡ ya

h h wi + 1 = wi + hf[ti + , wi + f(ti, wi )] 2 2

i = 0,L, N - 1

¿Es válida esta técnica para órdenes superiores?. La respuesta es (afortunadamente) sí. Por ejemplo si pasamos a cuarto orden [T(4)(x,y)] puede repetirse un proceso similar (aunque bastante más largo) al efectuado con T(2)(x,y) y así obtener el algoritmo de Runge-Kutta de orden 4 que es frecuentemente utilizado y que, de manera compacta se muestra a continuación:

wo ≡ ya

k1 ≡ hf(ti, wi )

k h k2 ≡ hf(ti + , wi + 1 ) 2 2 k h k3 ≡ hf(ti + , wi + 2 ) 2 2 k4 ≡ hf(ti + 1, wi + k3 )

1 wi + 1 = wi + (k1 + 2k2 + 2k3 + k4 ) 6

i = 0,L, N - 1

El método de Runge-Kutta mejora si se utilizan algoritmos más refinados como el control de paso (h) adaptativo (que pueden estudiarse en algunos de los libros recomendados en la bibliografía). Aunque no es el método definitivo puede usarse siempre que no se requiera una gran eficiencia computacional, algún método más sofisticado falle o no se conozca nada mejor. Los algoritmos de resolución de una ecuación diferencial ordinaria de primer orden con condiciones iniciales, pueden generalizarse directamente a un sistema de ecuaciones:

du (t) 1

= f1[t, u1 (t), u2(t),Lum (t)]

2

= f2[t, u1 (t), u2(t),Lum (t)]

dt du (t) dt

M

du (t) m

= fm[t, u1 (t), u2(t),Lum (t)] dt u1 (ta ) = u1a , u2 (ta ) = u2a ,L, um (ta ) = uma ; ta < t < tb resolviendo secuencialmente cada una de ellas. En particular debe tenerse en cuenta que las ecuaciones diferenciales de orden superior: y(m) (t) ≡

dm y(t) dt

m

= f[t, y(t), y'(t),L y(m - 1) (t)]

y(ta ) = ya, y'(ta ) = y'a ,L, y(m - 1) (ta ) = ya(m - 1) ; ta < t < tb

pueden ser convertidas en sistemas de ecuaciones de primer orden como sigue:

du (t)

dy(t) = u2(t) dt dt du (t) dy'(t) 2 = = u3(t) dt dt 1

=

M

du

m -1

(t)

dt du (t) m

dt

=

=

dy

dy

(m - 2)

(t)

dt

(m - 1)

dt

(t)

= um (t)

= ym (t) = f[t, y(t), y'(t),L y(m - 1) (t)] = f[t, u1 (t), u2(t),Lum - 1 (t)]

u1 (ta ) = y(ta ) = ya, u2(ta ) = y'(ta ) = y'a ,L, um (ta ) = y(m - 1) (ta ) = ya(m - 1) ; ta < t < tb Siendo:

u1 (t) ≡ y(t); u2 (t) ≡ y'(t); L ;um (t) ≡ y(m - 1) (t) b) Problemas con condiciones de contorno. Para resolver este tipo de problemas pueden emplearse diferentes métodos. El que va a esbozarse aquí es conceptualmente simple y de aplicación general. Es el denominado método del disparo. Considérese una ecuación diferencial elemental como y’’(x)=6x cuya solución debe satisfacer como condiciones y(0)=0 e y(1)=2. Los métodos de condiciones iniciales basados en desarrollos de Taylor, propagan

los valores de la función y(x) a partir de la condición inicial y(0). Es decir 2 nos permiten obtener funciones que satisfacen una de las y''(x)=6x condiciones de contorno aunque 3 y(x)=x +ax y(0)=0 y(1)=2 sin garantía de que la otra sea también satisfecha. En el caso simple que nos ocupa, el método de Runge-Kutta nos proporcionaría una familia de funciones y(x)=x3+ax de las que sólo una (la que corresponde a a=1) satisfaría la condición y(1)=2. En la figura se representa esta situación. La solución podríamos obtenerla 0 “disparando” desde y(0) con 0 1 diferente valores de a [es decir y’(0)] hasta “dar en el blanco” alcanzando el valor requerido por la condición de contorno en x=1. Obviamente esto puede hacerse manualmente por tanteo, pero también puede automatizarse. Analicemos paso a paso el algoritmo: 1) Dado que la ecuación es de 2º orden, obtengamos el sistema asociado de ecuaciones de primer orden. Definiendo u1≡y(x) y u2≡y’(x), las ecuaciones buscadas son: u’1(x)=u2(x) u’2(x)=6x con u1(0)=y(0)=0 y u2(0)=y’(0)=ap. 2) Resolvamos el sistema para un cierto de valor de a (ap) utilizando un algoritmo de condiciones iniciales (por simplicidad de notación el método de Euler, aunque sea poco recomendable): u1(xi+1;ap)=u1(xi;ap)+hu2(xi;ap) u2(xi+1;ap)=u2(xi;ap)+6hx

h≡xi+1-xi; i=0,…,N-1

3) Si u1(xN=1;ap)=y(1) hemos resuelto el problema. Normalmente este no será el caso y deberemos cambiar el valor del parámetro a. Este

cambio puede hacerse muy fácilmente si se tiene en cuenta que el valor correcto del parámetro de disparo es una raíz de la ecuación u1(xN=1;a)-y(1)=0. Así podemos obtener un nuevo valor ap+1 utilizando el método de Newton-Raphson:

(

) ( )

u x = 1; ap − y( 1 ) ap+1 = ap − 1 N du1 a da p du

( )

El problema puede ser calcular 1 ap , pero puede obtenerse con da relativa facilidad resolviendo un sistema de ecuaciones diferenciales ordinarias adicionales. 4) Con el nuevo valor se repiten los pasos 1) a 3) hasta que alcanza la convergencia deseada. Por supuesto el método funcionará con una buena primera aproximación para a1. El algoritmo general se esquematiza en el siguiente diagrama.