Metodos Num´ericos en Ecuaciones Diferenciales Ordinarias

Tema 4

M´etodos Num´ericos en Ecuaciones Diferenciales Ordinarias

4.1 Introducci´on

Estudiaremos en este Tema algunos m´etodos num´ericos para resolver problemas de valor inicial en ecuaciones diferenciales ordinarias y en sistemas de e.d.o.

4.2 M´etodo de Euler

El M´etodo de Euler o de las Tangentes constituye el primer y m´as sencillo ejemplo de m´etodo num´erico para la resoluci´on de un problema de valor inicial:

y0 = f(x,y), y(x0) = y0

donde suponemos adem´as que se verifican las hip´otesis del Teorema de Picard1, y en consecuencia existe soluci´on u´nica para el problema. Interpretando la e.d.o. y0 = f(x,y) como un campo de direcciones en el plano x−y y la condici´on inicial y(x0) = y0 como un punto (x0,y0) de dicho plano, podemos aproximar la funci´on soluci´on y(x) por medio de la recta tangente a la misma que pasa por ese punto: y(x)∼ = y0 +f(x0,y0)(x−x0) 1Consideraremos en general que la funci´on f(x,y) es diferenciable en un entorno del punto (x0,y0). Si bien es cierto que se trata de una condici´on m´as restrictiva de lo estrictamente necesario, en la pr´actica trabajaremos siempre con funciones de ese tipo.

59

60 TEMA 4

donde se ha utilizado que la pendiente de dicha tangente es: m = y0(x0) y, en conse- cuencia: m = f(x0,y0). Calculamos as´ı de manera aproximada el valor de la soluci´on y en el punto de abscisa x1 como: y(x1)∼ = y1 = y0 +f(x0,y0)(x1 −x0) y con este punto (aproximado) ya calculado, podemos repetir el m´etodo para obtener otro punto aproximado (x2,y2) de la forma: y(x2)∼ = y2 = y1 +f(x1,y1)(x2 −x1) y as´ı sucesivamente. Es habitual en este m´etodo tomar abscisas equiespaciadas, es decir, calcular la soluci´on aproximada en puntos de la forma: xn = xn−1 + h = x0 + nh, siendo h el paso del m´etodo. De esta forma se obtienen las f´ormulas que nos determinan la soluci´on aproximada en la forma: xn = xn−1 +h; yn = yn−1 +f(xn−1,yn−1)h Desde el punto de vista geom´etrico, tenemos en definitiva que el M´etodo de Euler aproxima a la funci´on soluci´on por medio de una l´ınea poligonal, la aproximaci´on ser´a tanto peor cuanto mayor sea en nu´mero de pasos, es decir, cuanto m´as “lejos” nos encontremos del punto inicial (x0,y0). Por otro lado, el error ser´a evidentemente tanto mayor cuanto m´as grande sea el “paso” del m´etodo, h. Ejemplo: Resolveremos el problema de valor inicial ( y0 = x√y y(1) = 4 por el m´etodo de Euler con h = 0.1 para los puntos x = 1.1,1.2,1.3,1.4 y 1.5. En este problema tenemos h = 0.1, (x0,y0) = (1,4) y la funci´on f(x,y) es: f(x,y) = x√y. Por tanto: yn = yn−1 +xn−1√yn−1 h Dado que el problema se puede resolver tambi´en de forma exacta, presentamos en la tabla y gr´afica siguientes los resultados:

i xi yi Sol. Exacta 0 1 4 4 1 1.1 4.2 4.21276 2 1.2 4.42543 4.45210 3 1.3 4.67787 4.71976 4 1.4 4.95904 5.01760 5 1.5 5.27081 5.34766 1.1 1.2 1.3 1.4 1.5 4.2 4.4 4.6 4.8 5 5.2

M´ETODOS NUM´ERICOS EN EDO 61

4.3 M´etodos de Taylor

El M´etodo de Euler que acabamos de describir no es m´as que un caso particular de los m´etodos de Taylor, que consisten de manera general en aproximar la soluci´on por su polinomio de Taylor de un orden determinado. Partiendo por tanto del P.V.I.:

y0 = f(x,y) y(x0) = y0 ) talquepresentasoluci´onu´nicay(x)enunentornode x0 (soluci´onquesuponemosadem´as derivable n veces en dicho entorno), aproximaremos dicha funci´on por su polinomio de Taylor de orden N:

y(x)'y(x0)+y0(x0)(x−x0)+

1 2

y00(x0)(x−x0)2+

1 3!

y000(x0)(x−x0)3+...+y(N) N!

(x−x0)N y el error de aproximaci´on viene determinado por el resto de orden N + 1, de manera que el error es proporcional a (x−x0)N+1. Si fijamos una sucesi´on de puntos equiespaciados : x0,x1,x2,..., con xn+1 = xn +h, y denominamos (de manera similar a lo hecho en la secci´on anterior) y0,y1,... a los valores paroximados correspondientes de y(x), tendremos que:

yn+1 = yn +y0(xn)h+

1 2

y00(xn)h2 +...+

1 N!

y(N)(xn)hN

con un error en cada paso (error local2) proporcional a hN+1. Parapoderaplicarelm´etodonecesitamosconocerlasderivadasdelasoluci´on(recorde- mos que desconocida), pero teniendo en cuenta la propia ecuaci´on diferencial:

y0(xn) = f(xn,yn)

mientras que y00(xn) puede ser calculada derivando: y00 = d dx f(x,y) = ∂f ∂x + ∂f ∂y dy dx =

∂f ∂x

+

∂f ∂y

f(x,y)

y as´ı sucesivamente. Es evidente que el m´etodo de Taylor de orden uno no es m´as que el M´etodo de Euler antes descrito. Analicemos con un ejemplo el M´etodo de Taylor de orden dos: Ejemplo: Apliquemos el m´etodo de Taylor de orden dos a la ecuaci´on y0 = cos(xy), con la condici´on inicial: y(0) = 1. La expresi´on a considerar ser´a: y(xn+1)'yn+1 = yn +y0(xn)h+ 1 2 y00(xn)h2 2En general se llama “error local” al cometido en cada uno de los pasos de un m´etodo num´erico. El error acumulado tras n pasos recibe el nombre de “error global” y si el local es proporcional a hp, entonces el global ser´a proporcional a hp−1.

62 TEMA 4

tendremos entonces:

y0(xn) = f(x,yn) = cos(xnyn)

y00 = dcos(xy) dx

=−sen(xy)y−sen(xy)y0 =−sen(xy)(y +cos(xy)) Y as´ı los primeros pasos de la resoluci´on ser´an (tomaremos h = 0.5):

x1 = x0 +h = 0.5 y1 = y0 +hcos(x0y0)+

h2 2

(−sen(x0y0)(y0 +cos(x0y0))) =

= 1+0.5cos0+

0.52 2

(−sen0(1+cos(0))) = 1.5

x2 = x1 +h = 1 y2 = y1 +hcos(x1y1)+

h2 2

(−sen(x1y1)(y1 +cos(x1y1))) = = 1.5+0.5cos(0.51.5)+ 0.52 2 (−sen(0.51.5)(1.5+cos(0.51.5))) = 1.67569 y as´ı sucesivamente. El error local en esta aproximaci´on ser´a proporcional a h3 y por tanto el error global lo ser´a a h2.

4.4 M´etodos de Runge-Kutta

La idea general de los M´etodos de Runge-Kutta es sustituir el Problema de Valor Inicial:

y0 = f(x,y) y(x0) = y0 )

por la ecuaci´on integral equivalente: Zy y0 dy =Zx x0 f(x,y(x))dx⇒y = y0 +Zx x0

f(x,y(x))dx

para proceder a aproximar esta u´ltima integral mediante un m´etodo num´erico adecuado (recordemos que y(x) es desconocida). Si nuevamente planteamos el problema “paso a paso” tendremos: yn+1 = yn +Zxn+1 xn f(x,y(x))dx

4.4.1 M´etodo de Runge-Kutta de segundo orden

La primera opci´on que podemos aplicar es integrar mediante el m´etodo de los trapecios, es decir tomando: Zxn+1 xn f(x,y(x))dx' 1 2 h (f(xn,yn)+f(xn+1,yn+1))

M´ETODOS NUM´ERICOS EN EDO 63

Al ser desconocida yn+1 en la expresi´on anterior, lo aproximaremos por ¯ yn+1, donde

¯ yn+1 es la estimaci´on de yn+1 que resultar´ıa aplicando el m´etodo de Euler. Tendremos as´ı: Zxn+1

...