Método de Euler-Cromer

Lembrando do que vimos no Método de Euler, o sistema de equações para o sistema massa-mola era:

$${\begin{aligned}{\frac {dv}{dt}}&=-\omega ^{2}x\left(t\right)\\{\frac {dx}{dt}}&=v\left(t\right)\end{aligned}}$$

Aplicando o método de Euler então:

$${\begin{aligned}v\left(t+\Delta t\right)&=v\left(t\right)-\omega ^{2}x\left(t\right)\Delta t\\x\left(t+\Delta t\right)&=x\left(t\right)+v\left(t\right)\Delta t\end{aligned}}$$

Em notação matricial temos:

$${\begin{aligned}\left({\begin{array}{c}x\left(t+\Delta t\right)\\v\left(t+\Delta t\right)\end{array}}\right)&=\left({\begin{array}{cc}1&\Delta t\\-\omega ^{2}\Delta t&1\end{array}}\right)\left({\begin{array}{c}x\left(t\right)\\v\left(t\right)\end{array}}\right)\\{\boldsymbol {u}}\left(t+\Delta t\right)&={\overline {M}}{\boldsymbol {u}}\left(t\right)\end{aligned}}$$

Porém a matriz ${\textstyle {\overline {M}}}$ transforma o vetor ${\textstyle {\boldsymbol {u}}\left(t\right)}$ no vetor ${\textstyle {\boldsymbol {u}}\left(t+\Delta t\right)}$, representando então a evolução no espao de fases e seu determinante representa a variação dovolume no espaço de fases. Para um problema conservativo, logo o determinante deve ser ${\textstyle 1}$, uma vez que essevolume deve se manter constante. Para o método de Euler temos:

$$\det \left({\overline {M}}\right)=1+\left(\omega \Delta t\right)^{2}>1$$

Outra forma de analisar o caso da oscilção quando usado o método explícito de Euler, é abrindo as contas. Escrevendo então a energia como:

$${\begin{aligned}E\left(t+\Delta t\right)&={\frac {1}{2}}mv^{2}\left(\Delta t+t\right)+{\frac {1}{2}}kx^{2}\left(\Delta t+t\right)\\&={\frac {1}{2}}m\left(v^{2}\left(\Delta t+t\right)+\omega ^{2}x^{2}\left(\Delta t+t\right)\right)\end{aligned}}$$

Onde fazemos ${\textstyle \omega ^{2}=k/m}$. Usando então:

$${\begin{aligned}v\left(\Delta t+t\right)&=v\left(t\right)-\omega ^{2}x\left(t\right)\Delta t\\x\left(\Delta t+t\right)&=x\left(t\right)+v\left(t\right)\Delta t\end{aligned}}$$

Temos:

$${\begin{aligned}E\left(t+\Delta t\right)&={\frac {1}{2}}m\left(\left(v\left(t\right)-\omega ^{2}x\left(t\right)\Delta t\right)^{2}+\omega ^{2}\left(x\left(t\right)+v\left(t\right)\Delta t\right)^{2}\right)\\&={\frac {1}{2}}m\left[v^{2}\left(t\right)+\omega ^{4}x^{2}\left(t\right)\Delta t^{2}-2v\left(t\right)\omega ^{2}x\left(t\right)\Delta t\right.\\&\left.+\omega ^{2}x^{2}\left(t\right)+\omega ^{2}v^{2}\left(t\right)\Delta t^{2}+2v\left(t\right)\omega ^{2}x\left(t\right)\Delta t\right]\\&={\frac {1}{2}}m\left[\left(v^{2}\left(t\right)+\omega ^{2}x^{2}\left(t\right)\right)+\omega ^{2}\left(v^{2}\left(t\right)+\omega ^{2}x^{2}\left(t\right)+\right)\Delta t^{2}\right.\\&\left.v\left(t\right)\omega ^{2}x\left(t\right)\Delta t-v\left(t\right)\omega ^{2}x\left(t\right)\Delta t\right]\end{aligned}}$$

FIcamos então apenas:

$${\begin{aligned}E\left(t+\Delta t\right)&=E\left(t\right)+E\left(t\right)\omega ^{2}\Delta t^{2}\end{aligned}}$$

Ou ainda:

$${\begin{aligned}E\left(t+\Delta t\right)&=E\left(t\right)\left(1+\omega ^{2}\Delta t^{2}\right)\end{aligned}}$$

Então a cada passo, a energia aumenta com um fator ${\textstyle \left(1+\omega ^{2}\Delta t^{2}\right)}$.

O método de Euler-Crome propõe usar ${\textstyle v\left(t+\Delta t\right)}$ no lugar de ${\textstyle v\left(t\right)}$ para calcular ${\textstyle x\left(t+\Delta t\right)}$. Manipulando temos, lembrando que podemos substituir o valor de ${\textstyle v\left(t+\Delta t\right)}$:

$${\begin{aligned}x\left(t+\Delta t\right)&=x\left(t\right)+v\left(t\right)\Delta t\\&=x\left(t\right)+v\left(t+\Delta t\right)\Delta t\\&=x\left(t\right)+\left[v\left(t\right)-\omega ^{2}x\left(t\right)\Delta t\right]\Delta t\\&=x\left(t\right)\left(1-\omega ^{2}\Delta t^{2}\right)+v\left(t\right)\Delta t\end{aligned}}$$

Atualizando então a notação matricial temos:

$${\begin{aligned}\left({\begin{array}{c}x\left(t+\Delta t\right)\\v\left(t+\Delta t\right)\end{array}}\right)&=\left({\begin{array}{cc}1-\omega ^{2}\Delta t^{2}&\Delta t\\-\omega ^{2}\Delta t&1\end{array}}\right)\left({\begin{array}{c}x\left(t\right)\\v\left(t\right)\end{array}}\right)\\{\boldsymbol {u}}\left(t+\Delta t\right)&={\overline {M}}{\boldsymbol {u}}\left(t\right)\end{aligned}}$$

Calculando então o novo determinante, temos:

$$\det \left({\overline {M}}\right)=1-\omega ^{2}\Delta t^{2}+\left(\omega \Delta t\right)^{2}=1$$

Algumas observações que podem ser feitas: a primeira é que também podemos fazer diferente e usar ${\textstyle x\left(t+\Delta t\right)}$ no lugar de ${\textstyle x\left(t\right)}$ para calcular ${\textstyle v\left(t+\Delta t\right)}$. E a segunda é que quando olhamos para nossa aproximação, temos um intervalo de tempo ${\textstyle \Delta t}$ entre ${\textstyle v\left(t+\Delta t\right)}$ e ${\textstyle v\left(t\right)}$. No método de Euler original, usamo a velocidade no começo intervalo (${\textstyle v\left(t\right)}$) para calcular a nova posição (${\textstyle x\left(t+\Delta t\right)}$, no de Euler-Cramer usamos no fim do intervalo (${\textstyle v\left(t+\Delta t\right)}$), mas de certa forma tem a mesma natureza de aproximação. Como para uma equação tivemos o método de Euler-implícito, porém agora trabalhamos com um sistema de equações. Esse método também é chamado de ’semi-implícito.

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 ; w=w2**(1/2)
#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
for it  in range(Np):
  x.append(x[it]+dt*v[it])  
  v.append(v[it]-dt*x[it+1]*w2) #Usamos x[it+1] ao invés de x[it]
  E.append(k*x[it+1]**2/2+m*v[it+1]**2/2)
  t.append(dt+it*dt)

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