Ya vimos que la discretización de las ecuaciones de conservación en forma Lagrangiana para SPH quedaba:
$latex frac{d}{dt}boldsymbol{v}_a = – sum_b m_b (frac{P_a}{rho_a^2} + frac{P_b}{rho_b^2}) nabla_a W{ab} = boldsymbol{F}_a$
$latex frac{d}{dt} e_a = frac{1}{2} sum_b m_b (frac{P_a}{rho_a^2} + frac{P_b}{rho_b^2}) boldsymbol{v}_{ab}cdot nabla_a W{ab} = E_a$
con $latex rho_a = sum_b m_b W_{ab}$ y $latex P_a$ de las EoS.
Para pensar en terminos de ODE solvers, nos fijamos solo en:
$latex frac{d}{dt} boldsymbol{v}_a = boldsymbol{F}_a(t,boldsymbol{v}_a)$
$latex frac{d}{dt} e_a = E_a(t,e_a)$
$latex frac{d}{dt} boldsymbol{r}_a = boldsymbol{v}_a(t,boldsymbol{r}_a)$
Para simplificar, suponiendo como metodo explicito un Euler para el predictor y el método del trapecio como implícito para el corrector, tendriamos:
$latex boldsymbol{tilde{v}}_a^{n+1} = boldsymbol{v}_a^{n} + h boldsymbol{F}_a(t^n,boldsymbol{v}_a^{n})$
$latex tilde{e}_a^{n+1} = e_a^{n} + h E_a(t^n,e_a^{n})$
$latex boldsymbol{tilde{r}}_a^{n+1} = boldsymbol{r}_a^{n} + h boldsymbol{v}_a(t^n,boldsymbol{r}_a^{n})$
para la primera y:
$latex boldsymbol{v}_a^{n+1} = boldsymbol{v}_a^n + frac{1}{2}h(boldsymbol{F}_a(t^n,boldsymbol{v}_a^n)+boldsymbol{F}_a(t^{n+1},boldsymbol{tilde{v}}_a^{n+1}))$
$latex e_a^{n+1} = e_a^n + frac{1}{2}h(E_a(t^n,e_a^n)+E_a(t^{n+1},tilde{e}_a^{n+1}))$
$latex boldsymbol{r}_a^{n+1} = boldsymbol{r}_a^n + frac{1}{2}h(boldsymbol{v}_a(t^n,boldsymbol{r}_a^n)+boldsymbol{v}_a(t^{n+1},boldsymbol{tilde{r}}_a^{n+1}))$
para la segunda.
De esta manera, sustituyendo, tenemos:
$latex boldsymbol{tilde{v}}_a^{n+1} = boldsymbol{v}_a^n + h boldsymbol{F}_a(t^n,boldsymbol{v}_a^{n}) =$
$latex = boldsymbol{v}_a^n – h sum_b m_b (frac{P_a^n}{(rho_a^n)^2} + frac{P_b^n}{(rho_b^n)^2}) nabla_a W(boldsymbol{r}_a^n-boldsymbol{r}_b^n,h)$