Métodos numéricos para ODEs

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 $latex n$ ecuaciones acopladas de la forma:

$latex frac{d}{dt}y_i = f_i(t,y_1,ldots,y_n)$ con $latex 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 $latex y_i^n$ en $latex t^n$ a partir de la condición inicial $latex y_i^0$ en $latex t^0$.

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

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

donde $latex vec{y}^{n+1} = vec{y}(t^{n+1})$ y $latex 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:

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

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

La posibilidad contraria es:

$latex 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:

$latex 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:

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

$latex 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 $latex frac{Delta t}{2}$, calculamos $latex 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 $latex vec{f}^{mid} = vec{f}(t^{frac{n+1}{2}},vec{y}^{mid})$
  2. Tomamos el paso entero con esta derivada: $latex vec{y}^{n+1} = vec{y}^n + Delta t vec{f}^{mid}$

En el caso de posiciones, velocidades y aceleraciones tenemos:

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

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

Runge-Kutta de cuarto orden:

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

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

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

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

$latex 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:

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

$latex 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.

Únete a la conversación

5 comentarios

  1. The Euler method you are commenting here is EXPLICIT, not implicit (one gets the $latex y^{n+1}$ data ONLY from values at $latex t^{n}$.

  2. The fact that the Störmer-Verlet method provides better conservation properties than a standard RK method, is only proven in NEWTONIAN dynamics, but this is not necessarily true in Relativistic Dynamics (Rosswog 2009 refers to Hairer et al 2006). I do not know whether something similar has been found in Relativistic Dynamics (it would be a good exercice for you to try to find out if one can extend the resutls of Hairer et al. to the relativistic case). In any instance, I recall you that A. Bauswein is using an RK4 method to integrate in time in his relativistic SPH…

    1. Ok. Como te comenté en el último correo, había visto que Bauswein utilizaba RK4 y de ahí la duda con lo que comentaba Rosswog, pero no había caido en la cuenta de que, como bien dices, en un caso estamos en dinámica relativista y en el otro caso lo comenta para dinámica newtoniana… Miro lo de Hairer.

Dejar un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *


¡IMPORTANTE! Responde a la pregunta: ¿Cuál es el valor de 3 5 ?