5. Método de Euler

Resolução numérica de equações diferenciais

A resolução numérica de uma equação de primeira ordem:

eq-aula5-1.png

consiste em calcular o valor da variável de estado numa sequência discreta de instantes

eq-aula5-2.png

usando alguma estimativa do valor médio das derivada durante cada intervalo de tempo

eq-aula5-3.png

a partir da função eq-aula5-4.png, que é a derivada instantânea. Podemos usar uma sequência de instantes igualmente espaçados entre si, com incremento de tempo h:

eq-aula5-5.png

assim, substituiremos a variável contínua eq-aula5-6.png por uma variável discreta:

eq-aula5-7.png

Método de Euler

Por definição, a derivada eq-aula5-8.png é

eq-aula5-9.png

se eq-aula5-10.png, então eq-aula5-11.png. Admitindo que eq-aula5-12.png é suficientemente pequeno, obtemos a aproximação

eq-aula5-13.png

e substituindo na equação do sistema:

eq-aula5-14.png

A condição inicial eq-aula5-15.png permite calcular sequencialmente eq-aula5-16.png, eq-aula5-17.png, etc.

Do ponto de vista gráfico, o que estamos a fazer é deslocarmos-nos, desde o ponto eq-aula5-18.png, uma distância eq-aula5-12.png, segundo o eixo eq-aula5-19.png, na direcção do campo de direcções:

fig-5-1.png

Como mostra a figura, a direcção do campo no ponto i já não é a mesma no ponto i+1 e, assim, a curva obtida não segue perfeitamente o campo de direcções. Mas se h for suficiente pequena, obtém-se uma boa aproximação.

Exemplo 1

Usando o método de Euler, encontre a carga em função do tempo, no circuito RC com equação (carga em micro-coulombs e tempo em mili-segundos)

eq-aula5-20.png

durante os primeiros 6 ms, se a carga inicial for nula.

Resolução: Convém guardar os resultados em listas para poder fazer gráficos. Vamos começar por usar incrementos de tempo de 0.1 ms. Teremos que iniciar 3 listas para as variáveis t, Q e a derivada de Q:

(%i1) h1: 0.1$

(%i2) t1: [0]$

(%i3) Q1: [0]$

(%i4) Qt1: [1.25]$

Será utilizado o comando endcons para acrescentar números no fim de cada lista. Para avançar 6 ms, precisamos 60 iterações:

(%i5) for n:1 thru 60 do
          (t1: endcons(last(t1) + h1, t1),
           Q1: endcons(last(Q1) + h1*last(Qt1), Q1),
           Qt1: endcons(1.25 - last(Q1), Qt1))$

Podemos ver o último elemento de cada lista, usando o comando last

(%i6) last(t1);
(%o6)       5.999999999999995

(%i7) last(Q1);
(%o7)       1.247753737125107

(%i8) last(Qt1);
(%o8)       .002246262874892935

Convém repetir o processo com um valor menor de h e comparar os resultados:

(%i9) h2: 0.01$

(%i10) t2: [0]$

(%i11) Q2: [0]$

(%i12) Qt2: [1.25]$

(%i13) for n:1 thru 600 do
          (t2: endcons(last(t2) + h2, t2),
           Q2: endcons(last(Q2) + h2*last(Qt2), Q2),
           Qt2: endcons(1.25 - last(Q2), Qt2))$

(%i14) last(Q2);
(%o14)            1.246993738385861

Este problema pode ser resolvido analiticamente. A solução exacta é:

eq-aula5-21.png

Para comparar com os resultados do método de Euler, vamos criar uma lista com a solução exacta:

(%i15) Q3: makelist(float(1.25*(1  - exp(-n*h2))), n, 0, 600)$

O gráfico pode ser feito com a função graph2d (disponível em http://fisica.fe.up.pt/pub/maxima/)

(%i16) load("graph2d")$

(%i17) graph2d([label,"h=0.1"],  t1, Q1,
               [label,"h=0.01"], t2, Q2,
               [label,"exacta"], t2, Q3);
fig-5-2.png

A solução obtida com h=0.01 já parece idêntica à solução exacta.

Pode ser conveniente também criar tabelas de dados a partir dos resultados. Por exemplo, os resultados das últimas 5 iterações:

(%i18) for n:597 thru 601
          do print(t2[n], Q2[n], Qt2[n])$
5.959999999999917 1.246870420465166 .003129579534833837
5.969999999999917 1.246901716260514 .003098283739485508
5.979999999999917 1.246932699097909 0.00306730090209073
5.989999999999917 1.24696337210693 .003036627893069799
5.999999999999917 1.246993738385861 .003006261614139083

É possível também guardar os resultados num ficheiro, para serem usados com outros programas. Por exemplo, para guardar as 3 listas num ficheiro dados.txt, usaremos o comando with_stdout:

(%i19) display2d: false$

(%i20) with_stdout("dados.txt",
         for n thru 601
             do print(t2[n], Q2[n], Qt2[n]))$

foi dado um valor falso à variável display2d para que não sejam escritos espaços entre os sinais negativos e os números.

Um tipo de ficheiro de dados bastante útil é o designado por CSV (Comma Separated Values), usualmente, que pode ser importado pelos programas de folha de cálculo


(%i31) with_stdout("dados.csv",
         for n thru 601
             do print(t2[n], ",", Q2[n], ",", Qt2[n]))$

Jaime Villate
Data: 2009-10-20 19:03:50 +0000