Método de Verlet
Para o método de Euler implícito havíamos utilizado a derivada a esquerda:
![{\displaystyle f'\left(t\right)\approx {\frac {f\left(t\right)-f\left(t-\Delta t\right)}{\Delta t}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fab187067de51e6d251ad4350dc8ee85a3020c09)
Então se
a segunda derivada é
, pela definição, da derivada a direita:
![{\displaystyle y'\left(t\right)=\lim _{\Delta t\rightarrow 0}{\frac {y\left(t+\Delta t\right)-y\left(t\right)}{\Delta t}}=\lim _{\Delta t\rightarrow 0}{\frac {f'\left(t+\Delta t\right)-f'\left(t\right)}{\Delta t}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3806ec04ff74e6f71d25bba9243b0855bdddb545)
Logo utilizando as aproximações:
![{\displaystyle {\begin{aligned}f''\left(t\right)&\approx {\frac {1}{\Delta t}}\left(f'\left(t+\Delta t\right)-f'\left(t\right)\right)\\&\approx {\frac {1}{\Delta t}}\left({\frac {f\left(t+\Delta t\right)-f\left(t+\Delta t-\Delta t\right)}{\Delta t}}-{\frac {f\left(t\right)-f\left(t-\Delta t\right)}{\Delta t}}\right)\\&\approx {\frac {f\left(t+\Delta t\right)+f\left(t-\Delta t\right)-2f\left(t\right)}{\left(\Delta t\right)^{2}}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/648b0a4b184478612392d2e478b5998a5d198c98)
Isolando então
:
![{\displaystyle {\begin{aligned}f\left(t+\Delta t\right)&=\left(\Delta t\right)^{2}f''\left(t\right)-f\left(t-\Delta t\right)+2f\left(t\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d319a5ad580959de12efa77e2f81b0e98f2cd18c)
Temos o método de Verlet. Podemos notar que precisamos conhecer
em dois tempos anteriores. Podemos utilizar outro algoritmo para o primeiro passo. Se
é posição, então
logo podemos reescrever:
![{\displaystyle {\begin{aligned}x\left(t+\Delta t\right)&=a\left(x,t\right)\left(\Delta t\right)^{2}-x\left(t-\Delta t\right)+2x\left(t\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2989662230fe4eb0121cf80a0e469c43dc5d835d)
Para calcular a energia, podemos obter a velocidade então utilizando a derivada centrada:
![{\displaystyle v\left(t\right)={\dot {x}}\left(t\right)={\frac {x\left(t+\Delta t\right)-x\left(t-\Delta t\right)}{2\Delta t}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c883f55aef6ccd7dbd4a5c316e5ab1d64ad86cb8)
Alternativamente podemos obter o mesmo resultado em termos da expansão de Taylor:
![{\displaystyle {\begin{aligned}x\left(t+\Delta t\right)&=x\left(t\right)+v\left(t\right)\Delta t+{\frac {1}{2}}a\left(t\right)\Delta t^{2}+{\frac {1}{6}}a'\left(t\right)\Delta t^{3}+\dots \\x\left(t-\Delta t\right)&=x\left(t\right)-v\left(t\right)\Delta t+{\frac {1}{2}}a\left(t\right)\Delta t^{2}-{\frac {1}{6}}a'\left(t\right)\Delta t^{3}+\dots \end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/006173a43823bdf9d4e0901a481b4d9337d95a1e)
Somando os dois termos, ficamos então com:
![{\displaystyle x\left(t+\Delta t\right)+x\left(t-\Delta t\right)=2x\left(t\right)+a\left(t\right)\Delta t^{2}+{\mathcal {O}}\left(\Delta t^{4}\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/13f0a2c9caa25e5d01c3d7b74ced104994a74721)
Obtemos então não só o algoritmo de Verlet, além de que sabemos que é uma expansão até a terceir ordem. Então o erro envolvido na truncação é
, e este é o erro local, associao a um único passo.
Além disso, se fizermos a diferença, obtemos o algoritmo da velocidade:
![{\displaystyle x\left(t+\Delta t\right)+x\left(t-\Delta t\right)=2v\left(t\right)+{\frac {1}{3}}a'\left(t\right)\Delta t^{3}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b1f453fc3d910a29e0f6d31ca6a169f49f3b6816)
Então:
![{\displaystyle v\left(t\right)={\frac {x\left(t+\Delta t\right)+x\left(t-\Delta t\right)}{\Delta t}}-{\frac {1}{3}}a'\left(t\right)\Delta t^{2}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/be84902c01105f220dacc939657e9c50e6660a31)
Logo temos um erro
na velocidade. Além do erro de truncação associado ao método de dierenças finita e que decai com o decaimento de
, também podemos lembrar que um erro de arredondamento, pois o computador usa uma quantidade finita de memória para representar os números. Isto, é, existe um número
em que para qualquer número
então
.
é o maior número que pode ser somado a
sem alterar o resultado.
Exemplo
Aplicando o algoritmo para o sistema massa-mola visto no método de Euler-Cromer:
![{\displaystyle {\frac {d^{2}x}{dt^{2}}}=-{\frac {k}{m}}x=-\omega ^{2}x}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e1d1fabd87d796e63c17f0a8d1c31ac438b82ef4)
Podemos ressaltar ainda que
e
.
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
#Valores iniciais
x=[1]; v=[0]; t=[0] ; E=[k*(x[0]**2)/2+m*(v[0]**2)/2]
#Parâmetros
dt = 0.1 ; tau = 2*np.pi; tf=4*tau ; Np= int(tf/dt)
#Método de Euler-Cormer para obter o primeiro passo:
x.append(x[0]+dt*v[0])
t.append(dt)
#Método de Verlet:
for it in range(1,Np):
x.append(-w2*x[it]*dt*dt-x[it-1]+2*x[it]) #Método de Verlet
v.append((x[it+1]-x[it-1])/(2*dt))
E.append(k*x[it]**2/2+m*v[it]**2/2)
t.append(dt+it*dt)
#plt.plot(t,x)
#plt.plot(t[:len(t)-1],v) #Velocidade tem um elemento a menos
#plt.plot(t[:len(t)-1],E)
plt.plot(x[:len(x)-1],v)
Método Velocidade de Verlet
O método de verlet é similar ao leapfrog, mas é síncrono e não exige inicialização. Pela derivada à direita com uma variação de
, obtemos uma equação para velocidade:
![{\displaystyle a\left(t\right)\approx {\frac {v\left(t+\Delta t/2\right)-v\left(t\right)}{\Delta t/2}}\longrightarrow v\left(t+\Delta t/2\right)=v\left(t\right)+a\left(t\right){\frac {\Delta t}{2}}\qquad \left(1\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/f48d0d31c9d4b070adab6514c6ae366860f5c821)
Fazendo o mesmo processo do leapfrog para a posição:
![{\displaystyle v\left(t+\Delta t/2\right)\approx {\frac {x\left(t+\Delta t\right)-x\left(t\right)}{\Delta t}}\longrightarrow x\left(t+\Delta t\right)=x\left(t\right)+v\left(t+\Delta t/2\right)\Delta t}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ee29a47257c165d808972ff54993e207a9c63498)
Ainda podemos reescrever, substituindo:
![{\displaystyle x\left(t+\Delta t\right)=x\left(t\right)+v\left(t\right)\Delta t+a\left(t\right){\frac {\Delta t^{2}}{2}}\qquad \left(2\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/24fd010a4f5b6a79b037928958d8ef831749c618)
Isto é, partindo de um tempo
onde conhecemos posição e velocidade, calculamos a velocidade em
e usamos essa velociade para calcular a poição em
. Agora para a velocidade no instante
, começamo com uma derivada centrada:
![{\displaystyle v\left(t\right)={\frac {x\left(t+\Delta t\right)-x\left(t-\Delta t\right)}{2\Delta t}}\longrightarrow v\left(t+\Delta t\right)={\frac {x\left(t+2\Delta t\right)-x\left(t\right)}{2\Delta t}}\qquad \left(3\right)}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b46bce7c5f1cba3f40d006eae30a2f3839af2b0c)
Então pegamos o método de Verlet e avançamos
:
![{\displaystyle {\begin{aligned}x\left(t+\Delta t\right)&=a\left(t\right)\Delta t^{2}-x\left(t-\Delta t\right)+2x\left(t\right)\\x\left(t+2\Delta t\right)&=a\left(t+\Delta t\right)\Delta t^{2}-x\left(t\right)+2x\left(t+\Delta t\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8fd3a5c9d54cd2aa268564efc55f67c0ea262f8c)
E substituímos na velocidade (equação 3):
![{\displaystyle {\begin{aligned}v\left(t+\Delta t\right)&={\frac {a\left(t+\Delta t\right)\Delta t^{2}-x\left(t\right)+2x\left(t+\Delta t\right)-x\left(t\right)}{2\Delta t}}\\=&a\left(t+\Delta t\right){\frac {\Delta t}{2}}-+{\frac {x\left(t+\Delta t\right)-x\left(t\right)}{\Delta t}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/918d25c088a830e7eed57b5a47c93ed77656afc3)
E usando então a equação para posição que encontramos (equação 2) para substituir
:
![{\displaystyle {\begin{aligned}v\left(t+\Delta t\right)&=a\left(t+\Delta t\right){\frac {\Delta t}{2}}+{\frac {x\left(t+\Delta t\right)-x\left(t\right)}{\Delta t}}\\=&v\left(t\right)+\left(a\left(t+\Delta t\right)+a\left(t\right)\right){\frac {\Delta t}{2}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a7e0dcf4ef7d679591c833e0c4103f0efe37f0fb)
Então agora utilizamos a posição em
para encontrarmos a velocidade no mesmo instante. Temos então:
![{\displaystyle {\begin{aligned}x\left(t+\Delta t\right)&=x\left(t\right)+v\left(t\right)\Delta t+a\left(t\right){\frac {\Delta t^{2}}{2}}\\v\left(t+\Delta t\right)&=v\left(t\right)+\left(a\left(t+\Delta t\right)+a\left(t\right)\right){\frac {\Delta t}{2}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/d88dd7d88ce5d9c8f19fd6f3f636b638bf1ec245)
Ou ainda simplesmente escrevendo explicitamente a equação 1:
![{\displaystyle {\begin{aligned}v\left(t+{\frac {\Delta t}{2}}\right)&=v\left(t\right)+a\left(t\right){\frac {\Delta t}{2}}\\x\left(t+\Delta t\right)&=x\left(t\right)+v\left(t+{\frac {\Delta t}{2}}\right)\Delta t\\v\left(t+\Delta t\right)&=v\left(t+{\frac {\Delta t}{2}}\right)+a\left(t+\Delta t\right){\frac {\Delta t}{2}}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/798327c82f375181cc5391954a8e7cf050c058f0)
Exemplo
Resolvendo o mesmo exemplo anterior, temos:
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
#Valores iniciais
x=[1]; v=[0]; t=[0] ; E=[k*(x[0]**2)/2+m*(v[0]**2)/2]
#Parâmetros
dt = 0.1 ; tau = 2*np.pi; tf=4*tau ; Np= int(tf/dt)
#Método de Velocidade Verlet:
for it in range(Np):
x.append(x[it]+v[it]*dt-w2*x[it]*dt**2/2)
v.append(v[it]+(-w2*x[it+1]-w2*x[it])*dt/2)
E.append(k*x[it]**2/2+m*v[it]**2/2)
t.append(dt+it*dt)
#plt.plot(t,x)
#plt.plot(t,v)
#plt.plot(t,E)
plt.plot(x,v)
Principais materiais utilizados
- The Second Derivative (Matthew Boelkins, David Austin & Steven Schlicker; LibreTexts)
- Verlet Method (Brasnislav K. Nikolic, Universidae de Delaware)