1

Deducción de las Ecuaciones del Método de Runge-Kutta El problema consiste en encontrar una solución numérica a la ecuación diferencial dy ordinaria de primer orden: = f ( x, y ) sujeta a la condición inicial y ( x0 ) = y 0 . El dx objetivo consiste en encontrar aproximaciones satisfactorias para los valores de la solución y ( x ) en un conjunto especificado de valores de x denotados como x1 , x 2 , x3 , x 4 " . Los valores exactos los denotaremos como y ( x1 ), y ( x 2 ), y ( x3 )" y sus valores aproximados los denotaremos como y1 , y 2 , y3 , y 4 " Empezamos aproximando el valor de y en el punto x1 = x0 + ∆x , esto es, y ( x1 ) = y 0 + ∆y . La manera mas simple de hacer esto consiste en aproximar ∆y por la estimación acostumbrada para el verdadero incremento, o sea ∆y ≅ dy = y ′( x0 )∆x . La ecuación diferencial misma nos da el valor de la derivada en el punto ( x0 , y 0 ) , esto es, y ′( x0 ) = f ( x0 , y 0 ) ; y entonces ∆y ≅ f ( x0 , y 0 )∆x y por lo tanto

y ( x1 ) = y 0 + ∆y ≅ y 0 + f ( x0 , y 0 )∆x = y1 Una vez que y1 ha sido obtenida como la aproximación a y ( x1 ) , el mismo procedimiento puede repetirse en ( x1 , y1 ) para obtener y ( x 2 ) = y ( x1 ) + ∆y ≅ y1 + f (x1 , y1 )∆x = y 2 ; y así sucesivamente hasta donde se necesite. Este método se conoce con el nombre de “Método de Euler” Ahora bien, habiendo obtenido y1 como una primera aproximación a y ( x1 ) por el método de Euler, podemos usar ahora la ecuación diferencial para calcular y ′ en el nuevo punto P1 : ( x1 , y1 ) y usar entonces el promedio de las derivadas en los puntos P0 : ( x0 , y 0 ) y P1 : ( x1 , y1 ) para obtener una mejor aproximación de ∆y , y por tanto de y ( x1 ) antes de calcular la siguiente aproximación y ( x 2 ) . Este método nos da el 1 valor ∆y ≅ [ y ′( x0 ) + y ′( x1 )]∆x , y la estimación mejorada del siguiente punto es 2

1 [ f (x0 , y0 ) + f (x1 , y1 )]∆x = ( y1 )2 . Este proceso se conoce 2 con el nombre de “Método Modificado de Euler”

entonces: y 0 + ∆y ≅ y 0 +

Otra posibilidad adicional, después de haber obtenido y1 como una primera aproximación a y ( x1 ) por el método de Euler, consiste en reaproximar ∆y y y ( x1 ) usando la derivada en el punto medio de P0 : ( x0 , y 0 ) y P1 : ( x1 , y1 ) en lugar de usar el promedio de las derivadas, esto es en el punto:

2

 x + x1 y 0 + y1  M : 0 ,  . Esto nos da la estimación mejorada 2   2 1 ∆x  x + x1 y 0 + y1    , , y 0 + f ( x0 , y 0 )∆x ∆x , y esto nos da una ∆y ≅ f  0 ∆x = f  x0 + 2  2 2    2 1 ∆x   , y 0 + f ( x0 , y 0 )∆x ∆x . tercera aproximación a y ( x1 ) como: ( y1 )3 = y 0 + f  x0 + 2 2   Este proceso se conoce con el nombre de “Método de Runge” El método de Runge-Kutta es básicamente una generalización de esos tres procedimientos simples en el que en cada paso se calculan tres o mas estimaciones de ∆y . El valor de ∆y que se usa entonces para calcular el siguiente valor de y es una combinación lineal de esas estimaciones en la cual las constantes de combinación se escogen para hacer el error tan pequeño como sea posible. En el método de Runge-Kutta de tercer orden se toman los siguientes tres estimados de ∆y : (∆y )1 = f (x0 , y 0 )∆x que es la estimación del método de Euler; (∆y )2 = f [x0 + p∆x, y 0 + p(∆y )1 ]∆x ; 0 < p < 1 ; que es parecido al estimado del 1 método de Runge excepto que en lugar de evaluar en el punto medio (con p = ) la 2 derivada se calcula en un punto P : [x0 + p∆x, y 0 + p(∆y )1 ] , que todavía no se ha determinado; y (∆y )3 = f [x0 + q∆x, y 0 + r (∆y )2 + (q − r )(∆y )1 ]∆x ; 0 < q , r < 1 en donde q y r deben de calcularse. Finalmente el valor de ∆y que se usa finalmente para calcular y1 se toma como: (∆y )4 = a(∆y )1 + b(∆y )2 + c(∆y )3 en donde a,b,c son parámetros que igual que los parámetros p,q,r deben escogerse para dar la mas alta precisión al estimar ∆y . Los detalles de este cálculo pueden consultarse en el libro “Advanced Engineering Mathematics” de C.R.Wylie; Third Edition, McGraw-Hill, 1966. Dicho procedimiento nos lleva a un sistema de cuatro ecuaciones con seis 1 1 1 ; p 2b + q 2 c = ; prc = incógnitas: a + b + c = 1 ; pb + qc = de donde 2 3 6 pueden despejarse cuatro de los parámetros en función de los otros dos, obteniendo: 6 pq − 3( p + q ) + 2 2 − 3q 2 − 3p q(q − p ) ; b= ; c= a= ; r= 6 pq 6 p( p − q ) 6q(q − p ) p(2 − 3 p ) Puesto que p y q son arbitrarias, tenemos entonces una familia de fórmulas biparámetricas que pueden usarse para resolver la ecuación diferencial de primer 4 orden con una precisión del orden de (∆x ) .

3

Dos casos especiales del “Método de Runge-Kutta de Tercer Orden” vale la pena anotar. Para listarlos usaremos la notación convencional: ∆x = h ; (∆y )1 = k1 ; (∆y )2 = k 2 ; (∆y )3 = k 3

1 3 ; b=0 ; c= ; 4 4 1 ∆y ≅ (∆y )4 = (k1 + 3k 3 ) 4 k1 = f ( x0 , y 0 )h

a=

CASO 1:

p=

1 2 ; q=r= 3 3

1 1   k 2 = f  x0 + h, y 0 + k1 h 3 3   2 2   k 3 = f  x0 + h, y 0 + k 2 h 3 3   1 3 2 ; b=c= ; p=q=r= 4 8 3 1 ∆y ≅ (∆y )4 = (2k1 + 3k 2 + 3k 3 ) 8 k1 = f ( x0 , y 0 )h

a=

CASO 2:

2 2   k 2 = f  x0 + h, y 0 + k1 h 3 3   2 2   k 3 = f  x0 + h, y 0 + k 2 h 3 3   La discusión anterior puede extenderse sin dificultad, (excepto los detalles), para llegar a procedimientos de solución en los cuales el error es del orden de h 5 = (∆x )5 En particular los dos conjuntos de fórmulas siguientes son bastante útiles:

1 (k1 + 2k 2 + 2k 3 + k 4 ) 6 k1 = f ( x0 , y 0 )h

∆y ≅ (∆y )5 =

CASO 3:

1  1  k 2 = f  x0 + h, y 0 + k1 h 2  2  1 1   k 3 = f  x0 + h, y 0 + k 2 h 2 2   k 4 = f ( x0 + h, y 0 + k 3 )h

4

1 (k1 + 3k 2 + 3k 3 + k 4 ) 8 k1 = f ( x0 , y 0 )h ∆y ≅ (∆y )5 =

CASO 4:

1  1  k 2 = f  x0 + h, y 0 + k1 h 3  3  2 1   k 3 = f  x0 + h, y 0 − k1 + k 2 h 3 3   k 4 = f ( x0 + h, y 0 + k1 − k 2 + k 3 )h

El proceso de solución basado en el CASO 3 se conoce usualmente como “El Método de Runge-Kutta”

5

Ejemplo de solución de una Ecuación Diferencial de Primer Orden por el Método de Runge-Kutta: y'=2xy n 0 1 2 3 4 5 6 7 8 9 10

Xn 1.00 1.10 1.20 1.30 1.40 1.50 1.60 1.70 1.80 1.90 2.00

Yn 1.0000 1.2337 1.5527 1.9937 2.6116 3.4902 4.7586 6.6188 9.3923 13.5969 20.0813

k1 0.20000 0.27141 0.37265 0.51836 0.73126 1.04706 1.52274 2.25040 3.38121 5.16682 8.03251

k2 k3 0.23100 0.23426 0.31496 0.31997 0.43475 0.44252 0.60827 0.62041 0.86341 0.88257 1.24426 1.27483 1.82157 1.87088 2.71041 2.79091 4.10066 4.23375 6.31032 6.53331 9.87998 10.25872

k4 kprom 0.27154 0.23367 0.37287 0.31902 0.51876 0.44099 0.73195 0.61795 1.04826 0.87858 1.52481 1.26834 2.25401 1.86028 3.38751 2.77342 5.17788 4.20465 8.05208 6.48436 12.74279 ########

6

k1 = hf ( x n , y n )

Método de Runge – Kutta ejemplo resuelto con detalles y ′ = 2 xy ;

0

f ( x, y ) = 2 xy ;

x n +1 = x n + h

y (1) = 1 ; h = 0.10

y n +1 = y n + k prom

x1 = 1.0000

y1 = 1.0000

y1 = 1.0000 f = 2.0000 k1 = 0.2000

x 2 = 1.1000 2

x 2 = 1.1000

y 2 = 1.2337

y 2 = 1.2337 f = 2.714083 k1 = 0.27141

x3 = 1.2000 3

x3 = 1.2000

y 3 = 1.5527

4

x 4 = 1.3000

y 4 = 1.9937

k prom =

y 3 = 1.5527 f = 3.72648 k1 = 0.37265

1 (k1 + 2k 2 + 2k 3 + k 4 ) 6

k2

k3

h = 1.0500 2 k y1 + 1 = 1.1000 2 f = 2.3100 k 2 = 0.2310

h = 1.0500 2 k y1 + 2 = 1.1155 2 f = 2.34255 k 3 = 0.234255

h = 1.1500 2 k y 2 + 1 = 1.3694 2 f = 3.14962 k 2 = 0.314962

h = 1.1500 2 k y 2 + 2 = 1.3911 2 f = 3.199716 k 3 = 0.319971

h = 1.2500 2 k y 3 + 1 = 1.7390 2 f = 4.34756 k 2 = 0.434756

h = 1.2500 2 k y 3 + 2 = 1.7701 2 f = 4.42519 k 3 = 0.44252

k1

x1 = 1.0000 1

k  h  k 2 = hf  x n + , y n + 1  2 2 

x1 +

x2 +

x3 +

x1 +

x2 +

x3 +

k  h  k 3 = hf  x n + , y n + 2  2 2  k 4 = hf ( x n + h, y n + k 3 )

k4

k prom

x1 + h = 1.1000 y1 + k 3 = 1.2342 f = 2.715361 k 4 = 0.271536

k prom = 0.233674

x 2 + h = 1.2000 y 2 + k 3 = 1.5537 f = 3.7288 k 4 = 0.37288

k prom = 0.319024

x3 + h = 1.3000 y 3 + k 3 = 1.99522 f = 5.187572 k 4 = 0.51876

k prom = 0.440994