Metodos num,

Tema 4 Métodos Numéricos en Ecuaciones Diferenciales Ordinarias 4.1 Introducción Estudiaremos en este Tema algunos métodos numéricos para resolver problemas de valor inicial en ecuaciones diferenciales ordinarias y en sistemas de e.d.o. 4.2 Método de Euler El Método de Euler o de las Tangentes constituye el primer y más sencillo ejemplo de método numérico para la resolución de un problema de valor inicial: y ′ = f (x, y) , y(x0 ) = y0 donde suponemos además que se verifican las hipótesis del Teorema de Picard1 , y en consecuencia existe solución única para el problema. Interpretando la e.d.o. y ′ = f (x, y) como un campo de direcciones en el plano x − y y la condición inicial y(x0 ) = y0 como un punto (x0 , y0 ) de dicho plano, podemos aproximar la función solución y(x) por medio de la recta tangente a la misma que pasa por ese punto: y(x) ∼ = y0 + f (x0 , y0 )(x − x0 ) 1 Consideraremos en general que la función f (x, y) es diferenciable en un entorno del punto (x0 , y0 ). Si bien es cierto que se trata de una condición más restrictiva de lo estrictamente necesario, en la práctica trabajaremos siempre con funciones de ese tipo. 59 60 TEMA 4 donde se ha utilizado que la pendiente de dicha tangente es: m = y ′ (x0 ) y, en consecuencia: m = f (x0 , y0 ). Calculamos ası́ de manera aproximada el valor de la solución 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étodo 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étodo tomar abscisas equiespaciadas, es decir, calcular la solución aproximada en puntos de la forma: xn = xn−1 + h = x0 + nh, siendo h el paso del método. De esta forma se obtienen las fórmulas que nos determinan la solución aproximada en la forma: xn = xn−1 + h; yn = yn−1 + f (xn−1 , yn−1 ) h Desde el punto de vista geométrico, tenemos en definitiva que el Método de Euler aproxima a la función solución por medio de una lı́nea poligonal, la aproximación será tanto peor cuanto mayor sea en número de pasos, es decir, cuanto más “lejos” nos encontremos del punto inicial (x0 , y0 ). Por otro lado, el error será evidentemente tanto mayor cuanto más grande sea el “paso” del método, h. Ejemplo: Resolveremos el problema de valor inicial ( √ y′ = x y y(1) = 4 por el método 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ón 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én de forma exacta, presentamos en la tabla y gráfica siguientes los resultados: i 0 1 2 3 4 5 xi 1 1.1 1.2 1.3 1.4 1.5 yi 4 4.2 4.42543 4.67787 4.95904 5.27081 Sol. Exacta 4 4.21276 4.45210 4.71976 5.01760 5.34766 5.2 5 4.8 4.6 4.4 4.2 1.1 1.2 1.3 1.4 1.5 61 MÉTODOS NUMÉRICOS EN EDO 4.3 Métodos de Taylor El Método de Euler que acabamos de describir no es más que un caso particular de los métodos de Taylor, que consisten de manera general en aproximar la solución por su polinomio de Taylor de un orden determinado. Partiendo por tanto del P.V.I.: ) y′ = f (x, y) y(x0 ) = y0 tal que presenta solución única y(x) en un entorno de x0 (solución que suponemos además derivable n veces en dicho entorno), aproximaremos dicha función por su polinomio de Taylor de orden N : 1 1 y (N ) y(x) ≃ y(x0 )+y ′ (x0 )(x−x0 )+ y ′′ (x0 )(x−x0 )2 + y ′′′ (x0 )(x−x0 )3 +. . .+ (x−x0 )N 2 3! N! y el error de aproximación 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ón de puntos equiespaciados : x0 , x1 , x2 , . . ., con xn+1 = xn + h, y denominamos (de manera similar a lo hecho en la sección anterior) y0 , y1 , . . . a los valores paroximados correspondientes de y(x), tendremos que: 1 (N ) 1 y (xn )hN yn+1 = yn + y ′ (xn )h + y ′′ (xn )h2 + . . . + 2 N! con un error en cada paso (error local2 ) proporcional a hN +1 . Para poder aplicar el método necesitamos conocer las derivadas de la solución (recordemos que desconocida), pero teniendo en cuenta la propia ecuación diferencial: y ′ (xn ) = f (xn , yn ) mientras que y ′′ (xn ) puede ser calculada derivando: y ′′ = d ∂f ∂f dy ∂f ∂f f (x, y) = + = + f (x, y) dx ∂x ∂y dx ∂x ∂y y ası́ sucesivamente. Es evidente que el método de Taylor de orden uno no es más que el Método de Euler antes descrito. Analicemos con un ejemplo el Método de Taylor de orden dos: Ejemplo: Apliquemos el método de Taylor de orden dos a la ecuación y ′ = cos(xy), con la condición inicial: y(0) = 1. La expresión a considerar será: 1 y(xn+1 ) ≃ yn+1 = yn + y ′ (xn )h + y ′′ (xn )h2 2 2 En general se llama “error local” al cometido en cada uno de los pasos de un método numérico. El error acumulado tras n pasos recibe el nombre de “error global” y si el local es proporcional a hp , entonces el global será proporcional a hp−1 . 62 TEMA 4 tendremos entonces: y ′ (xn ) = f (x, yn ) = cos(xn yn ) y ′′ = d cos(xy) = − sen(xy) y − sen(xy) y ′ = − sen(xy)(y + cos(xy)) dx Y ası́ los primeros pasos de la resolución serán (tomaremos h = 0.5): x1 = x0 + h = 0.5 y1 = y0 + h cos(x0 y0 ) + = 1 + 0.5 cos 0 + h2 (− sen(x0 y0 )(y0 + cos(x0 y0 ))) = 2 0.52 (− sen 0(1 + cos(0))) = 1.5 2 x2 = x1 + h = 1 y2 = y1 + h cos(x1 y1 ) + h2 (− sen(x1 y1 )(y1 + cos(x1 y1 ))) = 2 0.52 (− sen(0.5 1.5)(1.5 + cos(0.5 1.5))) = 1.67569 = 1.5 + 0.5 cos(0.5 1.5) + 2 y ası́ sucesivamente. El error local en esta aproximación será proporcional a h3 y por tanto el error global lo será a h2 . 4.4 Métodos de Runge-Kutta La idea general de los Métodos de Runge-Kutta es sustituir el Problema de Valor Inicial: ) y′ = f (x, y) y(x0 ) = y0 por la ecuación integral equivalente: Z Z x Z y f (x, y(x)) dx ⇒ y = y0 + dy = y0 x f (x, y(x)) dx x0 x0 para proceder a aproximar esta última integral mediante un método numérico adecuado (recordemos que y(x) es desconocida). Si nuevamente planteamos el problema “paso a paso” tendremos: Z xn+1 f (x, y(x)) dx yn+1 = yn + xn 4.4.1 Método de Runge-Kutta de segundo orden La primera opción que podemos aplicar es integrar mediante el método de los trapecios, es decir tomando: Z xn+1 1 f (x, y(x)) dx ≃ h (f (xn , yn ) + f (xn+1 , yn+1 )) 2 xn 63 MÉTODOS NUMÉRICOS EN EDO Al ser desconocida yn+1 en la expresión anterior, lo aproximaremos por ȳn+1 , donde ȳn+1 es la estimación de yn+1 que resultarı́a aplicando el método de Euler. Tendremos ası́: Z xn+1 1 f (x, y(x)) dx ≃ h (f (xn , yn ) + f (xn+1 , ȳn+1 )) 2 xn con ȳn+1 = yn + h f (xn , yn ) y llegaremos a la expresión del método: yn+1 = yn + h (f (xn , yn ) + f (xn+1 , ȳn+1 )) 2 Lo normal es presentar el método con las expresiones siguientes: k1 = h f (xn , yn ) ; k2 = h f (xn+1 , yn + k1 ) 1 (k1 + k2 ) 2 Comparando este método con el método de Taylor de segundo orden, es posible demostrar que el error local es también proporcional a h3 y, por tanto, el global lo es a h2 . yn+1 = yn + 4.4.2 Método de Runge-Kutta de tercer orden Se trata de la misma idea pero integrando por el Método de Simpson, entonces: Z xn+1 xn f (x, y(x)) dx ≃ h 2 3 ³ f (xn , yn ) + 4f (xn+ 1 , ȳn+ 1 ) + f (xn+1 , yn+1 ) 2 2 ´ donde ȳn+1 e ȳn+ 1 son estimaciones, puesto que yn+1 e yn+ 1 no son conocidos. 2 2 La estimación de ȳn+ 1 se hace por el método de Euler: 2 ȳn+ 1 = yn + 2 h f (xn , yn ) 2 mientras que la estimación de ȳn+1 se pueden considerar varias opciones, por ejemplo: ȳn+1 = yn + h f (xn , yn ) es decir el Método de Euler de nuevo, o por ejemplo: ȳn+1 = yn + h f (xn+ 1 , ȳn+ 1 ) 2 2 que consiste en variar el Método de Euler tomando la pendiente de la recta tangente en el punto medio en vez de la tangente en el punto propiamente dicha. Finalmente, lo más 64 TEMA 4 usual es tomar una combinación de las dos opciones3 : ³ ´ ȳn+1 = yn + h 2 f (xn+ 1 , ȳn+ 1 ) − f (xn , yn ) 2 2 Podemos entonces resumir el Método de Runge-Kutta de tercer orden en la forma: k1 = h f (xn , yn ) h 1 k2 = h f (xn + , yn + k1 ) 2 2 k3 = h f (xn + h, yn − k1 + 2k2 ) 1 yn+1 = yn + (k1 + 4k2 + k3 ) 6 h4 Finalmente, añadir que el error local en el Método de tercer orden es proporcional a y en consecuencia el global lo es a h3 . 4.4.3 Método de Runge-Kutta de cuarto orden Los Métodos de Runge-Kutta de cuarto orden se deducen de una manera similar a la expuesta en la sección anterior para el caso de tercer orden. Ahora se introduce un nuevo paso intermedio en la evaluación de la derivada. Una vez más se presentan varias opciones en la evaluación y es posible ajustar de tal manera que se garantice el error local de manera proporcional a h5 (es decir garantizando exactitud en el cuarto orden en el polinomio de Taylor), lo cual lleva a un error global proporcional a h4 . El Método de cuarto orden más habital es el determinado por las fórmulas siguientes k1 = h f (xn , yn ) h k1 k2 = h f (xn + , yn + ) 2 2 k2 h k3 = h f (xn + , yn + ) 2 2 k4 = h f (xn + h, yn + k3 ) 1 yn+1 = yn + (k1 + 2k2 + 2k3 + k4 ) 6 que al igual que el método de tercer orden está basado en el método de interación de Simpson. Los errores local y global son en este caso proporcionales a h5 y h4 respectivamente. 3 El razonamiento que lleva a tomar esta combinación concreta, que no reproduciremos aquı́, se basa en tomar una combinación genérica: ³ ´ ȳn+1 = yn + h θ f (xn , yn ) + (1 − θ)f (xn+ 1 , ȳn+ 1 ) 2 2 y optimizar el valor de θ imponiendo exactitud en la fórmula de Taylor al tercer orden, ello lleva a que necesariamente θ = −1. 65 MÉTODOS NUMÉRICOS EN EDO Ejemplo: Resolver por un método de Runge-Kutta de cuarto orden el problema de valor inicial: y ′ = x2 − 3y ; y(0) = 1 en el intervalo 0 ≤ x ≤ 0.4, con h = 0.1. Tenemos que x0 = 0, y0 = 1, y f (x, y) = x2 − 3y. Para x1 = 0.1 l ordenada correspondiente será: 0.1 y1 = 1 + (k1 + 2k2 + 2k3 + k4 ) 6 con k1 = f (x0 , y0 ) = 02 − 3 1 = −3 h 0.1 2 0.1 h , y0 + k1 ) = (0 + ) − 3(1 + (−3)) = −2.5475 2 2 2 2 k2 = f (x0 + k3 = f (x0 + h h 0.1 2 0.1 , y0 + k2 ) = (0 + ) − 3(1 + (−2.5475)) = −2.61538 2 2 2 2 k4 = f (x0 + h, y0 + h k3 ) = (0 + 0.1)2 − 3(1 + 0.1(−2.61538)) = −2.20539 y ası́: y1 = 0.741148 De manera análoga se determinan los puntos: (x2 , y2 ) = (0.2, 0.551151) La solución exacta es: ; (x3 , y3 ) = (0.3, 0.413894) ; 25 −3x 1 y= e + 27 3 (x4 , y4 ) = (0.4, 0.317435) µ ¶ 2 2 2 x − x+ 3 9 de manera que: y(0.1) = 0.741127 ; 4.5 y(0.2) = 0.551121 ; y(0.3) = 0.413860 ; y(0.4) = 0.317402 Estabilidad Para finalizar esta sección, analizaremos el comportamiento de Método de Euler en algunas situaciones, que motivarán la necesidad de introducir métodos más precisos. Además del error propio del método (local en h2 y global en h), lo cual no es en principio excesivamente grave si es posible reducir h a voluntad (aunque ello conlleva errores considerables de redondeo y una sobrecarga en el tiempo de computación), existe un problema serio asociado al método y es su inestabilidad. Veamos el siguiente ejemplo: Consideremos la ecuación de Malthus con constante negativa: y ′ = −α y ; y(0) = y0 66 TEMA 4 con α e y0 ambas positivas. La solución exacta es sencilla: y(t) = y0 e−αt El Método de Euler para este caso conduce a la expresión: yn+1 = yn − α yn h = (1 − αh) yn Es fácil comprobar que si αh < 1, entonces la solución numérica es decreciente y positiva, como ocurre en la solución real exacta, pero si αh > 1, entonces el signo de la solución se va alternando y, más aún, si αh > 2, entonces la magnitud de la solución se va incrementando en cada paso (o sea resulta creciente en valor absoluto), y la solución oscila. Este es un comportamiento tı́pico de lo que se conoce en Matemáticas como inestabilidad del método. 4.5.1 Método de Euler Modificado Veremos a continuación una variante del Método de euler, llamado habitualmente Método de Euler Modificado. Se trata de un método más preciso que el de Euler y además más estable. La idea fundamental del método modificado es usar el método de los trapecios para integrar la ecuación y ′ = f (x, y). Tomaremos ası́, en el intervalo [x0 , x − 1]: Z y1 Z x1 (x1 − x0 ) ′ y = f (x, y(x)) ⇔ dy = f (x, y(x))dx ≃ (f (x0 , y0 ) + f (x1 , y1 )) 2 y0 x0 Repitiendo el razonamiento, tendremos yn+1 = yn + h (f (xn+1 , yn+1 ) + f (xn , yn )) 2 Si la función f es lineal en la variable y, entonces es posible despejar fácilmente yn+1 en la ecuación anterior. Si no es lineal, entonces necesitamos un método numérico para calcular la correspondiente yn+1 , tı́picamente se utiliza el método de las sustituciones sucesivas. Veamos un par de ejemplos para aclarar estas ideas: Ejemplo 1: Si consideramos la ecuación lineal: y ′ = ay + cos x entonces: yn+1 = yn + 1+ h (ayn+1 + cos xn+1 + ayn + cos xn ) ⇒ yn+1 = 2 1− ah 2 ah 2 + h 2 1− ah 2 (cos xn+1 + cos xn )) 67 MÉTODOS NUMÉRICOS EN EDO Ejemplo 2: Tomemos ahora el siguiente problema de valor inicial basado en una ecuación no lineal: 3 y ′ = −y 2 + 1 ; y(0) = 10 El método requerirá ahora la resolución de: ´ 3 3 h³ yn+1 = yn + −(yn+1 ) 2 − (yn ) 2 + 2 2 4.6 Método Numéricos para Sistemas Los métodos de Euler y de Runge-Kutta que hemos planteado se aplican de manera sencilla a sistemas de ecuaciones de primer orden, y, en consecuencia, a ecuaciones de orden superior al primero. Comencemos con el Método de Euler aplicado a un sistema de dos ecuaciones de primer orden:  dx = f (t, x, y)   dt   dy = g(t, x, y) dt  x(0) = x0    y(0) = y0 la extensión natural del método ya expuesto nos conduce a las expresiones: tn+1 = tn + h xn+1 = xn + h f (tn , xn , yn ) yn+1 = yn + h g(tn , xn , yn ) Ejemplo: Dado el problema de valor inicial: y ′′ − 0.05 y ′ + 0.15 y = 0, y ′ (0) = 0, y(0) = 1 calcular y(1) e y ′ (1) utilizando el Método de Euler con h = 0.5. Si introducimos una nueva variable dependiente: z = y ′ , y pasamos por tanto de un problema de segundo orden en variables (t, y) a uno con tres variables (t; y, z) tendremos el sistema de ecuaciones de primer orden:   y′ = z   z′ = 0.05z − 0.15y   y(0) = 1    z(0) = 0 Las expresiones del método serán entonces: yn+1 = yn + h f (tn , yn , zn ) = yn + h zn , zn+1 = zn + h g(tn , yn , zn ) = yn + h (0.05zn − 0.15yn ) Y ası́, para t1 = t0 + h = 0 + 0.5 = 0.5: y1 = y0 + h z0 = 1 , z1 = z0 + h (0.05z0 − 0.15y0 ) = −0.075 68 TEMA 4 y para t2 = 1: y2 = y1 + h z1 = 0.96250 z2 = z1 + h (0.05z1 − 0.15y1 ) = −0.15187 , Finalmente, plantearemos en un ejemplo la aplicación de un método de Runge-Kutta de cuarto orden a un sistema de ecuaciones, siendo extensible dicho planteamiento a otros métodos y por supuesto a otro tipo de sistemas. Ejemplo: Planteemos la resolución del sistema: dx dt dy dx = = −4x + 3y + 6 −2x + y + 3 ) con las condiciones iniciales: x(0) = y(0) = 0 por un método de Runge-Kutta de cuarto orden: t0 = 0 x1 = x0 + , x0 = y0 = 0 h h (k1 + 2k2 + 2k3 + k4 ) ; y1 = y0 + (l1 + 2l2 + 2l3 + l4 ) 6 6 k1 = f (t0 , x0 , y0 ) = −4x0 + 3y0 + 6 l1 = g(t0 , x0 , y0 ) = −2x0 + y0 + 3 h h h , x0 + k1 , y0 + l1 ) 2 2 2 h h h l2 = g(t0 + , x0 + k1 , y0 + l1 ) 2 2 2 h h h k3 = f (t0 + , x0 + k2 , y0 + l2 ) 2 2 2 h h h l3 = g(t0 + , x0 + k2 , y0 + l2 ) 2 2 2 k4 = f (t0 + h, x0 + hk3 , y0 + hl3 ) k2 = f (t0 + l4 = g(t0 + h, x0 + hk3 , y0 + hl3 )