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$$

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{array}{rl}E(t+\mathrm{\Delta }t)& =\frac{1}{2}m{v}^2(\mathrm{\Delta }t+t)+\frac{1}{2}k{x}^2(\mathrm{\Delta }t+t)\\ & =\frac{1}{2}m(v^2(\mathrm{\Delta }t+t)+{\omega }^2{x}^2(\mathrm{\Delta }t+t))\end{array}$$

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

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

Temos:

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

FIcamos então apenas:

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

Ou ainda:

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

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

O mé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)