Método de Euler: mudanças entre as edições

De Física Computacional
Ir para navegação Ir para pesquisar
Sem resumo de edição
Sem resumo de edição
 
(14 revisões intermediárias pelo mesmo usuário não estão sendo mostradas)
Linha 39: Linha 39:


<span id="exemplo-1"></span>
<span id="exemplo-1"></span>
== Exemplo ==
== Exemplo: Decaimento ==


O primeiro exemplo de aplicação é o '''decaimento radiativo''', cuja equação diferencial é: <math display="block">\frac{dN}{dt}=-\lambda N</math>
O primeiro exemplo de aplicação é o '''decaimento radiativo''', cuja equação diferencial é: <math display="block">\frac{dN}{dt}=-\lambda N</math>
Linha 63: Linha 63:
<pre>
<pre>
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
import numpy as np                #Biblioteca com funções matemáticas
N =[10**6];nt=6;lam=0.1;dt=0.1    #Parâmetros
N =[10**6];nt=6;lam=0.1;dt=0.1    #Parâmetros


Linha 92: Linha 93:


<span id="exemplo"></span>
<span id="exemplo"></span>
== Exemplo ==
== Exemplo: Decaimento ==


Trabalhando novamente com o decaimento radioativo:
Trabalhando novamente com o decaimento radioativo:
Linha 125: Linha 126:
<pre>
<pre>
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
import numpy as np                #Biblioteca com funções matemáticas
N =[10**6];nt=6;lam=0.1;dt=0.1    #Parâmetros
N =[10**6];nt=6;lam=0.1;dt=0.1    #Parâmetros


Linha 138: Linha 140:


= Avaliação de erro =
= Avaliação de erro =
Para o cálculo de erros deste exemplo, utilizamos bastante a meia-vida como um valor caratectístico do sistema, principalmente para nos basearmos até definir até onde a simulação vai rodar. A meia vida é simplesmente o tempo necessário para a quantidade de núcleos cair para a metade do valor inicial. Isto é, primeiro temos que a solução analítica é simplesmente:
<math display="block">N=N_{0}e^{-\lambda t}</math>
Então a meia vida é quando <math display="inline">t=\tau</math>:
<math display="block">\frac{N_{0}}{2}=N_{0}e^{-\lambda\tau}\longrightarrow\tau=\frac{\ln\left(2\right)}{\lambda}</math>


== Erros computados ==  
== Erros computados ==  
Linha 143: Linha 153:


<pre>
<pre>
import numpy as np                #Biblioteca com recursos matemáticos
import numpy as np                #Biblioteca com recursos matemáticos
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
Linha 154: Linha 163:
T  = [0]                        #Lista de tempos no qual vamos calcular N
T  = [0]                        #Lista de tempos no qual vamos calcular N


tau=np.log(2)/lam                   #Tempo final em que vamos obter o erro
tau=np.log(2)/lam                 #Tempo final em que vamos obter o erro


#Solução numérica
#Solução numérica
Linha 163: Linha 172:


#Solução analítica
#Solução analítica
t = np.arange(0.0,nt*tau, 1)       #Tempo
t = np.arange(0.0,nt*tau, 1)     #Tempo
NA =N0*np.exp(-lam*t)            #Solução
NA =N0*np.exp(-lam*t)            #Solução


Linha 179: Linha 188:
Erro 2: Curva de erro pra diferentes dt
Erro 2: Curva de erro pra diferentes dt


<pre>import numpy as np                #Biblioteca com recursos matemáticos
<pre>
import numpy as np                #Biblioteca com recursos matemáticos
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
import matplotlib.pyplot as plt  #Biblioteca para plotar gráficos
N0 =10**6;lam=0.001;nt=6          #Parâmetros
N0 =10**6;lam=0.001;nt=6          #Parâmetros
Linha 195: Linha 205:
   T=[0]            #Tempo
   T=[0]            #Tempo
   frac = (1-lam*dt) #Funçao calculada
   frac = (1-lam*dt) #Funçao calculada
   tau=np.log(2)/lam                   #Meia-vida
   tau=np.log(2)/lam                   #Meia-vida
   for i in range(1,int(nt*tau/dt)+1):   #Vamo calcular um múltiplo da meia vida
   for i in range(1,int(nt*tau/dt)+1): #Vamo calcular um múltiplo da meia vida
     N.append(frac*N[i-1])              #Solução numérica
     N.append(frac*N[i-1])              #Solução numérica
     T.append(i*dt)                    #Tempo no qual achamos a solução
     T.append(i*dt)                    #Tempo no qual achamos a solução
     err=abs(N[i]-N0*np.exp(-lam*T[i])) #Erro
     err=abs(N[i]-N0*np.exp(-lam*T[i])) #Erro
    print(N[i],N0*np.exp(-lam*T[i]))
     D.append(err)                      #Guardamos o erro
     D.append(err)                      #Guardamos o erro
   return(D,T)                          #Retorno o erro
   return(D,T)                          #Retorno o erro
Linha 216: Linha 225:
plt.xlabel('Tempo [s]');plt.ylabel('Erro[abs(N-Nmétodo)]')            #Legenda dos eixos
plt.xlabel('Tempo [s]');plt.ylabel('Erro[abs(N-Nmétodo)]')            #Legenda dos eixos
plt.ticklabel_format(style='plain')                                    #Desativar a notação científica
plt.ticklabel_format(style='plain')                                    #Desativar a notação científica
plt.show()                         
plt.show()                                
</pre>
 
Erro 3: Erro em relação a dt com um t fixo
<pre>
import numpy as np                        #Biblioteca com recursos matemáticos
import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
N0 =10**6;lam=0.001                       #Parâmetros
dt=[1/10**(i) for i in np.arange(2.9,-0.1,-0.3)]
tfin=[(1,'ob','-k'),(4,'sr','-m')]        #Tuplas com os tempos finais e os marcadores principais e do fit
y="log";x="log"                            #Se os eixos vão ser "linear" ou "log"
 
def sol_metodo(dtT,N0,lam,t):
  """Função para calcular a solução de acordo com o método de Euler e
  retornar o erro comparado a solução analítica"""
  #dtT  - Intervalo do passo em unidade de tau
  #N0  - Número inicial de Núcleos
  #lam  - Lambda
  #t    - Tempo em que vamos pegar a solução, múltiplo da meia-vida
  N=N0
  T=np.log(2)/lam                  #Tempo final em que vamos obter o erro
  dt=dtT*T                        #Intervalo de tempo em segundos
  frac = (1-lam*dt)                #Funçao calculada
  for i in range(1,int(t*T/dt)+1): #Vamo calcular os passos necessários para 4000s
    N=frac*N                      #Solução numérica
  err =abs(N-N0*np.exp(-lam*i*dt)) #Erro
  return(err)                      #Retorno o erro
 
for t in tfin:                                  #Percorremos todos os tempos finais
  E=[]                                          #Lista para guardar todos os erros com
  for d in dt:
    E.append (sol_metodo(d,N0,lam,t[0]))              #Computamos o erro para cada dt
  plt.plot(dt,E,t[1], label='t = '+str(t[0])+' τ')    #Plot dos marcadores
  plt.plot(dt,E,'--'+t[1][1])                          #Plot Tracejado
  m, b = np.polyfit(np.log(np.array(dt)),np.log(np.array(E)), 1) #Uma curva para fitar os logaritmos
  label='ln(erro) = {:.2f} + {:.2f} ln(Δt)'.format(b,m)
  plt.plot(dt, np.e**b*np.array(dt)**m,t[2],label=label)        #Tiramos a exponencial do logaritmo
 
#Configuramos o gráfico
plt.title('Erro X dt')                              #Título
plt.legend()                                        #Legenda das curvas
plt.xlabel('Δt  [τ]')                                #Legenda dos eixos
plt.ylabel('Erro[abs(N-Nmétodo)]') 
plt.ticklabel_format(style='plain')                  #Desativar a notação científica
plt.xscale(x);plt.yscale(y)                          #Escala
plt.show()                                          #Exibir o resultado
</pre>
</pre>
== Erro teórico ==
== Erro teórico ==
Outro caminho que podemos tomar para obtermos o método de Euler, é através da expansão em série de Taylor de uma função <math display="inline">y\left(t+\Delta t\right)</math> em torno do ponto <math display="inline">y\left(t\right)</math>. Lembrando que a expansão em série de Taylor de uma função <math display="inline">f\left(x\right)</math> em torno de <math display="inline">a</math> é dada por:
<math display="block">f\left(x\right)=\sum_{n=0}^{\infty}\frac{f^{\left(n\right)}}{n!}\left(a\right)\left(x-a\right)^{n}</math>
Então fazendo esta expansão e escrevendo <math display="inline">t'=t+\Delta t</math>
<math display="block">\begin{align}
y\left(t'\right) & =\sum_{n=0}^{\infty}\frac{y^{\left(n\right)}}{n!}\left(t\right)\left(t'-t\right)^{n}\\
& =y\left(t\right)+\Delta ty'\left(t\right)+\sum_{n=2}^{\infty}\frac{y^{\left(n\right)}}{n!}\left(t\right)\Delta^{n}\end{align}</math>
Podemos olhar <math display="inline">f^{\left(n\right)}\left(a\right)</math> com mais atenção uma vez que tecnicamente nosssa variável da função que estamos expandindo é <math display="inline">t'</math>.
<math display="block">\frac{d}{dt'}y\left(t'\right)\approx\frac{y\left(t'+\Delta t'\right)-y\left(t'\right)}{\Delta t'}</math>
Considerando que <math display="inline">\Delta t</math> é constante, então: <math display="inline">\Delta t'=\Delta t</math>, e sendo o ponto <math display="inline">a=t</math>, temos:
<math display="block">\frac{d}{dt'}y\left(t'\right)\approx\frac{y\left(t'+\Delta t\right)-y\left(t'\right)}{\Delta t}=\frac{d}{dt}f\left(t'\right)\longrightarrow\left.\frac{d}{dt}y\left(t'\right)\right|_{t'=t}=y'\left(t\right)</math>
De forma que ficamos então apenas com:
<math display="block">\begin{align}
y\left(t+\Delta t\right) & =y\left(t\right)+f\left(t ,y \right) \Delta t+\sum_{n=2}^{\infty}\frac{y^{\left(n\right)}\left(t\right)}{n!}\Delta t^{n}\end{align}</math>
Os primeiros dois termos então correspondem ao método de Euler, e o somatório se torna o erro ao jogarmos fora. O menor grau que desprezamos é <math display="inline">n=2</math>, então o erro local é <math display="inline">Err\left(\Delta t\right)\propto\Delta t^{2}</math>. E como o número de passos é dado por <math display="inline">N=t/\Delta t</math> então após <math display="inline">N</math> passos temos uma estimativa de erro global dada por <math display="inline">Err\left(\Delta t,t\right)\propto N \cdot \Delta t^{2} </math> ou apenas <math display="inline">Err\left(\Delta t,t\right)\propto\Delta t</math>. Cálculos mais detalhados da estimativa do erro podem ser encontrados em textos matemáticos<ref>[https://math.libretexts.org/Courses/Monroe_Community_College/MTH_225_Differential_Equations/3%3A_Numerical_Methods/3.1%3A_Euler%27s_Method  Euler's Method] (William F. Trench, LibreTexts)</ref>.
= Exemplo: SIR =
O modelo SIR é um modelo de epidemia compartimental no qual separamos a população em três compartimentos:
* <math display="inline">S</math>: População suscetível
* <math display="inline">I</math>: População infectada
* <math display="inline">R</math>: População de recuperados
O modelo SIR é descrito então pelo seguinte conjunto de equações:
<math display="block">\begin{align}
\frac{dS}{dt} & =-\beta SI\\
\frac{dI}{dt} & =\beta SI-\gamma I\\
\frac{dR}{dt} & =\gamma I\end{align}</math>
Onde temos as seguintes variáveis:
* <math display="inline">\beta</math>: transmissibilidade
* <math display="inline">\gamma</math>: taxa de recuperação
** É o inverso do tempo que a pessoa permanece infectada <math display="inline">\tau=\frac{1}{\gamma}</math>.
Neste modelo, as pessoas podem se transferir entre os compartimento apenas em um sentido <math display="inline">S\rightarrow I\rightarrow R</math>, e pode-se destacar que <math display="inline">R</math> comporta tanto pessoas que se curaram e tornaram-se imunes, quanto falecimentos, em todo, caso, são pessoas que não podem ser suscetíveis a infecção novamente. Como a população total <math display="inline">N</math> é constante, temos <math display="inline">\dot{N}=0</math> e por conservação:
<math display="block">\dot{S}+\dot{I}+\dot{R}=0</math>
Desta forma só precisamos resolver duas equações a cada passo. Além disso podemos considerar que a transmissibilidade de uma doença depende da população total, de forma que <math display="inline">\beta\left(N\right)=\widehat{\beta}/N</math>, sendo <math display="inline">\widehat{\beta}</math> um parâmetro próprio da doença e independente da população. Substituindo e reescrevendo cada compartimento como uma fração da população total, ao invés de números totais, isto é <math display="inline">x=X/N</math>, seja <math display="inline">x</math> um compartimento qualquer, então podemos reescrever o sistema como:
<math display="block">\begin{align}
\frac{ds}{dt} & =-\widehat{\beta}si\\
\frac{di}{dt} & =\widehat{\beta}si-\gamma i\\
\frac{dr}{dt} & =-\left(\frac{ds}{dt}+\frac{di}{dt}\right)\end{align}</math>
Um ponto interessante é a condição de epidemia. A segunda equação nos dá a variação de infectados, para termos epidemia precisamos que este valor positivo, isto é <math display="inline">\frac{di}{dt}>0</math>, logo:
<math display="block">\begin{align}
\left(\widehat{\beta}s-\gamma\right)i & >0\\
\widehat{\beta}s & >\gamma\\
\widehat{\beta}\tau & >\frac{1}{s}\end{align}</math>
Considerando como condição inicial que <math display="inline">S\approx N</math>, isto é, praticamente toda população é suscetível, logo <math display="inline">s=1</math> e temos então o seguinte limiar <math display="inline">\widehat{\beta}\tau\gtrapprox1</math> para uma condição inicial de epidêma. Podemos definir assim o número reprodutivo básico:
<math display="block">R_{0}\equiv\widehat{\beta}\tau>1</math>
Além disso, um modelo SIRS (no qual recuperados podem tornar-se suscetíveis novamente) e com atraso foi discutido [[Sistemas de equações diferenciais com atrasos fixos - SIRS | aqui]]
<pre>
import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
b  = 0.1  #Transmissibilidade
tau = 11  #Tempo que a pessoa permanece infectada
g  = 1/tau#Taxa de recuperação
N  = 1000  #População total
Np  = 5000  #Número de passos
dt  = 0.5  #Intervalo de tempo
#Valores iniciais
i=[1/N]; s=[1-i[0]]; r=[0]; t=[0]
T=[s[0]+i[0]+r[0]]  #Total
for it  in range(Np):
  s.append(s[it]-dt*b*s[it]*i[it])
  i.append(i[it]+dt*(b*s[it]*i[it]-g*i[it]))
  r.append(1-(i[it+1]+s[it+1]))
  t.append(dt+it*dt)
  T.append(s[it+1]+i[it+1]+r[it+1])
plt.plot(t,i) #Infectados
plt.plot(t,s) #Suscetíveis
plt.plot(t,r) #Recuperados
plt.plot(t,T) #Total
</pre>
= Exemplo: Sistema massa-mola =
Podemos escrever o problema geral da mecânica clássica como:
<math display="block">m\frac{d^{2}x}{dt^{2}}=f\left(x,v,t\right)</math>
O método de Euler é um método para equações de 1ª ordem, então transformando em um sistema de equações de primeira ordem, escrevendo <math display="inline">\frac{dx}{dt}=v</math>, temos:
<math display="block">\begin{align}
m\frac{dv}{dt} & =f\left(x,v,t\right)\\
\frac{dx}{dt} & =v\left(t\right)\end{align}</math>
E pela segunda lei de Newton <math display="inline">f=ma</math>, então ficamos simplesmente com:
<math display="block">\begin{align}
\frac{dv}{dt} & =a\left(t\right)\\
\frac{dx}{dt} & =v\left(t\right)\end{align}</math>
E aplicando o método de Euler:
<math display="block">\begin{align}
v\left(t+\Delta t\right) & =v\left(t\right)+a\left(t\right)\Delta t\\
x\left(t+\Delta t\right) & =x\left(t\right)+v\left(t\right)\Delta t\end{align}</math>
Agora tendo a lei de Hooke:
<math display="block">m\frac{d^{2}x}{dt^{2}}=-kx</math>
E reescrevendo <math display="inline">\omega^{2}=\frac{k}{m}</math>, temo que:
<math display="block">\begin{align}
\frac{dv}{dt} & =-\omega^{2}x\left(t\right)\\
\frac{dx}{dt} & =v\left(t\right)\end{align}</math>
Onde <math display="inline">-\omega^{2}x\left(t\right)=a\left(t\right)</math>. Então o método de Euler fica:
<math display="block">\begin{align}
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{align}</math>
A solução analítica do sistema é dada por:
<math display="block">
\begin{align}
v\left(t\right) & = - \omega x_{0}\sin\left(\omega t+\phi\right)\\
x\left(t\right) & =x_{0}\cos\left(\omega t+\phi\right)
\end{align}
</math>
Onde podemos ver que a energia e conserva:
<math display="block">E=V+T=\frac{1}{2}kx\left(t\right)^{2}+\frac{1}{2}mv\left(t\right)^{2}</math>
Sendo <math display="inline">v\left(0\right)=0</math>, então:
<math display="block">E=\frac{1}{2}kx_{0}^{2}=\text{constante}</math>
<pre>
import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
Np  = 200000              #Número de passos
dt  = 0.001              #Intervalo de tempo
m=1  ; k= 1.; w2= k/m    #Constantes
x=[1]; v=[0]; t=[0]; E=[k*(x[0]**2)/2+m*(v[0]**2)/2] #Valores iniciais
#Método de Euler
for it  in range(Np):
  x.append(x[it]+dt*v[it])
  v.append(v[it]-dt*x[it]*w2)
  E.append(k*x[it+1]**2/2+m*v[it+1]**2/2)
  t.append(dt+it*dt)
plt.plot(t,x) #Posição
plt.plot(t,v) #Velocidade
plt.plot(t,E) #Energia
</pre>
= Exemplo: Modelo de Lotka Volterra =
Além disso temos uma discussão sobre o [[Modelo de Lotka-Volterra]] e também sobre [[Modelo de Lotka-Volterra amortecido]], que além da discussão apresenta uma solução numérica utilizando o Método de Euler, com a única diferença de um termo de amortecimento. Uma vez que zerarmos este termo, retornamos ao Modelo de Lotka-Volterra clássico.
= Citações =
<references />


= Principais materiais utilizados =
= Principais materiais utilizados =


# [https://www.ufrgs.br/reamat/CalculoNumerico/livro-py/pdvi-erro_de_truncamento.html  Erro de truncamento] (REAMAT, UFRGS)
# [https://tutorial.math.lamar.edu/classes/de/eulersmethod.aspx Euler's Method] (Paul Dawkins, Universidade Lamar)
# [https://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node3.html Forward and Backward Euler Methods] (Michael Zeltkevic, Instituto de Tecnologia de Massachusetts)
# [https://web.mit.edu/10.001/Web/Course_Notes/Differential_Equations_Notes/node3.html Forward and Backward Euler Methods] (Michael Zeltkevic, Instituto de Tecnologia de Massachusetts)
# [https://tutorial.math.lamar.edu/classes/de/eulersmethod.aspx Euler's Method] (Paul Dawkins, Universidade Lamar)
# [https://www.ufrgs.br/reamat/CalculoNumerico/livro-py/pdvi-metodo_de_euler.html#x105-19100010.2  Método de Euler] (REAMAT, UFRGS)
# [https://www.ufrgs.br/reamat/CalculoNumerico/livro-py/pdvi-metodo_de_euler.html#x105-19100010.2  Método de Euler] (REAMAT, UFRGS)

Edição atual tal como às 15h16min de 15 de março de 2022

Euler explícito

A grande maioria das equações diferenciais de primeira ordem não podem ser resolvidas analiticamente. Para o comportamento a longo prazo de uma solução podemos tentar esboçar um campo de direções, mas se precisamos conhecer mais especificamente como uma solução se comporta, precisamos de outra ferramenta. Os métodos numéricos nos permitem obter soluções aproximadas para as equações diferenciais.

Começando com uma problema genérico de valor inicial:

Considerando que conhecemos a função e os valores na condição inicial, assumimos que é tudo contínuo de forma que sabemos que uma solução de fato vai existir. Temos então para :

Dessa forma podemos escrever uma reta tangente à curva no ponto usando a inclinação :

Para visualizar melhor esta equação, podemos fazer , ficmos então com . Desta forma, fica ainda mais evidente que esta é uma equação de reta com inclinação , e quando temos , ou seja, uma reta que passa pelo ponto . Para temos apenas um deslocamento no eixo.

Então se é perto o suficiente de , a equação da reta vai estar perto do valor atual da solução em . Então podemos escrever:

Podemos repetir o processo, usando agora como valor inicial, então:

Ou de maneira genérica:

Podemos ainda reescrever o passo como , de forma que ficamos com:

Outra forma de visualizar o resultado, é considerar a reta:

Como a solução aproximada para o intervalo . Então com um conjunto de retas podemos ter uma aproximação para a solução como um todo.

Exemplo: Decaimento

O primeiro exemplo de aplicação é o decaimento radiativo, cuja equação diferencial é:

Onde é a quantidade de partículas que sofrem o decaimento e a taxa no qual o decaimento ocorre.

  • Notem que a mesma equação pode descrever a diminuição de uma população estéril ( sendo a quantidade de indivíduos vivos e a taxa de mortalidade) ou a descarga de um circuito RC.
  • A aplicação do método a este exemplo de primeira ordem nos leva a seguinte relação de recorrência

Ou mais explicitamente:

Implementando:

import matplotlib.pyplot as plt   #Biblioteca para plotar gráficos
import numpy as np                #Biblioteca com funções matemáticas
N =[10**6];nt=6;lam=0.1;dt=0.1    #Parâmetros

tau=np.log(2)/lam                 #Tempo final em que vamos obter o erro
fac = 1-lam*dt                    #Função calculada
for i in range(int(nt*tau/dt)+1): #Vamo calcular um múltiplo da meia vida
  N.append(fac*N[i])              #Salvamos o novo valor
  print(i*dt,N[i])                #printamos o resultado

plt.plot(N)                       #Construimos o gráfico
plt.show()                        #Plotamos

Euler implícito

A equação da reta obtida no euler explícito pode ser obtida a partir da definição da derivada:

Mas também podemos escrever a derivada como:

Mantendo a notação:

O termo não é conhecido, por isso temos uma equação implícita para . Métodos implícitos podem ser usados quando temos restrições muito rigorosas no método explícito devido a condições de estabilidade.

Exemplo: Decaimento

Trabalhando novamente com o decaimento radioativo:

Vamo ter então:

Ou mais explicitamente:

Uma comparação entre os dois métodos de Euler para o caso do decaimento é simples. Lembando da fórmula recursiva de ambos os casos:

E fazendo uma expansão em série de Taylor em torno de , escrevendo , temos:

Então os termos são iguais até a primeira ordem, sendo assim uma boa aproximação.

Implementando o método de euler implícito, temos:

import matplotlib.pyplot as plt   #Biblioteca para plotar gráficos
import numpy as np                #Biblioteca com funções matemáticas
N =[10**6];nt=6;lam=0.1;dt=0.1    #Parâmetros

tau=np.log(2)/lam                 #Meia-vida
fac = 1+lam*dt                    #Função calculada
for i in range(int(nt*tau/dt)+1): #Vamo calcular um múltiplo da meia vida
  N.append(N[i]/fac)              #Salvamos o novo valor
  print(i*dt,N[i])                #printamos o resultado

plt.plot(N)                       #Construimos o gráfico
plt.show()                        #Plotamos

Avaliação de erro

Para o cálculo de erros deste exemplo, utilizamos bastante a meia-vida como um valor caratectístico do sistema, principalmente para nos basearmos até definir até onde a simulação vai rodar. A meia vida é simplesmente o tempo necessário para a quantidade de núcleos cair para a metade do valor inicial. Isto é, primeiro temos que a solução analítica é simplesmente:

Então a meia vida é quando :

Erros computados

Erro 1: curva numérica vs exata

import numpy as np                #Biblioteca com recursos matemáticos
import matplotlib.pyplot as plt   #Biblioteca para plotar gráficos
N0=10**6;nt=6;lam=0.001;dt=100    #Parâmetros

exp = 1-lam*dt                    #Função calculada para o método explícito
imp = 1+lam*dt                    #Função calculada para o método implícito
NE  = [N0]                        #Lista de valore para o método explícito
NI  = [N0]                        #Lista de valore para o método implícito
T   = [0]                         #Lista de tempos no qual vamos calcular N

tau=np.log(2)/lam                 #Tempo final em que vamos obter o erro

#Solução numérica
for i in range(int(nt*tau/dt)+1): #Vamo calcular um múltiplo da meia vida
  NE.append(exp*NE[i])            #Explícito
  NI.append(NI[i]/imp)            #Implícito
  T.append((i+1)*dt)              #Tempo

#Solução analítica
t = np.arange(0.0,nt*tau, 1)      #Tempo
NA =N0*np.exp(-lam*t)             #Solução

#Construimos o gráfico
plt.plot(T,NE,color='blue', label='Explícito')                             #Solução explícita
plt.plot(T,NI,color='orange', label='Implícito')                           #Solução implícita    
plt.plot(t,NA,color='green', label='Sol. Exata')                           #Solução exata           
plt.title('Decaimento radioativo para  λ = '+str(lam)+' e dt = '+str(dt))  #Título
plt.legend()                                                               #Legenda das curvas
plt.xlabel('Tempo [s]');plt.ylabel('Número de Núcleos [N]')                #Legenda dos eixos
plt.ticklabel_format(style='plain')                                        #Desativar a notação científica
plt.show()                                                                 #Exibir o resultado

Erro 2: Curva de erro pra diferentes dt

import numpy as np                #Biblioteca com recursos matemáticos
import matplotlib.pyplot as plt   #Biblioteca para plotar gráficos
N0 =10**6;lam=0.001;nt=6          #Parâmetros
dt1=100;dt2=200;dt3=50

def sol_metodo(dt,N0,lam,nt):
  """Função para calcular a solução de acordo com o método de Euler e
  retornar o erro comparado a solução analítica e o tempo equivalente"""
  #dt   - Intervalo do passo
  #N0   - Número inicial de Núcleos
  #lam  - Lambda
  #nt  - Quantidade de meia vidas que vamos caclular
  N=[N0]            #Solução
  D=[0]             #Erro
  T=[0]             #Tempo
  frac = (1-lam*dt) #Funçao calculada
  tau=np.log(2)/lam                    #Meia-vida
  for i in range(1,int(nt*tau/dt)+1):  #Vamo calcular um múltiplo da meia vida
    N.append(frac*N[i-1])              #Solução numérica
    T.append(i*dt)                     #Tempo no qual achamos a solução
    err=abs(N[i]-N0*np.exp(-lam*T[i])) #Erro
    D.append(err)                      #Guardamos o erro
  return(D,T)                          #Retorno o erro

(D1,T1) = sol_metodo(dt1,N0,lam,nt)    #Solução para dt1
(D2,T2) = sol_metodo(dt2,N0,lam,nt)    #Solução para dt2
(D3,T3) = sol_metodo(dt3,N0,lam,nt)    #Solução para dt3

#Construimos o gráfico
plt.plot(T1,D1,color='blue', label='dt = '+str(dt1)+' s')              #Solução para dt1
plt.plot(T2,D2,color='orange', label='dt = '+str(dt2)+' s')            #Solução para dt2
plt.plot(T3,D3,color='green', label='dt = '+str(dt3)+' s')             #Solução para dt3
plt.title('Comparação dos erros para dierentes dt (Euler explícito)')  #Título
plt.legend()                                                           #Legenda das curvas
plt.xlabel('Tempo [s]');plt.ylabel('Erro[abs(N-Nmétodo)]')             #Legenda dos eixos
plt.ticklabel_format(style='plain')                                    #Desativar a notação científica
plt.show()                                  

Erro 3: Erro em relação a dt com um t fixo

import numpy as np                         #Biblioteca com recursos matemáticos
import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos
N0 =10**6;lam=0.001                        #Parâmetros
dt=[1/10**(i) for i in np.arange(2.9,-0.1,-0.3)]
tfin=[(1,'ob','-k'),(4,'sr','-m')]         #Tuplas com os tempos finais e os marcadores principais e do fit
y="log";x="log"                            #Se os eixos vão ser "linear" ou "log" 

def sol_metodo(dtT,N0,lam,t):
  """Função para calcular a solução de acordo com o método de Euler e
  retornar o erro comparado a solução analítica"""
  #dtT  - Intervalo do passo em unidade de tau
  #N0   - Número inicial de Núcleos
  #lam  - Lambda
  #t    - Tempo em que vamos pegar a solução, múltiplo da meia-vida
  N=N0
  T=np.log(2)/lam                  #Tempo final em que vamos obter o erro
  dt=dtT*T                         #Intervalo de tempo em segundos
  frac = (1-lam*dt)                #Funçao calculada
  for i in range(1,int(t*T/dt)+1): #Vamo calcular os passos necessários para 4000s
    N=frac*N                       #Solução numérica
  err =abs(N-N0*np.exp(-lam*i*dt)) #Erro
  return(err)                      #Retorno o erro

for t in tfin:                                  #Percorremos todos os tempos finais
  E=[]                                          #Lista para guardar todos os erros com
  for d in dt:
    E.append (sol_metodo(d,N0,lam,t[0]))               #Computamos o erro para cada dt
  plt.plot(dt,E,t[1], label='t = '+str(t[0])+' τ')     #Plot dos marcadores
  plt.plot(dt,E,'--'+t[1][1])                          #Plot Tracejado
  m, b = np.polyfit(np.log(np.array(dt)),np.log(np.array(E)), 1) #Uma curva para fitar os logaritmos
  label='ln(erro) = {:.2f} + {:.2f} ln(Δt)'.format(b,m)
  plt.plot(dt, np.e**b*np.array(dt)**m,t[2],label=label)         #Tiramos a exponencial do logaritmo

#Configuramos o gráfico
plt.title('Erro X dt')                               #Título
plt.legend()                                         #Legenda das curvas
plt.xlabel('Δt  [τ]')                                #Legenda dos eixos
plt.ylabel('Erro[abs(N-Nmétodo)]')  
plt.ticklabel_format(style='plain')                  #Desativar a notação científica
plt.xscale(x);plt.yscale(y)                          #Escala
plt.show()                                           #Exibir o resultado

Erro teórico

Outro caminho que podemos tomar para obtermos o método de Euler, é através da expansão em série de Taylor de uma função em torno do ponto . Lembrando que a expansão em série de Taylor de uma função em torno de é dada por:

Então fazendo esta expansão e escrevendo

Podemos olhar com mais atenção uma vez que tecnicamente nosssa variável da função que estamos expandindo é .

Considerando que é constante, então: , e sendo o ponto , temos:

De forma que ficamos então apenas com:

Os primeiros dois termos então correspondem ao método de Euler, e o somatório se torna o erro ao jogarmos fora. O menor grau que desprezamos é , então o erro local é . E como o número de passos é dado por então após passos temos uma estimativa de erro global dada por ou apenas . Cálculos mais detalhados da estimativa do erro podem ser encontrados em textos matemáticos[1].

Exemplo: SIR

O modelo SIR é um modelo de epidemia compartimental no qual separamos a população em três compartimentos:

  • : População suscetível
  • : População infectada
  • : População de recuperados

O modelo SIR é descrito então pelo seguinte conjunto de equações:

Onde temos as seguintes variáveis:

  • : transmissibilidade
  • : taxa de recuperação
    • É o inverso do tempo que a pessoa permanece infectada .

Neste modelo, as pessoas podem se transferir entre os compartimento apenas em um sentido , e pode-se destacar que comporta tanto pessoas que se curaram e tornaram-se imunes, quanto falecimentos, em todo, caso, são pessoas que não podem ser suscetíveis a infecção novamente. Como a população total é constante, temos e por conservação:

Desta forma só precisamos resolver duas equações a cada passo. Além disso podemos considerar que a transmissibilidade de uma doença depende da população total, de forma que , sendo um parâmetro próprio da doença e independente da população. Substituindo e reescrevendo cada compartimento como uma fração da população total, ao invés de números totais, isto é , seja um compartimento qualquer, então podemos reescrever o sistema como:

Um ponto interessante é a condição de epidemia. A segunda equação nos dá a variação de infectados, para termos epidemia precisamos que este valor positivo, isto é , logo:

Considerando como condição inicial que , isto é, praticamente toda população é suscetível, logo e temos então o seguinte limiar para uma condição inicial de epidêma. Podemos definir assim o número reprodutivo básico:

Além disso, um modelo SIRS (no qual recuperados podem tornar-se suscetíveis novamente) e com atraso foi discutido aqui

import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos

b   = 0.1  #Transmissibilidade
tau = 11   #Tempo que a pessoa permanece infectada
g   = 1/tau#Taxa de recuperação
N   = 1000  #População total
Np  = 5000  #Número de passos
dt  = 0.5  #Intervalo de tempo
#Valores iniciais
i=[1/N]; s=[1-i[0]]; r=[0]; t=[0]
T=[s[0]+i[0]+r[0]]  #Total

for it  in range(Np):
  s.append(s[it]-dt*b*s[it]*i[it])
  i.append(i[it]+dt*(b*s[it]*i[it]-g*i[it]))
  r.append(1-(i[it+1]+s[it+1]))
  t.append(dt+it*dt)
  T.append(s[it+1]+i[it+1]+r[it+1])

plt.plot(t,i) #Infectados 
plt.plot(t,s) #Suscetíveis
plt.plot(t,r) #Recuperados
plt.plot(t,T) #Total

Exemplo: Sistema massa-mola

Podemos escrever o problema geral da mecânica clássica como:

O método de Euler é um método para equações de 1ª ordem, então transformando em um sistema de equações de primeira ordem, escrevendo , temos:

E pela segunda lei de Newton , então ficamos simplesmente com:

E aplicando o método de Euler:

Agora tendo a lei de Hooke:

E reescrevendo , temo que:

Onde . Então o método de Euler fica:

A solução analítica do sistema é dada por:

Onde podemos ver que a energia e conserva:

Sendo , então:


import matplotlib.pyplot as plt            #Biblioteca para plotar gráficos

Np  = 200000              #Número de passos
dt  = 0.001               #Intervalo de tempo
m=1  ; k= 1.; w2= k/m     #Constantes
x=[1]; v=[0]; t=[0]; E=[k*(x[0]**2)/2+m*(v[0]**2)/2] #Valores iniciais

#Método de Euler
for it  in range(Np):
  x.append(x[it]+dt*v[it])
  v.append(v[it]-dt*x[it]*w2)
  E.append(k*x[it+1]**2/2+m*v[it+1]**2/2)
  t.append(dt+it*dt)

plt.plot(t,x) #Posição 
plt.plot(t,v) #Velocidade
plt.plot(t,E) #Energia

Exemplo: Modelo de Lotka Volterra

Além disso temos uma discussão sobre o Modelo de Lotka-Volterra e também sobre Modelo de Lotka-Volterra amortecido, que além da discussão apresenta uma solução numérica utilizando o Método de Euler, com a única diferença de um termo de amortecimento. Uma vez que zerarmos este termo, retornamos ao Modelo de Lotka-Volterra clássico.


Citações

  1. Euler's Method (William F. Trench, LibreTexts)

Principais materiais utilizados

  1. Erro de truncamento (REAMAT, UFRGS)
  2. Euler's Method (Paul Dawkins, Universidade Lamar)
  3. Forward and Backward Euler Methods (Michael Zeltkevic, Instituto de Tecnologia de Massachusetts)
  4. Método de Euler (REAMAT, UFRGS)