¡Descarga Métodos Numéricos para la Resolución de Ecuaciones Diferenciales Ordinarias y más Apuntes en PDF de Métodos Matemáticos solo en Docsity!
Resolución numérica de
Ecuaciones Diferenciales de
Ordinarias (ODEs) de primer
orden: Problema de Valor Inicial
. Contextualización....... . Métodos numéricos de resolu- ción para PVIs........... Errores asociados a los métodos de aproximación de ODEs... . Método de Euler........ Análisis del error en el método de Euler............... Pseudocódigo del Método de Eu- ler................... . Métodos de Heun y del punto medio................ Método de Heun....... Método del punto medio.. . Método de Runge-Kutta IV.
. Contextualización
Numerosos fenómenos, tanto naturales como artificiales, consisten de relaciones entre magnitudes en las que el concepto tasa de cambio juega un papel fundamental. La representación matemática de estas relaciones son, a menudo, ecuaciones en las que dicha tasa de cambio se expresa a través de derivadas. Por tal motivo, este tipo de ecuaciones reciben el nombre de ecuaciones diferenciales o, a veces, ecuaciones de cambio. Así, fenómenos tales como:
I El movimiento de los fluidos. I E flujo de corriente en un circuito. I La disipación del calor en estructuras sólidas. I La propagación y detección de ondas sísmicas. I La dinámica de poblaciones. I ...
se modelan e investigan a través de ecuaciones diferenciales, cuyo estudio (propiedades y soluciones) se hace entonces indispensable.
A modo de ilustración se muestran a continuación dos ecuaciones difer- enciales:
. Ecuación diferencial que modela el flujo de carga, &(C) , en un circuito RLC sometido a una tensión alterna ⇢(C) ,
3 2 &(C)
3C 2
3&(C)
3C
&(C) = ⇢(C) (.)
En (.) la cantidad que cambia, &(C), se denomina variable de- pendiente , mientras que la cantidad con respecto a la cual cambia, C , se denomina variable independiente. Como se puede observar, (.) involucra sólo una variable inde- pendiente. A este tipo de ecuaciones diferenciales se les denomina ecuaciones diferenciales ordinarias.
. Ecuación que describe la propagación de una onda unidimen- sional, 2 2
% 2 D(G, C)
%G 2
% 2 D(G, C)
%C 2
En (.) la variable dependiente es la propiedad de la onda D(G, C), y las variables independientes son, G y C. Como se puede observar, (.) involucra dos variables independi- entes. Las ecuaciones diferenciales que involucran dos o más variables independientes reciben el nombre de ecuaciones diferenciales en derivadas parciales.
primer orden Problema de Valor Inicial
Las ecuaciones diferenciales (ordinarias y en derivadas parciales) se clasifican a su vez por su orden.
Definition .. Orden de una ecuación diferencial El orden de una ecuación diferencial es el orden de la derivada más alta que aparece en la ecuación.
Así, (.) es una ecuación diferencial ordinaria (EDO) de orden , dado que la derivada más alta que aparece es de orden . Sin embargo, la ecuación diferencial que modela el movimiento en caída libre de un cuerpo con resistencia,
3E(C) 3C
E(C) (.)
es una EDO de orden , dado que la derivada más alta involucrada es de orden . De igual forma, (.) es una ecuación en derivadas parciales (EDP) de orden , dado que la derivada parcial más alta que aparece en la ecuación es de orden . Por el contrario, la ecuación diferencial conocida con el nombre de ecuación de advección ,
%D(G, H) %G
%D(G, H)
%H
es una EDP de orden , dado que la derivada parcial más alta involucrada es de orden ^.
Definition .. EDO de orden = En general, la ecuación
[C, H(C), H^0 (C),... , H(=)(C)] = 0 (.)
se denomina EDO de orden = y representa una relación entre la variable independiente , C , y la variable dependiente , la función H , junto con sus derivadas H^0 , H^00 ,... , H(=). Otra manera de expresar (.) es mediante su forma normal o estándar , i.e., H(=)^ = 5 (C, H, H^0 , H^00 ,... , H(= ^1 )) (.)
donde se asume que 5 es una función real.
Usualmente, interesa determinar la solución de la ecuación diferencial (ordinaria o en derivadas parciales) que modela el fenómeno bajo estudio, dado que provee de información relevante acerca de este, para lo cual en primer lugar hay que definir qué se entiende por solución de una EDO.
Definition .. Solución de una EDO Una solución de la EDO (.) sobre el intervalo C 2 ( 0 , 1 ) es una función H(C) = )(C) tal que )^0 , )^00 , )^000 ,... , )(=)^ existen y satisfacen
)(=)^ = 5 (C, ), )^0 , )^00 ,... , )(= ^1 )), 8 C 2 ( 0 , 1 )
(^) A partir de ahora sólo se hablará de EDOs y los métodos numéricos involucrados en su resolución. Las EDPs tienen su tema correspondiente más adelante.
primer orden Problema de Valor Inicial
valor inicial (PVI) en el que la EDO involucrada es de orden , i.e.,
3H(C) 3C
= 5 (C, H), 0 < C < 1
H(C 0 ) = H 0
La razón de estudiar sólo este tipo de EDOs se debe a que, por un pro- cedimiento denominado reducción de orden , toda EDO de orden superior a puede expresarse como un sistema de EDOs de orden . A modo de ejemplo, la EDO de orden ,
3 2 G(C)
3C 2
3G(C)
3C
+ :G(C) = 0
que modela el movimiento de un sistema masa-resorte (<, :) con amor- tigüación ( 2 ), puede escribirse como el siguiente sistema de EDOs de orden (^) ( 3G(C) 3C =^ H(C) 3H(C) 3C =^ ^
: < G(C)^ ^
2 < H(C)
. Métodos numéricos de resolución para PVIs
Para el PVI,
3H(C) 3C
= 5 (C, H), 0 C 1 (.)
H(C 0 ) = H 0
y suponiendo que se cumplen las condiciones de existencia y unicidad de la solución ( Teorema .. ), la resolución analítica sólo es posible cuando la EDO involucrada pertenece a una de estas categorías:
. Lineal. . Separable. . Exacta. . Ninguna de las anteriores, pero puede transformarse en una de ellas.
En el resto de casos, no hay método analítocos para obtener la solución del problema y, por tanto, se recurrirán a los métodos numéricos.
En este sentido, pueden identificarse dos categorías de métodos numéri- cos orientados a la resolución de (.):
. Métodos de un paso. Este tipo de métodos permiten el cálculo de la aproximación H (^8) + 1 a partir de la aproximación anterior H 8. A este tipo pertenecen los llamados métodos o técnicas Runge-Kutta. Entre ellos cabe mencionar: a) El método de Euler ‡^. b) Los métodos de Heun y del punto medio. c) El método de Runge-Kutta IV.
‡ (^) No pertenece estrictamente a los métodos Runge-Kutta, aunque puede pensarse como una versión cero de los mismos.
primer orden Problema de Valor Inicial
. Métodos multipaso. Estos métodos calculan la aproximación H (^8) + 1 a partir de varias aproximaciones anteriores.
Errores asociados a los métodos de aproximación de
ODEs
Previo al estudio de los métodos numéricos que permiten aproximar la solución de un PVI, es necesario introducir los errores asociados a su implementación:
Error de truncamiento o discretización Este error es inherente a los métodos numéricos empleados para aproximar los valores de la solución real. Este a su vez, está compuesto de los siguientes errores: I Error de truncamiento local, ⇢ (^) ; Este error depende solamente del método numérico empleado y aparece al aplicar este en un sólo paso. Para conocer su magnitud y propiedades se recurre al desar- rollo de Taylor de la solución del PVI: Sea H(C) la solución real de la EDO 3H3C = 5 (C, H) asociada al PVI (.). Suponiendo que tiene derivadas continuas, su desarrollo de Taylor de orden = en un entorno del punto (C 8 , H 8 ) resulta
H(C 8 + 1 ) = H(C 8 +⌘) = H(C 8 )+H 0 (C 8 )⌘+
H 00 (C 8 )
H (=)^ (C 8 )
⌘ =^ +
H (=+^1 )^ (⇢)
⌘ =+^1 , ⇢ 2 [C 8 , C 8 + 1 ]
donde el término en rojo es el resto asociado a dicho desarrollo El símbolo “ O () ” significa “ términos de orden ().
. De hecho, el resto constituye una cota superior al error de truncamiento local, i.e.,
⌘ =+^1 , " max | 5 (=)^ (⇢, H(⇢))|, ◆ 2 (G 8 , G (^8) + 1 )
Una forma alternativa de escribir (.) involucra a la función 5 de la EDO, i.e.,
H(C 8 + 1 ) = H(C 8 +⌘) = H(C 8 )+ 5 (C 8 , H 8 )⌘+
5 0 (C 8 , H 8 )
5 (= ^1 )^ (C 8 , H 8 )
⌘ =^ +$(⌘ =+^1 )
donde $(⌘ =+^1 ) establece que el error local de trun- camiento es proporcional a ⌘ =+^1. Como se comprobará en cada uno de los métodos, el error de truncamiento aparece al aproximar la solución verdadera haciendo uso de un número finito de términos en el desarrollo de Taylor de la solución real, esto es, se trunca u omite parte de esta solución.
Si el error de truncamiento local tiende a 0 cuando ⌘! 0 , se dice que el método en cuestión es consistente. Específicamente, si ⇢ (^) ; = $(⌘?^ ) cuando ⌘! 0 , se dice que el método es consistente de orden?.
primer orden Problema de Valor Inicial
I Para obtener la siguiente aproximación, H 2 , se repite de nuevo el proceso. Como no se conoce el valor de la solución real en C 1 , )(C 1 ), se hace uso del valor aproximado H 1. Así, en el punto (C 1 , H 1 ) se traza la recta tangente con pendiente 5 (C 1 , H 1 ), i.e.,
H = H 1 + 5 (C 1 , H 1 )(C C 1 )
de forma que la aproximación al valor )(C 2 ) en torno al C = C (^2) resulta H 2 = H 1 + 5 (C 1 , H 1 )(C 2 C 1 ) I Procediendo de forma similar, se emplea el valor de cada aproxi- mación H en cada iteración para calcular la pendiente de la recta tangente con la que se calcula la aproximación siguiente. La expresión general para la aproximación H (^8) + 1 en términos de C 8 , C (^8) + 1 e H 8 es H (^8) + 1 = H 8 + 5 (C 8 , H 8 )(C (^8) + 1 C 8 ) (.) I Asumiendo que la longitud de paso , ⌘ , es uniforme entre los puntos C 0 , C 1 , C 2 ,... , i.e.,
C (^8) + 1 = C 8 + ⌘, 8 = 0 , 1 , 2 ,...
la fórmula de Euler puede reescribirse de la forma
H (^8) + 1 = H 8 + ⌘ 5 (C 8 , H 8 ), 8 = 0 , 1 , 2 ,... (.)
Así, mediante el método de Euler se genera una sucesión de valores {H 0 , H 1 , H 2 ,.. .} que aproximan los valores de la solución real )(C) en los puntos {C 0 , C 1 , C 2 ,.. .}. Desde un punto de vista geométrico, la gráfica de la solución real, H = )(C), se aproxima mediante la gráfica de una función lineal a trozos, donde cada tramo es una recta tangente.
Ejemplo .. Aplicación del Método de Euler
Sea el PVI (^) ( 3H 3G =^4 +^4
C 1
3 H,
H( 0 ) = 1
Hacer uso del método de Euler con cuatro nodos para aproximar los valores de la solución en el intervalo [ 0 , 0. 4 ]. Comparar el resultado con los valores de la solución real.
La solución exacta del PVI es
H(C) = )(C) =
4 C^ ( 3 19 4
(^23) C
. Dado el número de nodos, # = 4 , y el intervalo donde se quiere aproximar la solución [ 0 , 0. 4 ] , se calcula la longitud de paso, ⌘ , como ⌘ =
Luego, los nodos en los que se van a calcular las aproximaciones
primer orden Problema de Valor Inicial
Fig. .. ⌘ = 0. 1
Fig. .. ⌘ = 0. 05
son C 0 = 0 , C 1 = 0. 1 , C 2 = 0. 2 , C 3 = 0. 3 , C 4 = 0. 4
. Conocidos los nodos y la longitud de paso, procedemos a la aplicación de la fórmula de Euler para calcular cada una de las aproximaciones H 8 , 8 = 1 , 2 , 3 , 4 , donde 5 (C, H) = 4 + 4 C^ 13 H , junto con el valor de la solución exacta, H(C), y el error global, ⇢ 8 = |H 8 H(C 8 )|, 8 = 1 , 2 , 3 , 4 en cada nodo:
(# = 1 ) ) H 1 = H 0 + ⌘ 5 (C 0 , H 0 ) = 0 + 0. 1 5 ( 0 , 1 ) = 1. 466667 H(C 1 ) = )( 0. 1 ) = 1. 454191 ⇢ 6 ( 1 ) = H(C 1 ) H( 1 ) = 0. 012476 (# = 2 ) ) H 2 = H 1 + ⌘ 5 (C 1 , H 1 ) = 1. 466667 + 0. 1 5 ( 0. 1 , 1. 466667 ) = 1. 908262 H(C 2 ) = )( 0. 2 ) = 1. 884588 ⇢ 6 ( 2 ) = H(C 2 ) H 2 = 0. 023674 (# = 3 ) ) H 3 = H 2 + ⌘ 5 (C 2 , H 2 ) = 1. 908262 + 0. 1 5 ( 0. 2 , 1. 908262 ) = 2. 326526 H(C 3 ) = )( 0. 3 ) = 2. 292817 ⇢ 6 ( 3 ) = H(C 3 ) H 3 = 0. 033709 (# = 4 ) ) H 4 = H 3 + ⌘ 5 (C 3 , H 3 ) = 2. 326526 + 0. 1 5 ( 0. 3 , 2. 326526 ) = 2. 723057 H(C 4 ) = )( 0. 4 ) = 2. 680373 ⇢ 6 ( 4 ) = H(C 4 ) H 4 = 0. 353847
A partir de las figuras . y . del ejemplo, puede observarse lo sigu- iente:
. La solución numérica sigue la tendencia general de la solución exacta de forma correcta. . La solución numérica asociada a una longitud de paso más pequeña es más precisa que la correspondiente a un ⌘ más grande. . Como ⌘ 1 = 0. 1 y ⌘ 2 = 0. 05 si, por ejemplo, se calcula el cociente de los errores globales en C 4 = 0. 4 se observa que
⇢ 6 (⌘ = 0. 1 ) ⇢ 6 (⌘ = 0. 05 )
lo cual muestra que el método es de primer orden ( ⇢ 6 = $(⌘)). . Los errores en general son demasiado grandes. Esto se debe a que el error global es del orden de ⌘. . Los errores son todos negativos. Esto indica que la solución numérica “sigue” a la solución exacta. Esto ocurre dado que la primera derivada de 5 (C, H) decrece a medida que se C aumenta. . La solución numérica se acerca mejor a la solución exacta cuando ⌘ disminuye. Esta propiedad se conoce como convergencia del método.
Análisis del error en el método de Euler
Comparando las ecuaciones (.) y (.) se observa que el método de Euler se corresponde con el desarrollo de Taylor hasta el término lineal 5 (C 8 , H 8 )⌘. Así, el error de truncamiento local se debe a los términos omitidos de
primer orden Problema de Valor Inicial
Pseudocódigo del Método de Euler
A continuación se muestra el pseudocódigo asociado al método de Euler para su implementación en MATLAB.
Algoritmo : Pseudo-código del método de Euler input : f, ytrue, a, b, y, N h = (b - a)/N ; /* longitud de paso (^) / t = a:h:b ; / inicializacion de t (^) / y = zeros(N+,) ; / inicializacion de y (^) / y() = y; errg = zeros(N+,); / inicializacion del error (^) / for 8 = 1 , 2 ,... , N do k = f(t(i),y(i)); y(i+) = y(i) + kh; errg = ytrue(t(i)) - y(i); end output : x, y, ytrue, errg plot(t,y); /* representacion grafica (^) */ hold on; fplot(@(t) ytrue(t),[a b]);
Comentarios:
. En la práctica, se escogerá un valor la longitud de paso, ⌘ , su- ficientemente pequeño como para alcanzar una cierta precisión, pero no demasiado pequeño. Esto se debe a que una longitud de paso demasiado pequeña aumenta el número de cálculos y los enletece, de forma que se incrementa el coste computacional, lo que en muchos casos supone una pérdida en la precisión de los resultados. . Dado que, de entre todos los errores, el error de truncamiento global es el más accesible, será el que se calcula en cada uno de los métodos numéricos de este tema.
. Métodos de Heun y del punto medio
Dos de las principales desventajas del método de Euler son:
I Requiere una longitud de paso muy pequeña para producir resul- tados suficientemente precisos. I Emplea la misma derivada a lo largo de todo el intervalo de integración. Ambos “handicaps” puede resolverse mediante dos sencillas modi- ficaciones. Ambas pertenecen a un conjunto más amplio de técnicas numéricas conocidas como métodos de Runge-Kutta.
Método de Heun
El método de Heun o método de Euler mejorado perfecciona la esti- mación de la pendiente mediante el cálculo de dos derivadas: una en el punto inicial, C 8 , y otra en el punto final, C (^8) + 1. Después, considera el promedio de ambas para conseguir una estimación mejorada de la
primer orden Problema de Valor Inicial
Fig. .. Representación gráfica del método de Heun
pendiente a lo largo de todo el intervalo de integración, [C 8 , C (^8) + 1 ]. Los pasos son los siguientes:
. La pendiente al comienzo del intervalo es
H
0 8 =^5 (C^8 ,^ H^8 )
que se emplea para calcular la siguiente aproximación
H (^08) + 1 = H 8 + 5 (C 8 , H 8 )⌘ (.)
Si quisiera aplicarse el método de Euler, nos detendríamos aquí. Sin embargo, en el método de Heun, H (^08) + 1 es una estimación intermedia, de ahí el empleo del superíndice 0. . La ecuación (.) a veces se conoce como la ecuación predictora , dado que provee de una primera estimación de H (^8) + 1 que permite el cálculo final de la pendiente al final del intervalo, en G (^8) + 1 ,
H
0 8 + 1 =^5 (C^8 +^1 ,^ H^
0 8 + 1 )^ (.)
. Combinando las pendientes (.) y (.) se calcula una pendiente promedio
H ¯
0
H
0 8 +^ H^
0 8 + 1 2
5 (C 8 , H 8 ) + 5 (C 8 + 1 , H 08 + 1 )
. La pendiente (.) es la que se emplea para calcular la siguiente aproximación H (^8) + 1 a partir de H 8 del mismo modo que en el método de Euler,
H 8 + 1 = H 8 + ¯H
0 = H 8 +
5 (C 8 , H 8 ) + 5 (C 8 + 1 , H 80 + 1 )
donde (.) recibe el nombre de ecuación correctora.
Como puede observarse, el método de Heun es un método que consta de dos pasos: en primer lugar, calcula H (^08) + 1 mediante la fórmula de Euler; en segundo lugar calcula H (^8) + 1 mediante (.).
Ejemplo .. Aplicación del Método de Heun Hacer uso del método de Heun para aproximar la solución del PVI, ( (^) 3H 3C =^ HC^
2 1. 1 H
H( 0 ) = 1
en el intervalo [ 0 , 1 ] con # = 4.
La solución exacta del PVI es:
H(C) = 4 ^
11 10 C+^ 1 3 C^ 3 (.)
Por otro lado, dado el intervalo de integración y el número de nodos
(que no incluye al primero) calculamos la longitud de paso, ⌘ ,
primer orden Problema de Valor Inicial
Fig. .. Representación gráfica del método del punto medio
[C 8 , C (^8) + 1 ]. Esto resulta importante dado que, típicamente, casi todo el esfuerzo computacional en cada paso se emplea en la evaluación de 5. . En comparación con el método de Euler, el método de Heun es más eficiente§^ , proporcionando mejores aproximaciones o requiriendo menor esfuerzo computacional, o ambos.
A continuación se muestra el pseudo-código asociado al método de Heun:
Algoritmo : Pseudo-código del método de Heun input : f, ytrue, a, b, y, N h = (b - a)/N ; /* longitud de paso (^) / T = a:h:b ; / inicializacion de t (^) / Y = zeros(N+,); Y() = ; Y = zeros(N+,) ; / inicializacion de y (^) / Y() = y; errg = zeros(N+,); / inicializacion del error (^) / for 8 = 1 , 2 ,... , N do k = f(T(i),Y(i)); Y(i+) = Y(i)+hk; k = f(T(i+),Y(i+)); Y(i+) = Y(i) + .h(k+k); errg = ytrue(T(i+)) - Y(i+);
end output : T, Y, ytrue, errg plot(T,Y); /* representacion grafica (^) */ hold on; fplot(@(t) ytrue(t),[a b]);
Método del punto medio
Otra simple modificación en el método de Euler da lugar a lo que se conoce como el método del punto medio o método de la poligonal mejorado. Esta técnica emplea el método de Euler para predecir el valor de H en el punto medio del intervalo. Así:
. Calcula H (^8) + 1 2 , la aproximación en C 8 + ⌘ 2 (punto medio del intervalo),
H 8 + 12 = H 8 + 5 (C 8 , H 8 )
. Emplea el valor en (.) para calcular la pendiente en el punto medio, H
0 8 + 12 =^5
C 8 + 12 , H 8 + 12
que se asume como una correcta aproximación de la pendiente promedio en todo el intervalo [C 8 , C (^8) + 1 ]. . La pendiente en (.) se utiliza para calcular la aproximación en el nodo C (^8) + 1 , H (^8) + 1 = H 8 + 5
C 8 + 1
2
, H 8 + 1
2
§ (^) El método de Heun con una longitud de paso, ⌘ , requiere el mismo número de evaluaciones de 5 que el método de Euler con una longitud de paso ⌘ 2.
primer orden Problema de Valor Inicial
Ejemplo .. Aplicación del método del punto medio Emplear el método del punto medio para aproximar la solución del PVI (^) ( 3H 3C =^
H 2 + 2 C H 3 +C 2 H( 0 ) = 0. 5
en el intervalo [ 0 , 1 ] y # = 4.
La solución exacta del PVI es
H(C) =
3 C 2
C 6
Como el intervalo de integración es [ 0 , 1 ] y el número de nodos es
= 4 , la longitud de paso resulta
Luego, el intervalo [ 0 , 1 ] queda discretizado del siguiente modo,
{C 0 = 0 , C 1 = 0. 25 , C 2 = 0. 50 , C 3 = 0. 75 , C 4 = 1 }
y, por tanto, las aproximaciones a calcular son las correspondientes a
estos nodos con 5 (C, H) = H 2 + 2 C H 3 C 2. Así:
(# = 1 ) ) H 1
2
= H 0 +
5 (C 0 , H 0 ) = 0 + 0. 125 5 ( 0 , 0. 5 ) = 0. 510417
H 1 = H 0 + 0. 25 5 (C 12 , H 12 ) = 0. 532177
H(C 1 ) = )( 0. 25 ) = 0. 532609
⇢ 6 ( 1 ) = H(C 1 ) H 1 = 0. 00432
(# = 2 ) ) H 32 = H 1 +
5 (C 1 , H 1 )
H 2 = H 1 + 0. 25 (C 32 , H 32 ) = 0. 589771
H(C 2 ) = )( 0. 50 ) = 0. 590909
⇢ 6 ( 2 ) = H(C 2 ) H 2 = 0. 001138
(# = 3 ) ) H 5
2
= H 2 +
5 (C 2 , H 2 )
H 3 = H 2 + 0. 25 5 (C 52 , H 52 ) = 0. 676330
H(C 3 ) = )( 0. 75 ) = 0. 678571
⇢ 6 ( 3 ) = H(C 3 ) H 3 = 0. 002242
(# = 4 ) ) H 72 = H 3 +
5 (C 3 , H 3 )
H 4 = H 3 + 0. 25 5 (C 7
2
, H 7
2
H(C 4 ) = )( 1 ) = 0. 8
⇢ 6 ( 4 ) = H(C 4 ) H 4 = 0. 003909
primer orden Problema de Valor Inicial
. Método de Runge-Kutta IV
Existen numerosas versiones de los métodos de Runge-Kutta (RK) , sin embargo, todas ellas pueden generalizarse mediante la ecuación
H (^8) + 1 = H 8 + "(C 8 , H 8 , ⌘)⌘ (.)
donde "(C 8 , H 8 , ⌘) se denomina función de incremento , que puede inter- pretarse como la pendiente promedio a lo largo del intervalo [C 8 , C (^8) + 1 ]. La forma general de esta función es
" = 0 1 : 1 + 0 2 : 2 +... + (^0) = : (^) = (.)
donde los coeficientes { 0 8 } = 8 = 1 son constantes reales y los parámetros {: 8 } = 8 = 1 son las pendientes en diferentes puntos del intervalo, i.e.,
: 1 = 5 (C 8 , H 8 )
: 2 = 5 (C 8 +? 1 ⌘, H 8 + @ 11 : 1 ⌘)
: 3 = 5 (C 8 +? 2 ⌘, H 8 + @ 21 : 1 ⌘ + @ 22 : 2 ⌘)
: 4 = 5 (C 8 +? 3 ⌘, H 8 + @ 31 : 1 ⌘ + @ 32 : 2 ⌘ + @ 33 : 3 ⌘)
: = = 5 (C 8 +? = 1 ⌘, H 8 + @ = 1 , 1 : 1 ⌘ + @ = 1 , 2 : 2 ⌘, +... + @ = 1 ,= 1 : = 1 ⌘)
donde los coeficientes {? 8 } = 8 = 11 y {@ (^) = 1 , 9 } = 9 = 11 son constantes reales.
Es importante resaltar que las diferentes pendientes, : 8 , están definidas de forma recursiva: : 1 aparece en la ecuación de : 2 , que aparece en la ecuación de : 3 , y así sucesivamente. Además, como cada : consiste en una evaluación de la función 5 (C, H), dicha recurrencia hace que los métodos RK sean eficientes para los cálculos que ha de realizar el ordenador.
A partir de la expresión general, (.), pueden derivarse varias técnicas Runge-Kutta dando diferentes valores a los términos involucrados en la función incremento en (.).
I El método de Runge-Kutta de orden 1 ( = = 1 ) se corresponde con el método de Euler con
0 1 = 1 ,? 0 = 0 y @ (^0) , 0 = 0
de forma que H (^8) + 1 = H 8 + 5 (C 8 , H 8 )⌘
I El método de Runge-Kutta de orden 2 ( = = 2 ) con
se corresponde con el método de Heun,
H 8 + 1 = H 8 +
primer orden Problema de Valor Inicial
Fig. .. Representación gráfica del método RK IV
donde
: 1 = 5 (C 8 , H 8 ) : 2 = 5 (C 8 + ⌘, H 8 + : 1 ⌘)
I El método de Runge-Kutta de orden 2 ( = = 2 ) con
se corresponde con el método del punto medio,
H (^8) + 1 = H 8 + : 2 ⌘
donde
: 1 = 5 (C 8 , H 8 )
: 2 = 5
C 8 +
⌘, H 8 +
I El método de Runge-Kutta de orden 2 ( = = 2 ) con
se corresponde con el llamado método de Ralston ¶^ ,
H 8 + 1 = H 8 +
donde
: 1 = 5 (C 8 , H 8 )
: 2 = 5
C 8 +
⌘, H 8 +
I De entre todos, el más popular de los métodos RK es el de cuarto orden (= = 4 ). Al igual que para los métodos RK de orden , existen numerosas versiones de este método, sin embargo, la que se mostrará a continuación es la más comúnmente usada, y por ello recibe el nombre de método RK de cuarto orden ,
H 8 + 1 = H 8 +
donde
: 1 = 5 (C 8 , H 8 )
: 2 = 5
C 8 +
⌘, H 8 +
C 8 +
⌘, H 8 +
: 4 = 5 (C 8 + ⌘, H 8 + : 3 ⌘)
Puede observarse a partir de (.) que el método de RK de cuarto
¶ (^) El método de Ralston fue propuesto por los matemáticos Anthony Ralston y Philip Rabinowitz en .
primer orden Problema de Valor Inicial
= 4 , la longitud de paso resulta
Luego, el intervalo [ 0 , 2 ] queda discretizado del siguiente modo,
{C 0 = 0 , C 1 = 0. 5 , C 2 = 1. 5 , C 3 = 1. 5 , C 4 = 2 }
y, por tanto, las aproximaciones a calcular son las correspondientes a
primer orden Problema de Valor Inicial
Fig. .. ⌘ = 0. 25
Fig. .. ⌘ = 0. 125
estos nodos con 5 (C, H) = 12 C + 2 H. Así:
( 1 ) 1 =^5 (C^0 ,^ H^0 )^ =^2.^5
: 2 (^1 )= 5
C 0 +
⌘, H 0 +
: 1 (^1 )⌘
: 3 (^1 )= 5
C 0 +
⌘, H 0 +
: 2 (^1 )⌘
: 4 (^1 )= 5
C 0 + ⌘, H 0 + : ( 31 )⌘
H 1 = H 0 +
( 1 ) 1 +^2 :^
( 1 ) 2 +^2 :^
( 1 ) 3 +^ :^
( 1 ) 4 )⌘^ =^2.^958333 H(C 1 ) = )( 0. 5 ) = 2. 968282 ⇢ 6 ( 1 ) = H(C 1 ) H 1 = 0. 009948
(# = 2 ) ) : 1 (^2 )= 5 (C 1 , H 1 ) = 5. 916667
: 2 (^2 )= 5
C 1 +
⌘, H 1 +
: 1 (^2 )⌘
: 3 (^2 )= 5
C 1 +
⌘, H 1 +
: 2 (^2 )⌘
( 2 ) 4 =^5
C 1 + ⌘, H 1 + :
( 2 ) 3 ⌘
H 2 = H 1 +
(: 1 (^2 )+ 2 : ( 2 2 )+ 2 : 3 (^2 )+ : ( 42 ))⌘ = 7. 835069
H(C 2 ) = )( 1. 0 ) = 7. 889056
⇢ 6 ( 2 ) = H(C 2 ) H 2 = 0. 053987
(# = 3 ) ) : 1 (^3 )= 5 (C 2 , H 2 ) = 15. 170139
: 2 (^3 )= 5
C 2 +
⌘, H 2 +
: 1 (^3 )⌘
( 3 ) 3 =^5
C 2 +
⌘, H 2 +
( 3 ) 2 ⌘
: 4 (^3 )= 5
C 2 + ⌘, H 2 + : ( 33 )⌘
H 3 = H 2 +
(: 1 (^3 )+ 2 : ( 2 3 )+ 2 : 3 (^3 )+ : ( 43 ))⌘ = 20. 615813
H(C 3 ) = )( 1. 5 ) = 20. 835537
⇢ 6 ( 3 ) = H(C 3 ) H 3 =
(# = 4 ) ) : 1 (^4 )= 5 (C 3 , H 3 ) = 40. 231626
( 4 ) 2 =^5
C 3 +
⌘, H 3 +
( 4 ) 1 ⌘
: 3 (^4 )= 5
C 3 +
⌘, H 3 +
: 2 (^4 )⌘
: 4 (^4 )= 5
C 3 + ⌘, H 3 + : ( 34 )⌘
H 1 = H 3 +
(: 1 (^4 )+ 2 : ( 2 4 )+ 2 : 3 (^4 )+ : ( 44 ))⌘ = 54. 803244
H(C 4 ) = )( 2 ) = 55. 598150
⇢ 6 ( 4 ) = H(C 4 ) H 4 = 0. 794906
A partir de las figuras . y . del ejemplo, puede observarse lo siguiente:
