Solución numérica de ecuaciones diferenciales

Solución numérica de ecuaciones diferenciales (I)
Las leyes que gobiernan los fenómenos de la naturaleza se expresan habitualmente en forma de ecuaciones diferenciales. Las ecuaciones del movimiento de los cuerpos (la segunda ley de Newton) es una ecuación diferencial de segundo orden, como lo es la ecuación que describe los sistemas oscilantes, la propagación de las ondas, la transmisión delcalor, la difusión, el movimiento de partículas subatómicas, etc.
Pocas ecuaciones diferenciales tienen una solución analítica sencilla, la mayor parte de las veces es necesario realizar aproximaciones, estudiar el comportamiento del sistema bajo ciertas condiciones. Así, en un sistema tan simple como un péndulo, la amplitud de la oscilación ha de ser pequeña y el rozamiento ha de serdespreciable, para obtener una solución sencilla que describa aproximadamente su movimiento periódico.
Se estudia el procedimiento de Runge-Kutta que se aplica de forma directa a una ecuación diferencial de primer orden, pero veremos como se extiende a un sistema de ecuaciones de primer orden, a un ecuación diferencial de segundo orden y a un sistema de ecuaciones diferenciales de segundo orden.
Elprocedimiento de Runge-Kutta se puede programar fácilmente en los ordenadores y además, se emplea mucho en la práctica, debido a la su exactitud relativamente elevada de la solución aproximada de la ecuación diferencial. La justificación del procedimiento de Runge-Kutta no es sencilla, el lector interesado puede consultar algún libro de métodos numéricos de análisis.
Método de Euler
Vamos a resolverla ecuación diferencial de primer orden

Con la condición inicial de que en el instante t0 la posición es x0






La primera derivada nos permite conocer la posición xi+1 en el instante ti+1, a partir de la posición xi en el instante ti de acuerdo a la fórmula siguiente. La línea de color rojo es la tangente a la curva en el instante ti

El procedimiento de Euler produce un error que seacumula a cada paso h de integración, que es el segmento en color azul que une los dos puntos en la figura.
Escribimos una función denominada euler, a la que le pasaremos:
la función f(t,x),
la condición inicial de que en el instante t0 la posición es x0,
el instante final tf
el número de pasos de integración n
y nos devolverá un vector t y su correspondiente vector x.
Supongamos quequeremos integrar la ecuación diferencial
dx/dt= cos t
con las condición inicial t=0, x=0.

Tomamos un intervalo h=π/6, y construimos la siguiente tabla
t
dxdt=cost
x(Euler)
x=sin t
0
1
0
0
π/6
0.866
0.523
0.5
π/3
0.5
0.977
0.866
π/2
0
1.239
1
2π/3
-0.5
1.239
0.866
5π/6
-0.866
0.977
0.5
π
 
0.523
0
Esta tabla nos ilustra el modo de aplicar el método de Euler a unaecuación diferencial de primer orden. Para aplicar el método de Euler precisamos de un paso h pequeño, incluso así los errores se van acumulando y al cabo de cierto tiempo la diferencia entre el valor exacto y el calculado es grande.

Hay diferencia entre la solución exacta y la obtenida mediante integración numérica por el método de Euler
Método de Runge-Kutta
En esta sección vamos a estudiarla aplicación del método de Runge-Kutta a:
Una ecuación diferencial de primer orden
Un sistema de dos ecuaciones diferenciales de primer orden
Una ecuación difrencial de segundo orden
Un sistema de dos ecuaciones diferenciales de segundo orden
Ecuación diferencial de primer orden
Sea una ecuación diferencial de primer orden, con la condición inicial de que en el instante t0 el valor inicialde x es x0
Se elige una anchura de paso h y se calculan cuatro números k1, k2, k3, k4 de acuerdo con el procedimiento esquematizado en la tabla adjunta. Según el procedimiento ordinario de Runge-Kutta, a partir del valor de x en el instantet se determina el valor de x en el instante t+h mediante la fórmula que figura en la última fila de dicha tabla.

Definimos la función rk_1 que resuelve...