Método de Euler

De ComplexWiki

Tabela de conteúdo

O problema e o método

Dado um problema cuja dinâmica obedece uma equação diferencial ordinária e condição inicial da forma:

\frac{dx(t)}{dt}=f(x(t),t);\qquad x(t_0) = x_0 \quad (1)

queremos encontrar a solução x(t). O método de Euler parte de fazer uma expansão em série de Taylor (truncada em primeira ordem) da variável incógnita para obter:

x(t+\Delta t) = x(t) + \frac{dx(t)}{dt} \Delta t + ...

e usando a informação contida na equação diferencial obtemos a seguinte relação:

x(t+\Delta t) \simeq x(t) + f(x(t),t) \Delta t \quad (2)

a que pode, mediante iteração a partir da condição inicial x(t0) = x0, nos dar x(t) a qualquer tempo futuro t_1=t_0+\Delta t, \quad t_2=t_1+\Delta t, ...,\quad t_{n+1}=t_n+\Delta t

x_1 \equiv x(t_0+\Delta t) = x_0 + f(x_0,t_0) \Delta t\,

x_2 \equiv x(t_1+\Delta t) = x_1 + f(x_1,t_1) \Delta t\,

\cdots\,

x_{n+1} \equiv x(t_n+\Delta t) = x_n + f(x_n,t_n) \Delta t\,


Exemplo 1: decaimento radiativo

Teoria

Como exemplo de aplicação do método de Euler vamos considerar o decaimento radioativo, que é um problema que tem solução analítica simples. A ideia é, por um lado ganhar confianza com o método estudando uma equação simples, e ao mesmo tempo poder avaliar o erro cometido na integração numérica.

Trata-se de contabilizar o número ou fração de núcleos que vão ficando conforme avanza o tempo a partir de um certo número inicial. A equação diferencial que rege o decaimento radioativo é dada por:

\frac{dN}{N}=-\lambda dt \quad (3)

onde N é o número de núcleos radioativos no instante t e λ é uma constante, diferente para cada elemento. A interpretação física desta equação é simples: a perda relativa de núcleos dN / N é proporcional ao tempo transcorrido (dt), sempre que dt seja pequeno. A constante de proporcionalidade λ é a taxa do decaimento, mede a "velocidade" do decaimento.


Solução analítica

Como antecipamos, e também é fácil de ver da própria equação, a solução é simples pois procuramos por uma função cuja derivada seja proporcional a ela mesma A exponencial é justamente a função que quando derivada volta a aparecer: No caso em questão e usando a condição inicial a solução para N(t) resulta:

N(t) = N_0 \exp(-\lambda t), \quad ou \quad N(t) = N_0 \exp(-ln(2)t/\tau)

onde \quad \tau = ln(2)/\lambda é o tempo necessário para a quantidade de núcleos cair até a metade: \quad N(\tau)=N_0/2

Solução numérica

Rescrevemos primeiro a equação do decaimento radiativo nos termos da equação (1):

\frac{dN(t)}{dt}=-\lambda N(t), \quad N(0) = N_0

para agora aplicar o método de Euler assim:

N_1 \equiv N(\Delta t) = N_0 - \lambda N_0 \Delta t\,

N_2 \equiv N(2\Delta t) = N_1 - \lambda N_1 \Delta t\,

\cdots\,

N_{n+1} \equiv N((n+1)\Delta t) = N_n - \lambda N_n \Delta t\,

Neste caso particular chegamos então a seguinte expressão geral

N_{n+1} = N_n (1 - \lambda \Delta t) \quad (4)


Programa

A expressão (4) está numa forma muito conveniente para se programar o método no computador. Em FORTRAN por exemplo fica assim: 

!Decaimento Radioativo: dN = -lam*N*dt.  Método de Euler
Implicit None
Integer i, Nstep
Real    N, lam, tau, dt, tf

Read* N, tau, dt, tf
Nstep = tf/dt;  lam = log(2.)/tau

Print*, 0., N
Do i = 1, Nstep

   N = (1.0 - lam*dt)*N     !Notar que esta é a forma de implementar (4) em FORTRAN

   Print*, i*dt, N
End Do
End


Método de Euler nas equações de Newton

A equação de Newton, de um modo geral, trata-se de uma equação diferencial ordinâria de segunda ordem. O que quer dizer que aparece uma derivada segunda. Porem, o método de Euler, como fora formulado antes, envolve a primeira derivada. O problema se resolve re-escrevendo a equação de Newton como duas equações de primeira ordem. Vejamos:

\frac{d^2x}{dt^2} \equiv a = f/m

passa para:

\frac{dv}{dt} \equiv a; \qquad \frac{dx}{dt} \equiv v

Ou seja de uma equação de segunda ordem passamos a duas de primeira ordem representando duas formas equivalentes do mesmo problema.


Exemplo 2: oscilador harmônico

O oscilador harmônico é um problema clássico da mecânica Newtoniana e seu equivalente linear é o sistema massa-mola; um bloco de massa m ligado à parede por uma mola de constante elástica k.


Teoria

A equação de Newton para este sistema fica: m\frac{d^2x(t)}{dt^2} = -k x(t).

Solução analítica

Definindo \omega = \sqrt{k/m} a solução do problema está dada por:

x(t) = x_0 \sin({\omega t + \phi})\,, x(0) e φ saem das condiçãoes iniciais: x(0) = x_0 \sin(\phi)\, e \dot x(0) = \omega x_0 \cos (\phi)\,


Solução numérica

Como no Exemplo 1 usamos a expansão de Taylor a primeira ordem para escrever:

v(t+\Delta t) = v(t) + a(t) \Delta t \,

x(t+\Delta t) = x(t) + v(t) \Delta t \,

Essa a forma geral do Método ou Algoritmo de Euler para a equação de Newton, onde em cada passo de iteração tanto a posição quanto a velocidade devem ser atulizadas para gerar uma seqüência de x e v nos tempos discretos a partir de t = 0 e em passos Δt. No Exemplo 2 (o oscilador armônico) é só sustituir a(t)\, por -\omega^2 x(t)\,, a aceleração do bloco, e temos:

v(t+\Delta t) = v(t) - \omega^2 x(t) \Delta t \,

x(t+\Delta t) = x(t) + v(t) \Delta t \,

Programa

A base do código FORTRAN é a mesma do Exemplo anterior, a única diferença está no laço central envolvendo agora duas variáveis. Ha também um detalhe técnico justamente pelo fato de trabalhar agora com duas grandezas (posição e velocidade), em lugar de só uma como foi o N no Exemplo 1. A parte específica deste problema em FORTRAN ficaria a priori: 

...
v = v - omega2 * x * dt
x = x + v * dt
...

Acontece que ao fazer assim, na primeira linha o v passa de v(t) para v(t + Δt), sendo este último o usado na linha de atualização do x. Mas o algoritmo de Euler indica que nesse momento x deve usar v(t) para ser atualizado e não v(t + Δt). Este é o problema técnico que assinalamos antes. Como resolvé-lo na prática? Usando uma variável auxiliar onde guardamos o valor anterior, assim: 

...
vold = v
v = v - omega2 * x * dt
x = x + vold * dt
...

Desta forma x usa o v no tempo correto para fazer a sua atualização. Uma forma alternativa de codificar (um pouco mais elegante do posto de vista físico) seria: 

...
a = -omega2 * x
x = x + v * dt
v = v + a * dt
...

Onde ha duas alterações: a troca entre as linhas de x e v e o uso de variável auxiliar a (aceleração) para que v use a aceleração no tempo correto. È importante que se convençam da equivalencia dos dois códigos. Agora é so implementar o código completo para observar a solução numérica do oscilador armônico. Neste caso teremos uma surpresa, essa solução resulta diferente do esperado, o que nos fará pensar: O que ha de errado?


Erro do algoritmo

Integrando numericamente com o algoritmo de Euler o problema do oscilador armônico, se observa claramente que a amplitude e a velociade aumentam continuamante, quando deveriam oscilar entre valores bem definidos. Em fim, que a energia do sistema aumenta quando deveria ser estritamente constante.

Isso pode ser mostrado analiticamante escrevendo a evolução do sistema em termos matriciais, assim:


{x' \choose v'} = \begin{pmatrix} 1 & \delta t \\ -\omega^2 \delta t & 1 \end{pmatrix} {x \choose v}\,


O aumento da energia está representado pelo determinante da matriz de evolução:

||det|| = 1 + \omega^2 \delta t^2\,


quando aplicada n vezes para ir de \quad t=0 até \quad t = n\delta t provoca um aumento na energia que vai como e^{(\omega^2 \delta t) t}.

Ou seja, cresce exponencialmente, mesmo que o \quad \delta t seja muito pequeno. O \quad \delta t modifica o pre-fator (o que multiplica t na exponencial) porem, mais cedo o mais tarde, (no caso de um passo de tempo menor) sempre aumenta. Todo isso pode ser conferido numericamente estudando a evolução da energia com o tempo.

Obtido em "http://davinci.if.ufrgs.br/wiki/index.php/M%C3%A9todo_de_Euler"
Ferramentas pessoais