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
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
N =[10**6];Np=100;lam=0.1;dt=0.1 #Parâmetros
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
fac = 1-lam*dt                    #Função calculada
for i in range(Np):               #Vamo calcular Np passos
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
   N.append(fac*N[i])              #Salvamos o novo valor
   print(i*dt,N[i])                #printamos o resultado
   print(i*dt,N[i])                #printamos o resultado
Linha 123: Linha 124:


<pre>
<pre>
import matplotlib.pyplot as plt #Biblioteca para plotar gráficos
import matplotlib.pyplot as plt   #Biblioteca para plotar gráficos
N =[10**6];Np=100;lam=0.1;dt=0.1 #Parâmetros
N =[10**6];nt=6;lam=0.1;dt=0.1   #Parâmetros


fac = 1+lam*dt                   #Função calculada
tau=np.log(2)/lam                #Meia-vida
for i in range(Np):             #Vamo calcular Np passos
fac = 1+lam*dt                   #Função calculada
   N.append(N[i]/fac)             #Salvamos o novo valor
for i in range(int(nt*tau/dt)+1): #Vamo calcular um múltiplo da meia vida
   print(i*dt,N[i])               #printamos o resultado
   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.plot(N)                       #Construimos o gráfico
plt.show()                       #Plotamos
plt.show()                       #Plotamos
</pre>
</pre>


Linha 141: Linha 143:


<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
N0 =10**6;Np=40;lam=0.001;dt=100 #Parâmetros
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
exp = 1-lam*dt                    #Função calculada para o método explícito
Linha 150: Linha 153:
NI  = [N0]                        #Lista de valore para o método implícito
NI  = [N0]                        #Lista de valore para o método implícito
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


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


#Solução analítica
#Solução analítica
t = np.arange(0.0, 4000, 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 169: Linha 174:
plt.xlabel('Tempo [s]');plt.ylabel('Número de Núcleos [N]')                #Legenda dos eixos
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.ticklabel_format(style='plain')                                        #Desativar a notação científica
plt.show()                                                                #Exibir o resultado
plt.show()                                                                #Exibir o resultado                                                              #Exibir o resultado
</pre>
</pre>


Linha 180: Linha 185:
dt1=100;dt2=200;dt3=50
dt1=100;dt2=200;dt3=50


def sol_metodo(dt,N0,lam):
def sol_metodo(dt,N0,lam,nt):
   """Função para calcular a solução de acordo com o método de Euler e
   """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"""
   retornar o erro comparado a solução analítica e o tempo equivalente"""
Linha 186: Linha 191:
   #N0  - Número inicial de Núcleos
   #N0  - Número inicial de Núcleos
   #lam  - Lambda
   #lam  - Lambda
  #nt  - Quantidade de meia vidas que vamos caclular
   N=[N0]            #Solução
   N=[N0]            #Solução
   D=[0]            #Erro
   D=[0]            #Erro
   T=[0]            #Tempo
   T=[0]            #Tempo
   frac = (1-lam*dt) #Funçao calculada
   frac = (1-lam*dt) #Funçao calculada
   for i in range(1,int(4000/dt)+1):   #Vamo calcular os passos necessários para 4000s
  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
     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
Linha 197: Linha 204:
   return(D,T)                          #Retorno o erro
   return(D,T)                          #Retorno o erro


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


#Construimos o gráfico
#Construimos o gráfico
Linha 209: Linha 216:
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()                                                             #Exibir o resultado
plt.show()            
</pre>
</pre>
== Erro teórico ==
== Erro teórico ==



Edição das 14h10min de 31 de janeiro 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

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

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

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                                                               #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               #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,6)       #Solução para dt1
(D2,T2) = sol_metodo(dt2,N0,lam,6)       #Solução para dt2
(D3,T3) = sol_metodo(dt3,N0,lam,6)       #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 teórico

Principais materiais utilizados

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