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{array}{rl}\frac{dv}{dt}& =-{\omega }^{2}x\left(t\right)\\ \frac{dx}{dt}& =v\left(t\right)\end{array}$$

Aplicando o método de Euler então:

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

Em notação matricial temos:

$$\begin{array}{rl}\left(\begin{array}{c}x(t+\mathrm{\Delta }t)\\ v(t+\mathrm{\Delta }t)\end{array}\right)& =\left(\begin{array}{cc}1& \mathrm{\Delta }t\\ -{\omega }^2\mathrm{\Delta }t& 1\end{array}\right)\left(\begin{array}{c}x\left(t\right)\\ v\left(t\right)\end{array}\right)\\ \mathit{u}(t+\mathrm{\Delta }t)& =\bar{M}\mathit{u}\left(t\right)\end{array}$$

Porém a matriz $\bar{M}$ transforma o vetor $\mathit{u}\left(t\right)$ no vetor $\mathit{u}(t+\mathrm{\Delta }t)$, 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 $1$, uma vez que essevolume deve se manter constante. Para o método de Euler temos:

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

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

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

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

$$\begin{array}{rl}\left(\begin{array}{c}x(t+\mathrm{\Delta }t)\\ v(t+\mathrm{\Delta }t)\end{array}\right)& =\left(\begin{array}{cc}1-{\omega }^2\mathrm{\Delta }{t}^2& \mathrm{\Delta }t\\ -{\omega }^2\mathrm{\Delta }t& 1\end{array}\right)\left(\begin{array}{c}x\left(t\right)\\ v\left(t\right)\end{array}\right)\\ \mathit{u}(t+\mathrm{\Delta }t)& =\bar{M}\mathit{u}\left(t\right)\end{array}$$

Calculando então o novo determinante, temos:

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

Algumas observações que podem ser feitas: a primeira é que também podemos fazer diferente e usar $x(t+\mathrm{\Delta }t)$ no lugar de $x\left(t\right)$ para calcular $v(t+\mathrm{\Delta }t)$. E a segunda é que quando olhamos para nossa aproximação, temos um intervalo de tempo $\mathrm{\Delta }t$ entre $v(t+\mathrm{\Delta }t)$ e $v\left(t\right)$. No método de Euler original, usamo a velocidade no começo intervalo ($v\left(t\right)$) para calcular a nova posição ($x(t+\mathrm{\Delta }t)$, no de Euler-Cramer usamos no fim do intervalo ($v(t+\mathrm{\Delta }t)$), 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)