Lembrando do que vimos no Método de Euler, o sistema de equações para o sistema massa-mola era:
![{\displaystyle {\begin{aligned}{\frac {dv}{dt}}&=-\omega ^{2}x\left(t\right)\\{\frac {dx}{dt}}&=v\left(t\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/67f17e5ca40758bf64617048cd5ff246b5d280dd)
Aplicando o método de Euler então:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/667aa006ded2f1459068fab294dc4954b1fcaf56)
Em notação matricial temos:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/82bc138838c6ef2f7e19b0dbbd90e93dda23de33)
Porém a matriz
transforma o vetor
no vetor
, 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
, uma vez que essevolume deve se manter constante. Para o método de Euler temos:
![{\displaystyle \det \left({\overline {M}}\right)=1+\left(\omega \Delta t\right)^{2}>1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/134653a1791bfdbe98339251af02521b64927f34)
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:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/64a7421f2442f662a95b619d041db253c2d80191)
Onde fazemos
. Usando então:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/be2c2f9c256604f6a4a555d32365530d456381b0)
Temos:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8f0fd807a412052c0652711e82a21365f13c713b)
FIcamos então apenas:
![{\displaystyle {\begin{aligned}E\left(t+\Delta t\right)&=E\left(t\right)+E\left(t\right)\omega ^{2}\Delta t^{2}\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a8c36618d79905d9fe7b0de3882f21662efe793a)
Ou ainda:
![{\displaystyle {\begin{aligned}E\left(t+\Delta t\right)&=E\left(t\right)\left(1+\omega ^{2}\Delta t^{2}\right)\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e6b93cfdfa938582d4b09b4dcb5c68c9bf7adf1d)
Então a cada passo, a energia aumenta com um fator
.
O método de Euler-Crome propõe usar
no lugar de
para calcular
. Manipulando temos, lembrando que podemos substituir o valor de
:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/fc86afcc2bd2549f50c03b98bf1177984cbcaf39)
Atualizando então a notação matricial temos:
![{\displaystyle {\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6e235f60c434ea62d17da180a65fd2cd46719360)
Calculando então o novo determinante, temos:
![{\displaystyle \det \left({\overline {M}}\right)=1-\omega ^{2}\Delta t^{2}+\left(\omega \Delta t\right)^{2}=1}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8b4ddabcf200b0c9a82970ac0315384507cf185a)
Algumas observações que podem ser feitas: a primeira é que também podemos fazer diferente e usar
no lugar de
para calcular
. E a segunda é que quando olhamos para nossa aproximação, temos um intervalo de tempo
entre
e
. No método de Euler original, usamo a velocidade no começo intervalo (
) para calcular a nova posição (
, no de Euler-Cramer usamos no fim do intervalo (
), 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
#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-Cromer
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)