Método de Verlet

O método

Para o método de Euler implícito havíamos utilizado a derivada a esquerda:

$$f^{\prime }\left(t\right)\approx \frac{f\left(t\right)-f(t-\mathrm{\Delta }t)}{\mathrm{\Delta }t}$$

Então se $y\left(t\right)=f^{\prime }\left(t\right)$ a segunda derivada é $y^{\prime }\left(t\right)=f^{″}\left(t\right)$, pela definição, da derivada a direita:

$$y^{\prime }\left(t\right)=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{y(t+\mathrm{\Delta }t)-y\left(t\right)}{\mathrm{\Delta }t}=\underset{\mathrm{\Delta }t\to 0}{\lim }\frac{f^{\prime }(t+\mathrm{\Delta }t)-f^{\prime }\left(t\right)}{\mathrm{\Delta }t}$$

Logo utilizando as aproximações:

$$\begin{array}{rl}{f}^{″}\left(t\right)& \approx \frac{1}{\mathrm{\Delta }t}(f^{\prime }(t+\mathrm{\Delta }t)-f^{\prime }\left(t\right))\\ & \approx \frac{1}{\mathrm{\Delta }t}(\frac{f(t+\mathrm{\Delta }t)-f(t+\mathrm{\Delta }t-\mathrm{\Delta }t)}{\mathrm{\Delta }t}-\frac{f\left(t\right)-f(t-\mathrm{\Delta }t)}{\mathrm{\Delta }t})\\ & \approx \frac{f(t+\mathrm{\Delta }t)+f(t-\mathrm{\Delta }t)-2f\left(t\right)}{{\left(\mathrm{\Delta }t\right)}^2}\end{array}$$

Isolando então $f(t+\mathrm{\Delta }t)$:

$$\begin{array}{rl}f(t+\mathrm{\Delta }t)& ={\left(\mathrm{\Delta }t\right)}^2{f}^{″}\left(t\right)-f(t-\mathrm{\Delta }t)+2f\left(t\right)\end{array}$$

Temos o método de Verlet. Podemos notar que precisamos conhecer $f\left(t\right)$ em dois tempos anteriores. Podemos utilizar outro algoritmo para o primeiro passo. Se $f\left(t\right)=x\left(t\right)$ é posição, então $f^{″}\left(t\right)=\ddot{x}\left(t\right)=a(x,t)$ logo podemos reescrever:

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

Para calcular a energia, podemos obter a velocidade então utilizando a derivada centrada: $$v\left(t\right)=\dot{x}\left(t\right)=\frac{x(t+\mathrm{\Delta }t)-x(t-\mathrm{\Delta }t)}{2\mathrm{\Delta }t}$$ Alternativamente podemos obter o mesmo resultado em termos da expansão de Taylor:

$$\begin{array}{rl}x(t+\mathrm{\Delta }t)& =x\left(t\right)+v\left(t\right)\mathrm{\Delta }t+\frac{1}{2}a\left(t\right)\mathrm{\Delta }{t}^2+\frac{1}{6}{a}^{\prime }\left(t\right)\mathrm{\Delta }{t}^3+\dots \\ x(t-\mathrm{\Delta }t)& =x\left(t\right)-v\left(t\right)\mathrm{\Delta }t+\frac{1}{2}a\left(t\right)\mathrm{\Delta }{t}^2-\frac{1}{6}{a}^{\prime }\left(t\right)\mathrm{\Delta }{t}^3+\dots \end{array}$$

Somando os dois termos, ficamos então com:

$$x(t+\mathrm{\Delta }t)+x(t-\mathrm{\Delta }t)=2x\left(t\right)+a\left(t\right)\mathrm{\Delta }{t}^2+{𝒪}\left(\mathrm{\Delta }{t}^4\right)$$

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 é ${𝒪}\left(\mathrm{\Delta }{t}^4\right)$, e este é o erro local, associao a um único passo.

Além disso, se fizermos a diferença, obtemos o algoritmo da velocidade:

$$x(t+\mathrm{\Delta }t)+x(t-\mathrm{\Delta }t)=2v\left(t\right)+\frac{1}{3}{a}^{\prime }\left(t\right)\mathrm{\Delta }{t}^3$$

Então:

$$v\left(t\right)=\frac{x(t+\mathrm{\Delta }t)+x(t-\mathrm{\Delta }t)}{\mathrm{\Delta }t}-\frac{1}{3}{a}^{\prime }\left(t\right)\mathrm{\Delta }{t}^2$$

Logo temos um erro ${𝒪}\left(\mathrm{\Delta }{t}^2\right)$ na velocidade. Além do erro de truncação associado ao método de dierenças finita e que decai com o decaimento de $\mathrm{\Delta }t$, 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 $\alpha \le ϵ$ então $1+\alpha =1$. $ϵ$ é o maior número que pode ser somado a $1$ sem alterar o resultado.

Exemplo

Aplicano o algoritmo para: $$\frac{d^{2}x}{d{t}^2}=-\frac{k}{m}x={\omega }^{2}x$$

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)

Principais materiais utilizados

1. The Second Derivative (Matthew Boelkins, David Austin & Steven Schlicker; LibreTexts)