En [Rosswog 2009] comenta algunos métodos para ODEs, ya que una de las características relevantes de SPH es que transforma las PDEs de la hidrodinámica en un sistema de ODEs y hacer una simulación en SPH no es mas que resolver dicho sistema.

Las ODEs en astrofísica suelen ser de segundo orden pero podemos transformarlas en un sistema de ODEs de primer orden introduciendo como nuevas variables a las velocidades. De esta manera obtenemos n ecuaciones acopladas de la forma:

\frac{d}{dt}y_i = f_i(t,y_1,\ldots,y_n) con i = 1,\ldots,n

Tenemos, en este caso, dos tipos de problemas: los problemas de contorno y los problemas de valor inicial. Este último caso, que es el que nos interesa, se trata de determinar el valor de y_i^n en t^n a partir de la condición inicial y_i^0 en t^0.

El método de Euler es el método de integración pedagógico. La solución avanza de t^n a t^{n+1} \equiv t^n + \Delta t según:

\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}(t^n,\vec{y}^n)

donde \vec{y}^{n+1} = \vec{y}(t^{n+1}) y \Delta t es el paso de tiempo escogido. Es un método explícito de primer orden y es asimétrico en tiempo, pues la derivada es correcta al principio del paso de tiempo pero no al final. De ahí el nombre de Euler hacia adelante. En el caso de posiciones, velocidades y aceleraciones, tenemos:

\vec{x}^{n+1} = \vec{x}^n + \Delta t \vec{v}^n

\vec{v}^{n+1} = \vec{v}^n + \Delta t \vec{a}^n

La posibilidad contraria es:

\vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}(t^{n+1},\vec{y}^{n+1})

llamado el método de Euler hacia atrás, que es un método implícito. Todos los siguientes métodos que comentaremos son explícitos.

El siguiente método es el Euler modificado que utiliza el valor de la derivada al principio y al final del paso de tiempo:

\vec{y}^{n+1} = \vec{y}^n + \Delta t \frac{\vec{f}^n + \vec{f}^{n+1,p}}{2}

que es un método de segundo orden. En el caso de posiciones, velocidades y aceleraciones tenemos:

\vec{x}^{n+1} = \vec{x}^n + \Delta t \frac{\vec{v}^n + \vec{v}^{n+1,p}}{2}

\vec{v}^{n+1} = \vec{v}^n + \Delta t \frac{\vec{a}^n + \vec{a}^{n+1,p}}{2}

En el método de Runge-Kutta de segundo orden tomamos la derivada en el punto medio de los intervalos de tiempo tomados:

  1. Tomamos como paso \frac{\Delta t}{2}, calculamos \vec{y}^{mid} = \vec{y}^n + \frac{1}{2} \Delta t \vec{f}(t^n,\vec{y}^n) y evaluamos la derivada en ese punto \vec{f}^{mid} = \vec{f}(t^{\frac{n+1}{2}},\vec{y}^{mid})
  2. Tomamos el paso entero con esta derivada: \vec{y}^{n+1} = \vec{y}^n + \Delta t \vec{f}^{mid}

En el caso de posiciones, velocidades y aceleraciones tenemos:

\vec{x}^{n+1} = \vec{x}^n + \Delta t \vec{v}^{mid}

\vec{v}^{n+1} = \vec{v}^n + \Delta t \vec{a}^{mid}

Runge-Kutta de cuarto orden:

\vec{k_1} = \Delta t \vec{f}(t^n,\vec{y}^n)

\vec{k_2} = \Delta t \vec{f}(t^n + \frac{\Delta t}{2},\vec{y}^n + \frac{\vec{k}_1}{2})

\vec{k_3} = \Delta t \vec{f}(t^n + \frac{\Delta t}{2},\vec{y}^n + \frac{\vec{k}_2}{2})

\vec{k_4} = \Delta t \vec{f}(t^n + \Delta t,\vec{y}^n + \vec{k}_3)

\vec{y}^{n+1} = \vec{y}^n + \frac{1}{6}(\vec{k}_1+ 2 \vec{k}_2 + 2 \vec{k}_3 + \vec{k}_4)

Velocity Stormer-Verlet:

En el caso de posiciones, velocidades y aceleraciones tenemos:

\vec{x}^{n+1} = \vec{x}^n + \vec{v}^n \Delta t + \frac{1}{2} \vec{a}^n \Delta t^2

\vec{v}^{n+1} = \vec{v}^n + \Delta t \frac{\vec{a}^n + \vec{a}^{n+1}}{2}

Dado que la conservación exacta es en nuestro puede ser mas importante  que el orden alto, aconseja el último método.

Anuncios