Método de Leapfrog

Partindo da ideia que uma inclinação entre dois em uma curva é uma aproximação muito melhor da derivada no ponto médio em alguma das extremidades., e tendo as seguintes equações do movimento:

$$\begin{array}{rl}\frac{dx}{dt}& =v\\ \frac{dv}{dt}& =a\end{array}$$

Sendo a derivada numérica:

$$\begin{array}{rl}v\left(t\right)& \approx \frac{x(+\mathrm{\Delta }t)-x(t-\mathrm{\Delta }t)}{2\mathrm{\Delta }t}\\ & \approx \frac{x(t+\mathrm{\Delta }t/2)-x(t-\mathrm{\Delta }t/2)}{\mathrm{\Delta }t}\end{array}$$

Então para a equação da velocidade temos que:

$$v(t+\mathrm{\Delta }t/2)\approx \frac{x(t+\mathrm{\Delta }t)-x\left(t\right)}{\mathrm{\Delta }t}$$

Ou ainda:

$$x(t+\mathrm{\Delta }t)=x\left(t\right)+v(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t$$

E aplicando a mesma ideia para a aceleração:

$$a\left(t\right)\approx \frac{v(t+\mathrm{\Delta }t/2)-v(t-\mathrm{\Delta }t/2)}{\mathrm{\Delta }t}$$

Logo:

$$v(t+\mathrm{\Delta }t/2)=a\left(t\right)\mathrm{\Delta }t+v(t-\mathrm{\Delta }t/2)$$

Temos então:

$$\begin{array}{rl}v(t+\mathrm{\Delta }t/2)& =v(t-\mathrm{\Delta }t/2)+a\left(t\right)\mathrm{\Delta }t\\ x(t+\mathrm{\Delta }t)& =x\left(t\right)+v(t+\mathrm{\Delta }t/2)\mathrm{\Delta }t\end{array}$$

Exemplo

Aplicando o algoritmo para o sistema massa-mola visto no método de Euler-Cromer: $$\frac{d^{2}x}{d{t}^2}=-\frac{k}{m}x=-{\omega }^{2}x$$

Podemos ressaltar ainda que $a=-{\omega }^{2}x$ e $\frac{dv}{dt}=\frac{d^{2}x}{d{t}^2}$.

import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
import numpy as np                         #Biblitoeca de cálculos científicos

#Constantes
m=1  ; k= 1.; w2= k/m
#Parâmetros
dt  = 0.1 ; tau = 2*np.pi; tf=4*tau ; Np= int(tf/dt)
#Valores iniciais
x=[1]; v0=0; tx=[0]; tv=[dt/2]

#Usamo Euler-Cromer para calcular a velocidade em t=dt/2
x_temp=x[0]+(dt/2)*v0              #Posição    em t=dt/2
v=[v0-(dt/2)*x_temp*w2]            #Velocidade em t=dt/2

#Método Leapfrog
for it  in range(Np):
  x.append(x[it]+dt*v[it])  
  tx.append(dt+it*dt)
  v.append(v[it]-w2*x[it+1]*dt)
  tv.append(dt/2+(1+it)*dt)

#plt.plot(tx,x)
#plt.plot(tv,v)

#Para calcular a energia vamos pegar a média da velocidade
#para termos velocidade e posição no mesmo instante

vm=[v0]
E=[k*(x[0]**2)/2+m*(v[0]**2)/2] 
for it in range(len(x)-1):
  vm.append((v[it+1]+v[it])/2)
  E.append(k*x[it+1]**2/2+m*v[it+1]**2/2)

#plt.plot(tx,E)
plt.plot(x,v)