Método de Verlet

De Física Computacional
Ir para navegação Ir para pesquisar

Método de Verlet

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

f′(t)≈f(t)−f(t−Δt)Δt

Então se y(t)=f′(t) a segunda derivada é y′(t)=f″(t), pela definição, da derivada a direita:

y′(t)=limΔt→0y(t+Δt)−y(t)Δt=limΔt→0f′(t+Δt)−f′(t)Δt

Logo utilizando as aproximações:

f″(t)≈1Δt(f′(t+Δt)−f′(t))≈1Δt(f(t+Δt)−f(t+Δt−Δt)Δt−f(t)−f(t−Δt)Δt)≈f(t+Δt)+f(t−Δt)−2f(t)(Δt)2

Isolando então f(t+Δt):

f(t+Δt)=(Δt)2f″(t)−f(t−Δt)+2f(t)

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

x(t+Δt)=a(x,t)(Δt)2−x(t−Δt)+2x(t)

Para calcular a energia, podemos obter a velocidade então utilizando a derivada centrada: v(t)=x˙(t)=x(t+Δt)−x(t−Δt)2Δt Alternativamente podemos obter o mesmo resultado em termos da expansão de Taylor:

x(t+Δt)=x(t)+v(t)Δt+12a(t)Δt2+16a′(t)Δt3+…x(t−Δt)=x(t)−v(t)Δt+12a(t)Δt2−16a′(t)Δt3+…

Somando os dois termos, ficamos então com:

x(t+Δt)+x(t−Δt)=2x(t)+a(t)Δt2+𝒪(Δt4)

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 é 𝒪(Δt4), 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+Δt)+x(t−Δt)=2v(t)+13a′(t)Δt3

Então:

v(t)=x(t+Δt)+x(t−Δt)Δt−13a′(t)Δt2

Logo temos um erro 𝒪(Δt2) na velocidade. Além do erro de truncação associado ao método de dierenças finita e que decai com o decaimento de Δ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 α≤ϵ então 1+α=1. ϵ é o maior número que pode ser somado a 1 sem alterar o resultado.

Exemplo

Aplicando o algoritmo para o sistema massa-mola visto no método de Euler-Cromer: d2xdt2=−kmx=−ω2x

Podemos ressaltar ainda que a=−ω2x e dvdt=d2xdt2.

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 Δt/2, obtemos uma equação para velocidade:

a(t)≈v(t+Δt/2)−v(t)Δt/2⟶v(t+Δt/2)=v(t)+a(t)Δt2(1)

Fazendo o mesmo processo do leapfrog para a posição:

v(t+Δt/2)≈x(t+Δt)−x(t)Δt⟶x(t+Δt)=x(t)+v(t+Δt/2)Δt

Ainda podemos reescrever, substituindo:

x(t+Δt)=x(t)+v(t)Δt+a(t)Δt22(2)

Isto é, partindo de um tempo t onde conhecemos posição e velocidade, calculamos a velocidade em t+Δt2 e usamos essa velociade para calcular a poição em t+Δt. Agora para a velocidade no instante t+Δt , começamo com uma derivada centrada:

v(t)=x(t+Δt)−x(t−Δt)2Δt⟶v(t+Δt)=x(t+2Δt)−x(t)2Δt(3)

Então pegamos o método de Verlet e avançamos Δt:

x(t+Δt)=a(t)Δt2−x(t−Δt)+2x(t)x(t+2Δt)=a(t+Δt)Δt2−x(t)+2x(t+Δt)

E substituímos na velocidade (equação 3):

v(t+Δt)=a(t+Δt)Δt2−x(t)+2x(t+Δt)−x(t)2Δt=a(t+Δt)Δt2−+x(t+Δt)−x(t)Δt

E usando então a equação para posição que encontramos (equação 2) para substituir x(t+Δt)−x(t) :

v(t+Δt)=a(t+Δt)Δt2+x(t+Δt)−x(t)Δt=v(t)+(a(t+Δt)+a(t))Δt2

Então agora utilizamos a posição em t+Δt para encontrarmos a velocidade no mesmo instante. Temos então:

x(t+Δt)=x(t)+v(t)Δt+a(t)Δt22v(t+Δt)=v(t)+(a(t+Δt)+a(t))Δt2

Ou ainda simplesmente escrevendo explicitamente a equação 1:

v(t+Δt2)=v(t)+a(t)Δt2x(t+Δt)=x(t)+v(t+Δt2)Δtv(t+Δt)=v(t+Δt2)+a(t+Δt)Δt2


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

  1. The Second Derivative (Matthew Boelkins, David Austin & Steven Schlicker; LibreTexts)
  2. Verlet Method (Brasnislav K. Nikolic, Universidae de Delaware)