You are currently browsing the category archive for the ‘ODE’ category.

El caso mas sencillo es cuando tenemos cinco puntos:

(x_0,y_0), (x_1,y_1), (x_2,y_2), (x_3,y_3), (x_4,y_4),

de manera que:

x_1-x_0 = 2(x_2-x_1) = 2(x_3-x_2) = x_4-x_3.

Tal y como escribimos en el anterior post, el polinomio general para cinco puntos quedaria:

L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x) + y_3 l_3(x) + y_4 l_4(x),

donde:

l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4},

l_1(x) = \frac{x-x_0}{x_1-x_0} \frac{x-x_2}{x_1-x_2} \frac{x-x_3}{x_1-x_3} \frac{x-x_4}{x_1-x_4},

l_2(x) = \frac{x-x_0}{x_2-x_0} \frac{x-x_1}{x_2-x_1} \frac{x-x_3}{x_2-x_3} \frac{x-x_4}{x_2-x_4},

l_3(x) = \frac{x-x_0}{x_3-x_0} \frac{x-x_1}{x_3-x_1} \frac{x-x_2}{x_3-x_2} \frac{x-x_4}{x_3-x_4},

l_4(x) = \frac{x-x_0}{x_4-x_0} \frac{x-x_1}{x_4-x_1} \frac{x-x_2}{x_4-x_2} \frac{x-x_3}{x_4-x_3}.

Dado que tratamos de aproximar la segunda derivada, tenemos, en el caso de equidistancia quedaría:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-2}+16u_{i-1}-30u_i+16u_{i+1}-u_{i+2}}{12 h^2}.

Vamos a ver que queda ahora en nuestro caso. Reescribimos los l_i(x) en función de x_0 de manera que, por ejemplo, tenemos:

l_0(x) = \frac{x-x_1}{x_0-x_1} \frac{x-x_2}{x_0-x_2} \frac{x-x_3}{x_0-x_3} \frac{x-x_4}{x_0-x_4} =

=\frac{x-x_0-2h}{x_0-x_0-2h} \frac{x-x_0-3h}{x_0-x_0-3h} \frac{x-x_0-4h}{x_0-x_0-4h} \frac{x-x_0-6h}{x_0-x_0-6h} =

=\frac{(x-x_0-2h)(x-x_0-3h)(x-x_0-4h)(x-x_0-6h)}{144h^4},

que para x=x_2=x_0+3h (2h+h), la segunda derivada queda:

\frac{d^2}{dx^2} l_0(x) |_{x=x_2} = -\frac{1}{72h^2}.

Haciendo lo mismo para todas las derivadas segundas de todos los l_i(x), obtenemos la siguiente formula en diferencias finitas:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-3} + 81 u_{u-1} - 160 u_i + 81 u_{i+1} - u_{i+3}}{72h^2},

donde los índices hacen referencia a una malla inicial equiespaciada.

Para el mismo denominador, antes teniamos:

\frac{d^2}{dx^2}u_i = \frac{-6u_{i-2}+96u_{i-1}-180u_i+96u_{i+1}-6u_{i+2}}{72 h^2}.

Podemos hacer lo mismo para cada uno de los puntos x=x_0, x=x_1=x_0+2h, x=x_3=x_0+4h y x=x_4=x_0+6h:

\frac{d^2}{dx^2}u_{i-3} = \frac{40 u_{i-3} -243 u_{i-1} + 352 u_i -162 u_{i+1} + 13 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i-2} = \frac{7 u_{i-3} - 32 u_i + 27 u_{i+1} - 2 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i+1} =\frac{-2 u_{i-3} + 27 u_{i-1} - 32 u_i -+7 u_{i+3}}{36 h^2}

\frac{d^2}{dx^2}u_{i+3} =\frac{13 u_{i-3} -162 u_{i-1} + 352 u_i -243 u_{i+1} + 40 u_{i+3}}{36 h^2}

 En el caso de x=x_0+4h (3h+h), tenemos:

\frac{d^2}{dx^2}u_i = \frac{-u_{i-4} + 256_{u-1} - 510 u_i + 256 u_{i+1} - u_{i+4}}{140h^2}.

Dados n+1 puntos:

(x_0,y_0), \ldots, (x_n,y_n),

definimos el polinomio interpolador de Lagrange:

L(x) := \Sigma_{i=0}^n y_i l_i(x)

donde:

l_i(x) := \Pi_{0 \leq j \leq n,j \neq i} \frac{x-x_j}{x_i-x_j}.

Con esto, tenemos:

\frac{d^n}{dx^n}L(x) = \Sigma_{i=0}^n y_i \frac{d^n}{dx^n}l_i(x)

Supongamos que tenemos n=2:

(x_0,y_0), (x_1,y_1).

En este caso, tenemos:

L(x) = y_0 l_0(x) + y_1 l_1(x)

con:

l_0(x) = \frac{x-x_1}{x_0-x_1}

l_1(x) = \frac{x-x_0}{x_1-x_0}

Como tenemos dos puntos, únicamente podemos aproximar la primera derivada:

\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x)

con:

\frac{d}{dx} l_0(x) = \frac{1}{x_0-x_1}

\frac{d}{dx} l_1(x) = \frac{1}{x_1-x_0},

de manera que:

\frac{d}{dx}L(x) = \frac{y_1 - y_0}{x_1 - x_0}.

Con esto, tenemos las aproximaciones de primer orden:

u_x(x_0) = u_x(x_1) = \frac{y_1 - y_0}{x_1-x_0},

que, en índices, queda:

\frac{d}{dx}u_i \approx \frac{u_{i+1}-u_i}{h_x},

\frac{d}{dx}u_{i+1} \approx \frac{u_{i+1}-u_i}{h_x},

donde h_x = x_1 - x_0 (y que es lógico, ya que la derivada será una constante).

Supongamos ahora que tenemos tres puntos:

(x_0,y_0),(x_1,y_1),(x_2,y_2).

En este caso tenemos:

L(x) = y_0 l_0(x) + y_1 l_1(x) + y_2 l_2(x),

\frac{d}{dx} L(x) = y_0 \frac{d}{dx} l_0(x) + y_1 \frac{d}{dx} l_1(x) + y_2 \frac{d}{dx} l_2(x),

\frac{d^2}{dx^2} L(x) = y_0 \frac{d^2}{dx^2} l_0(x) + y_1 \frac{d^2}{dx^2} l_1(x) + y_2 \frac{d^2}{dx^2} l_2(x).

Ahora tenemos:

l_0(x) = \frac{(x-x_1)(x-x_2)}{(x_0-x_1)(x_0-x_2)} = \frac{x^2 + (-x_1-x_2)x + x_1x_2}{(x_0-x_1)(x_0-x_2)},

l_1(x) = \frac{(x-x_0)(x-x_2)}{(x_1-x_0)(x_1-x_2)} = \frac{x^2 + (-x_0-x_2)x + x_0x_2}{(x_1-x_0)(x_1-x_2)},

l_2(x) = \frac{(x-x_0)(x-x_1)}{(x_2-x_0)(x_2-x_1)} = \frac{x^2 + (-x_0-x_1)x + x_0x_1}{(x_2-x_0)(x_2-x_1)}.

\frac{d}{dx} l_0(x) = \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)},

\frac{d}{dx} l_1(x) = \frac{2x+(-x_0-x_2)}{(x_1-x_0)(x_1-x_2)},

\frac{d}{dx} l_2(x) = \frac{2x+(-x_0-x_1)}{(x_2-x_0)(x_2-x_1)}.

\frac{d^2}{dx^2} l_0(x) = \frac{2}{(x_0-x_1)(x_0-x_2)},

\frac{d^2}{dx^2} l_1(x) = \frac{2}{(x_1-x_0)(x_1-x_2)},

\frac{d^2}{dx^2} l_2(x) = \frac{2}{(x_2-x_0)(x_2-x_1)}.

De esta manera, por ejemplo, tenemos que:

\frac{d^2}{dx^2}L(x) = y_0 \frac{2}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2}{(x_2-x_0)(x_2-x_1)},

que, en el caso de tener los puntos equiespaciados (h_x:=x_{i+1}-x_{i} con 0 \leq i \leq 1), queda la aproximación de segundo orden:

\frac{d^2}{dx^2}u(x) = \frac{y_0 - 2 y_1 + y_2}{h_x^2},

que además, como era de esperar, es independiente de x:

\frac{d^2}{dx^2}u_{i-1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2},

\frac{d^2}{dx^2}u_i \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2},

\frac{d^2}{dx^2}u_{i+1} \approx \frac{u_{i-1}-2u_i+u_{i+1}}{h_x^2}.

¿Qué pasa con la primera derivada? En este caso, el resultado si que depende de x (por lo que tendremos un resultado diferente en función de si la x vale x_0, x_1 o x_2) y lo que obtenemos es:

\frac{d}{dx} L(x) = y_0 \frac{2x + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)},

por lo que:

\frac{d}{dx} L(x_0) = y_0 \frac{2x_0 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_0 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_0 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}

\frac{d}{dx} L(x_1) = y_0 \frac{2x_1 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_1 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_1 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)}

\frac{d}{dx} L(x_2) = y_0 \frac{2x_2 + (-x_1-x_2)}{(x_0-x_1)(x_0-x_2)} + y_1 \frac{2x_2 + (-x_0-x_2)}{(x_1-x_0)(x_1-x_2)} + y_2 \frac{2x_2 + (-x_1-x_2)}{(x_2-x_0)(x_2-x_1)},

que, con equiespaciado, quedan:

L'(x_0) = \frac{3 y_0 -4 y_1 + y_2}{2 h_x},

L'(x_1) = \frac{y_0 -2 y_1 + y_2}{2 h_x} y

L'(x_2) = \frac{y_0 -4 y_1 + 3 y_2}{2 h_x}.

Laplaciano en cartesianas:

\Delta u = \Sigma_i \frac{\partial^2}{\partial x_i^2}u

1d

\frac{u_{i-1}-2u_i+u_{i+1}}{h^2} = f_i

\frac{1}{h^2}u_{i-1} + \frac{1}{h^2}u_{i+1} +\frac{-2}{h^2}u_i= f_i

2d

\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{h_x^2} + \frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{h_y^2} = f_{i,j}

i fijo:

\frac{1}{h_y^2}u_{i,j-1} + \frac{1}{h_y^2}u_{i,j+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_x^2}u_{i-1,j} + \frac{-1}{h_x^2}u_{i+1,j})

j fijo:

\frac{1}{h_x^2}u_{i-1,j} + \frac{1}{h_x^2}u_{i+1,j} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2})u_{i,j}= g_{i,j}(:=f_{i,j} + \frac{-1}{h_y^2}u_{i,j-1} + \frac{-1}{h_y^2}u_{i,j+1})

3d

\frac{u_{i-1,j,k}-2u_{i,j,k}+u_{i+1,j,k}}{h_x^2} + \frac{u_{i,j-1,k}-2u_{i,j,k}+u_{i,j+1,k}}{h_y^2} + \frac{u_{i,j,k-1}-2u_{i,j,k}+u_{i,j,k+1}}{h_z^2} = f_{i,j,k}

i,j fijos:

\frac{1}{h_z^2}u_{i,j,k-1} + \frac{1}{h_z^2}u_{i,j,k+1} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k})

i,k fijos:

\frac{1}{h_y^2}u_{i,j-1,k} + \frac{1}{h_y^2}u_{i,j+1,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_x^2}u_{i-1,j,k} + \frac{-1}{h_x^2}u_{i+1,j,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})

j,k fijos:

\frac{1}{h_x^2}u_{i-1,j,k} + \frac{1}{h_x^2}u_{i+1,j,k} +(\frac{-2}{h_x^2}+\frac{-2}{h_y^2}+\frac{-2}{h_z^2})u_{i,j,k}= g_{i,j,k}

(g_{i,j,k}:=f_{i,j,k} + \frac{-1}{h_y^2}u_{i,j-1,k} + \frac{-1}{h_y^2}u_{i,j+1,k} + \frac{-1}{h_z^2}u_{i,j,k-1} + \frac{-1}{h_z^2}u_{i,j,k+1})

Dado el campo vectorial invariante por la izquierda

u(a) = (L_a)_* u del grupo de Lie G,

llamamos trayectorias del campo a la soluciones de la ODE:

\frac{d}{dt}a = u(a),

y que, fijada una condición inicial a(0) = a, la trayectoria queda unívocamente determinada en forma de serie exponencial:

a(t) = \exp (tu) a con \exp(tu) = I + tu + \frac{t^2}{2}u^2 + \ldots,

donde u^r(a) = u(a) \circ u^{r-1}(a).

Las trayectorias:

  1. son analíticas  y esta definidas para todo valor de t por el carácter anaítico y completo de los campos invariantes.
  2. que pasan por la unidad del grupo de Lie, a(0) = e, son grupos uniparamétricos de G, ya que cumplen la propiedad a(s+t)=a(s)a(t). En este caso:

\exp: \mathfrak{g} \longrightarrow G \, / \, u \mapsto \exp u = a(1).

En el caso de GL(n,\mathbb{K}), hemos visto que los campos invariantes a la izquierda tienen la forma:

U(A) = AU con A \in GL(n,\mathbb{K}) y U \in M(n,\mathbb{K}).

El subgrupo uniparamétrico A(t) es la solución a:

\frac{d}{dt}A = AU con A(0)=I,

y que es:

A(t) = \exp(tU) = \Sigma_{s=0}^{\infty}\frac{t^s}{s!}U^s.

De esta manera, tenemos:

\exp: \mathfrak{g}(n,\mathbb{K}) \longrightarrow GL(n,\mathbb{K}) \, / \, A = e^U,

es decir, la exponencial matricial habitual.

En la siguiente simulación tenemos 1000 partículas. Cada partícula a representa un péndulo de masa m_a y posición (x_a(t),y_a(t),z_a(t)), con z_a(t) constante, sometido a las fuerzas:

  1. F_1 = -m_ag, la fuerza de la gravedad, con g constante.
  2. F_2, la fuerza elática de la suspensión, de sentido opuesta a la posición del péndulo y de magnitud k_1d, con k_1 > 0 y d la variación de la longitud del péndulo respecto de su longitud en reposo l.
  3. F_3, la fuerza de fricción, con sentido opuesta a la velocidad del péndulo y magnitud k_2v, con v = \sqrt{(x_a')^2+(y_a')^2} y k_2 > 0.

La ley de Newton F = ma queda, para cada uno de los péndulos, es decir, para cada partícula, como el sistema de dos ecuaciones diferenciales de segundo orden:

m_a \begin{bmatrix} x_a''(t) \\ y_a''(t) \end{bmatrix} = F_1 + F_2 + F_3

que es:

m_ax_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})x_a - k_2\frac{x'_a}{\sqrt{x_a'^2+y_a'^2}}

m_ay_a'' = -k_1 (1-\frac{l}{\sqrt{x_a^2+y_a^2}})y_a - k_2\frac{y'_a}{\sqrt{x_a'^2+y_a'^2}}

Podemos expresar este sistema como un sistemad de ecuaciones de primer orden introduciendo incongnitas adicionales: u_a = x_a' y v_a = y_a', con lo que nos queda:

\left \{  \begin{array}{rcl}  x_a' &=& u_a \\  y_a' &=& v_a \\  m_au_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)x_a - k_2\frac{u_a}{\sqrt{u_a^2+v_a^2}} \\  m_av_a' & = & k_1 (\frac{l}{\sqrt{x_a^2+y_a^2}}-1)y_a - k_2\frac{v_a}{\sqrt{u_a^2+v_a^2}} - m_ag \\  \end{array}  \right .

En la siguiente simulación tomamos m_a aleatoria para cada partícula a, por lo que tendremos que resolver el sistema de EDOs numéricamente para cada una de ellas. Además, tendremos l=1, k_1 = 10000 y k_2=1 en todos los péndulos. Cada uno de éstos se distribuye a lo largo de una espilar cilíndrica y las velocidades iniciales son, para todas ellas, u_{a_0}=0, v_{a_0} = -0.7 y v_{z_0} = 0. El intervalo de tiempo es desde t_i = 0 hasta t_f=10 que se divide en 400 subintervalos.

El resultado es el siguiente:

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.

agosto 2017
L M X J V S D
« Jul    
 123456
78910111213
14151617181920
21222324252627
28293031