Solución numérica de las ecuaciones de movimento, varios métodos Runge-Kutta

En este capítulo vamos a discutir el método usado en nuestro trabajo para resolver las ecuaciones de movimiento. Hay muchos métodos numéricos para integrar ecuaciones diferenciales y su aplicabilidad depende del tipo de ecuación que se esté tratando. Uno de los métodos más simples y que es muy confiable, es el de Runge-Kutta en varias versiones como explicaremos a continuación.

2.1 Supóngase una ecuación (trataremos por el momento con una sola variable) de la forma

$$\frac{dy}{dt} = F(y,t)$$ (II-1-1)

donde $F$ es del tipo

$$\begin{eqnarray*} F(y,t) & = & m(t)\cdot y +g(t) \end{eqnarray*}$$

La gráfica de la función se representará como en la Figura 2.1. Las derivadas de dicha función se conocen en todo punto, ya que vienen especificadas por la ecuación diferencial. Si $h$ es el incremento de tiempo, entonces, a primera aproximación, se tendrá como valor de la función después de un tiempo $h$, la cantidad

$$y(t_0+h) = y_0 + y^{\prime}\cdot h$$ (II-1-2)

esto puede representarse en la Figura 2.1

Figure 2.1: derivada

En el caso no lineal que es el que nos interesa se escribe

$$y(t_0+h) = y(t_0) + h\cdot F(y,t)$$ (II-1-3)

aunque eso no es precisamente correcto debido a la no linealidad. Como se ve, este método consiste en tomar una nueva derivada después de cierto tiempo, y así sucesivamente, porque la tangente a la curva está cambiando constantemente su dirección. Puede entenderse a la función $F(y,t)$ como una superficie en tres dimensiones de la cual nos trasladamos a un campo de tangentes como indica el siguiente dibujo.

Figure 2.2: flechas

En cada punto habrá una flecha indicando en que dirección hay que moverse. Con este procedimiento no se sigue exactamente a la curva pero se puede confiar en que escogiendo $h$ suficientemente pequeño habrá un error $\varepsilon$, de tal manera que el resultado final diferirá del valor real por una cantidad del orden de $\varepsilon$.

El método que acabamos de exponer se conoce como ``Método de Euler'' y no es nuestro propósito discutirlo aquí con todo su detalle. La técnica puede mejorarse si en vez de una recta se toma por ejemplo una parábola tangente a cada punto de la curva.

2.2 Considere ahora la ecuación $y^{\prime}=-y$, si se toma $y=1$ en $t=0$, al desarrollar la función en serie de potencias alrededor de $t_0=0$ se tiene:

$$\begin{eqnarray*} y & = & y_0+hy^{\prime}_0 + \frac{h^2}{2!}y^{\prime\prime}... ... = 1 - h + \frac{h^2}{2!} - \frac{h^3}{3!} + \cdots = e^{-h} \end{eqnarray*}$$

eso debido a que

$$\begin{eqnarray*} y^{\prime}_0 & = & -1 \\ y^{\prime\prime}_0 & = & (y^{\pr... ..._0^{\prime\prime})^{\prime} = -1\\ \cdots & \cdots & \cdots \end{eqnarray*}$$

No hay ningún argumento que impida desarrollar ``$y$'' en serie de potencias pero si se toma la función en toda su generalidad se tendrá

$$\begin{eqnarray*} y^{\prime} & = & F(y,t) \\ y^{\prime\prime} & = & \frac{d... ... y^{\prime\prime\prime} & = & \frac{d}{dt}(F_y\cdot F + F_t) \end{eqnarray*}$$

Los subíndices significan derivación parcial con respecto a la variable indicada; entonces, para $F^{\prime\prime}$ se tendra

$$\begin{eqnarray*} y^{\prime\prime\prime} & = & \frac{d}{dt}F_y\cdot F + F_y\... ... F_y\cdot F^{\prime} + \frac{\partial}{\partial t}F^{\prime} \end{eqnarray*}$$

Se cambió el órden de la derivación. Sustituyendo $F^{\prime}$ resulta lo siguiente:

$$\begin{eqnarray*} y^{\prime\prime\prime} & = & \frac{\partial}{\partial y} (... ...\cdot F+F_t) + \frac{\partial}{\partial t}(F_y\cdot F + F_t) \end{eqnarray*}$$

Resumiendo, se tiene para las tres primeras derivadas:

$$\begin{array}{ccl} y^{\prime} & = & F(y,t) \\ y^{\prim... ...& = & F_{yy}F^2+2F_y^2F+2F_{yt}F+2F_yF_t+F_{tt} \end{array}$$ (II-2-1)

Por el mismo procedimiento pueden ser encontradas las derivadas de orden superior. Simbólicamente se resuelve la ecuación diferencial tomando todas las derivadas formando así una serie de Taylor lo que no es práctico porque hay que efectuar una cantidad de cálculos muy grande y tediosa, y eso para cada ecuación diferencial que se presente.

2.3 Nuestro deseo es conocer una técnica más general que permita resolver cualquier ecuación diferencial, sin necesidad de conocer su forma precisa, para obtener las derivadas; el método que emplearemos permite calcularlas de otra manera. La base de este concepto es que existen muchas funciones lineales, por ejemplo, los polinomios forman un espacio vectorial (espacio dual) y por lo tanto puede formarse una base en ese espacio de manera que toda función analítica puede expresarse como una combinación lineal de elementos de dicha base. Los polinomios tiene propiedades muy interesantes, por ejemplo, en un polinomio de segundo grado basta conocer tres valores diferentes de la función para conocer su comportamiento en cualquier otro punto y de esa manera puede encontrarse su derivada como combinación de los tres valores conocidos. Es bien sabido que tres puntos determinan una parábola, de la manera siguiente:

$$y = y_1\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)} + y_2\f... ...1)(x_2-x_3)} + y_3\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}$$

Esta expresión puede también escribirse como un determinante:

$$\left\vert \begin{array}{cccc} 1 & 1 & 1 & 1 \\ x & x_1 ... ...& x_3^2 \\ y & y_1 & y_2 & y_3 \end{array} \right\vert = 0$$

Esta situación se representa en la Figura 2.3:

Figure 2.3: parábola

Al hacer operaciones se obtiene un polinomio de segundo grado, así, puede calcularse fácilmente $dy/dt$ porque es una suma de productos y diferencias. Con el objeto de evitar una gran cantidad de pasos algebraicos se procede en una forma equivalente pero más sencilla como veremos en seguida:

Figure 2.4: curva

Se toma un punto, definido por $(y_0+\lambda y^{\prime}, t_0+\mu h)$, donde la función expresa la derivada de cualquier curva que pase por ese punto; posiblemente no la curva que nos interesa sino alguna que pase precisamente por ahí con esa derivada, Figura 2.4. Esperamos obtener las derivadas de F en diversos de tales puntos. Dibujemos esta función en tres dimensiones, representando un cambio como el que acabamos de mencionar, Figura 2.5.

Figure 2.5: superficie

Desarróllese la función en serie de potencias alrededor del punto escogido de la manera anterior y utilizando los parámetros indicados, eso es

$$\begin{array}{ccl} F(y_0+\lambda hF_0,t_0+\mu h) & = & F(... ...{\partial^2 F}{\partial t^2} + \cdots \right) \end{array}$$ (II-3-1)

La pendiente en este punto, entonces, queda expresada en términos de la pendiente en el punto original, teniendo en cuenta que F como función de dos variables puede desarrollarse en serie de potencias, como acaba de hacerse. Esta ecuación contiene las derivadas deseadas; la idea consiste en invertirla. Esto se hace por un procedimiento más sencillo, por ejemplo, escribiendo $y$ como:

$$y = y_0+hy^{\prime}+\frac{h^2}{2}y^{\prime\prime} + \cdots$$ (II-3-2)

Se considerará la hipótesis de suponer para $y$ un desarrollo de la forma:

$$y = y_0+\alpha_0\kappa_0+\alpha_1\kappa_1 + \cdots$$ (II-3-3)

donde $\alpha_0$ y $\alpha_1$, son coeficientes, mientras que $\kappa_0$ y $\kappa_1$ se definen como

$$\begin{eqnarray*} \kappa_0 & = & h F(y_0,t_0) = h F \\ \kappa_1 & = & h F(y_0+\lambda\kappa_0,t_0+\mu h) \end{eqnarray*}$$

se tienen dos valores de $F$, uno en cada punto; puede identificarse a $\kappa_0$ como el término $h$, $y'$ que se introdujo en (II-1-2). Al expresar $y$ en términos de esas dos cantidades se evita el problema de trabajar con las derivadas; el procedimiento es tan legítimo como lo es el obtener la función, la derivada y la segunda derivada usando tres puntos de la curva. Que dicho procedimiento sea o no correcto depende de qué tanto coinciden las series (II-3-2) y (II-3-3). Entonces, en lugar de tener que invertir el desarrollo (II-3-1) se supone que $y$ puede expresarse en la forma que hemos indicado y se verá que error resulta al hacer esa consideracion. El error es, naturalmente, la diferencia entre las expresiones (II-3-2) y (II-3-3). Al substituir $\kappa_0$ y $\kappa_1$, en la segunda de estas series se tiene

$$y = y_0+\alpha_0Fh+\alpha_1 h(F+\lambda h F F_y +\mu hF_t)+oh^3+\cdots$$ (II-3-4)

haciendo también una substitución en (II-3-2) y usando (II-3-1) resulta

$$y = y_0+hF+\frac{h^2}{2}(FF_y+F_t)+oh^3+\cdots$$ (II-3-5)

Por comparación de (II-3-4) y (II-3-5) se observa que son muy próximas una de otra si los coeficientes de las mismas potencias de $h$ son aproximadamente iguales, o sea, que

$$\begin{eqnarray*} Fh\alpha_0+Fh\alpha_1 & \approx & hF \\ h^2(\lambda \alpha_1FF_y+\alpha_1\mu F_t) & \approx & (FF_y+F_t)h^2 \end{eqnarray*}$$

Si la aproximación no representa a $y$, por lo menos, la diferencia entre aquella y el caso real tiende a cero a tercer orden en $h$. Como deseamos que el método tenga validez general, debemos pedir que se cumplan las relaciones

$$\alpha_0 + \alpha_1 = 1 \quad , \quad \lambda \alpha_1 = \frac{1}{2} \quad , \quad \mu \alpha_1 = \frac{1}{2}$$

Si esto se cumple, podemos confiar en nuestro método. Como tenemos cuatro incógnitas y sólo tres ecuaciones, el sistema esta indeterminado, pero si se toma $\alpha_1=c$ por ejemplo, entonces todas las soluciones quedan en términos del parámetro $c$ y se tendrá por consiguiente

$$\begin{eqnarray*} \alpha_1 = c \qquad & & \qquad \lambda = \frac{1}{2c} \\ \alpha_0 = 1 - c & & \qquad \mu = \frac{1}{2c}\\ \end{eqnarray*}$$

un valor muy favorecido para $c$ es $\frac 1 2$ porque de esa manera se tiene

$$\begin{eqnarray*} \alpha_1 = \frac{1}{2} & \qquad & \lambda = 1 \\ \alpha_0 = \frac{1}{2} & \qquad & \mu = 1 \\ \end{eqnarray*}$$

La selección de $c$ de ninguna manera aumenta el orden de aproximación (igualar términos en $h^3$) sino que simplifica el cómputo si se toma un valor adecuado. La técnica que acabamos de exponer es conocida como ``Método de Runge-Kutta'' y puede extenderse para incluir aproximaciones de orden mayor en las potencias de $h$. Escribiendo la fórmula de Runge-Kutta para cualquier punto se tiene

$$\begin{eqnarray*} y_{n+1} & = & y_n+\alpha_0\kappa_0+\alpha_1\kappa_1 \end{eqnarray*}$$

donde

$$\begin{eqnarray*} \kappa_0 & = & h F(y_n,t_n) \\ \kappa_1 & = & h F(y_n+\lambda\kappa_0,t_n+\mu h) \end{eqnarray*}$$

Al usar el valor convenido anteriormente para $c$, resulta

$$\begin{array}{rcl} y_{n+1} & = & y_n + \frac{1}{2}\kappa_... ..._n) \\ \kappa_1 & = & h F(y_n+\kappa_0,t_n+h) \end{array}$$ (II-3-6)

Gráficamente, el proceso puede extenderse de la manera siguiente: Primero hay que formar $\kappa_0$ que es el incremento vertical de la figura; con $\kappa_1$ se realiza el segundo paso de Euler en el punto $(y_n+\kappa,t_n+h)$. La fórmula (II-3-6) dice que hay que sumar los dos incrementos y el error obtenido es del orden de $\varepsilon^2$ cuando hay un error de $2\varepsilon$ al ajustar nuestro elemento de arco con la parábola.

Figure 2.6:

Cuando se hacen los cálculos a tercer orden el error irá como $\varepsilon^3$ y el punto resultante se toma como el peso de tres incrementos. Existen varios métodos de Runge-Kutta y su elaboración es más o menos la misma. Hay otros métodos y todos ellos suponen que son conocidas las derivadas. La ventaja de los métodos de Runge-Kutta está en que no hay que calcular derivadas, además, es posible modificar el intervalo (se aumenta o se reduce el incremento ) de acuerdo con la variación de la tangente. El inconveniente que representan estos métodos es el de que hay que realizar un gran número de cálculos para obtener $F$ y cada paso está basado en los anteriores. En general, supóngase que $y$ es un punto sobre la curva y que puede escribirse en términos de uno de sus valores anteriores más multiplos de $F$ en diferentes puntos. Se expresa $y$ en serie de potencias y se comparan los dos desarrollos. En otras palabras, suponemos que:

$$\begin{eqnarray*} y &=& y_0+\alpha F(y_1,t_1)+\beta F(y_2,t_2)+\gamma F(y_3,t_3)+\cdots \end{eqnarray*}$$

Aquí, la aproximación (el órden) depende del número de puntos que se consideren, uno para la primera aproximación, dos para la segunda, etc. Además, se supone, como hemos mencionado, que

$$y = y_0+hy'_0\frac{h^2}{2!}y^{\prime\prime}_0+ \frac{h^3}{3!}y^{\prime\prime\prime}_0+\cdots$$ (II-3-7)

El plan es el mismo que en caso anterior; hay que expresar $y(t)$ en términos de $(y_0,t_0)$ y elaborar $F$ en varios puntos para obtener una serie; los valores de $F$ se obtienen con ayuda de la ecuación diferencial. Hecho todo lo anterior, la tarea consiste en seleccionar los coeficientes de tal manera que haya una correspondencia, válida hasta la potencia deseada.

2.4 Presentamos ahora varios resultados sin entrar en más detalles, ya que el procedimiento es el mismo que se discutió en páginas anteriores. En (II-3-7) se escribe la fórmula hasta tener órden, lo que significa que el error es a cuarto órden en $h$. En ese caso se tiene

$$\begin{array}{rcl} y_{n+1} &=& y_n + \frac{1}{6}(\kappa_0... ...=& h F(y_n+\frac{3}{2}\kappa_0,t_n+\frac{3}{2}h) \end{array}$$ (II-4-1)

Esta fórmula se atribuye a Kutta. Hay otra que se atribuye a Heun y es la siguiente:

$$\begin{array}{rcl} y_{n+1} &=& y_n + \frac{1}{4}(\kappa_0... ...=& h F(y_n+\frac{3}{2}\kappa_1,t_n+\frac{3}{2}h) \end{array}$$ (II-4-2)

Nótese que en el incremento de la función no aparece $\kappa_1$. Este sólo se usa como un paso intermedio. Con un método de Runge-Kutta a cuarto orden se obtiene

$$y_{n+1} = y_n + \frac{1}{6}(\kappa_0+2\kappa_1+2\kappa_2+\kappa_3) + o(h^4)$$ (II-4-3)

$$\begin{array}{rcl} \kappa_0 & = & h F(y_n,t_n) \\ & & ... ...& & \\ \kappa_2 & = & h F(y_n+\kappa_2,t_n+h) \end{array}$$ (II-4-4)

Este es el método que usaremos para integrar nuestras ecuaciones diferenciales. Cuando $F$ no depende del tiempo se reduce al método de Simpson. Puede hacerse un diagrama semejante a los que se han hecho anteriormente pero resulta muy complicado, de modo que con los que ya se tienen dan una idea acerca del tipo de gráfica que se obtendría aquí. El único problema que nos queda es el de introducir las extrapolaciones convenientes. En nuestro caso resultan ecuaciones de movimiento vectoriales, lo que significa que existe un vector que depende del tiempo según la ecuación.

$$\begin{eqnarray*} {\bf y}^{\prime} & = & {\bf F}({\bf y},t) \end{eqnarray*}$$

que es equivalente a

$$\begin{eqnarray*} \frac{dy_1}{dt} & = & F_1(y_1,y_2,y_3,\ldots,y_6,t) \\ \f... ...vdots \\ \frac{dy_6}{dt} & = & F_6(y_1,y_2,y_3,\ldots,y_6,t) \end{eqnarray*}$$

Las componentes de ${\bf y}$ son tres momentos, tres coordenadas y el tiempo. Este es un sistema de ecuaciones con seis variables y el tiempo. Además, aquí las derivadas son parciales porque se tienen seis variables y el tiempo en cada ecuación; el hecho es que tenemos el mismo problema repetido seis veces; ahora bien, si se evita la derivación en todo su detalle, pueden considerarse las cosas de tal modo que las ecuaciones salgan exactamente iguales, tomando en cuenta que las ${\bf\kappa}$ son vectores, las ${\bf y}$ son vectores. Se tiene un argumento vectorial ${\bf y}$; ${\bf F}$ indica que existen seis funciones que dependen de seis variables. Si se considera que los seis argumentos son vectores y que la función es un vector, es entonces correcto pensar que las fórmulas para una variable expuestas en hojas anteriores tiene también validez aquí. Teniendo en cuenta la naturaleza vectorial de nuestra función, podemos ahora traducirla convenientemente. Otra cosa que puede verse, aunque no influye, es que en este problema, la funcion hamiltoniana no depende del tiempo. Nuestras ecuaciones, por consiguiente, sólo son dependientes de las $y_i$, por lo que no se escribirá explicitamente al tiempo como argumento. La discusión de los métodos que hemos presentado es solamente descriptiva. Para una presentación rigurosa puede consultarse cualquier libro de Análisis Numérico, por ejemplo, el de Kunz [10] y Hildebrand [11].