Skip to content

Instantly share code, notes, and snippets.

@fccoelho
Created November 4, 2011 13:03
Show Gist options
  • Save fccoelho/1339271 to your computer and use it in GitHub Desktop.
Save fccoelho/1339271 to your computer and use it in GitHub Desktop.
Exemplos de integração numérica de sistemas de equações diferenciais ordinárias.
# coding:utf8
"""
Exemplos de integração numérica de sistemas de equações diferenciais ordinárias.
ode : interface para a coleção de solvers VODE http://www.netlib.org/ode
Flávio Codeço Coelho
Licença: GPL
"""
from scipy. integrate import ode
import sympy #opcional
import pylab as P #para plotar resultados
import numpy as np #
def fun1_vode(t,y): #
"""
Primeiro exemplo: crescimento exponencial
dy/dt = ay
VODE requer função que define as odes
receba arqumentos na ondem inversa de odeint
"""
a = .5
return a*y
def fun2_vode(t,y):
"""
Modelo Logístico
Crescimento populacional dependente da densidade.
EDO não-linear de primeira order
dy/dt = a(1-y/k)y
"""
a = .5 #taxa de crescimento
k = 1000.0 # capacidade e suporte
return a*(1-y/k)*y
def fun3_vode(t,y): #
"""
Primeiro exemplo: crescimento exponencial
dy/dt = ay
VODE requer função que define as odes
receba arqumentos na ondem inversa de odeint
"""
a = .5
return a*np.sin(.2*t)*y
def solve_vode(r,tf,dt):
"""
Resolve o sistema de equações
chamando a integração passo a passo e agregando o resultado em uma matriz
"""
res = np.zeros(10000)
i = 0
while r.successful() and r.t < tf:
r.integrate(r.t+dt)
res[i] = r.y
i += 1
return res
t = np.arange(0,100,.01)
y0 = 1e-6
#========= Resolvendo com VODE ==============
#~ r = ode(fun2_vode).set_integrator('vode',method='bdf', with_jacobian=False)
r = ode(fun2_vode).set_integrator('dopri5')
#~ r = ode(fun2_vode).set_integrator('dop835')
r.set_initial_value(y0,t[0])
y = solve_vode(r,100,.01)
#Visualizando a solução
P.plot(t,y)
P.grid()
P.ylabel('$y(t)$')
P.xlabel('$t$')
P.show()
# coding:utf8
"""
Exemplos de integração numérica de sistemas de equações diferenciais ordinárias.
Alguns exemplos são apresentados e resolvidos usando solver
oferecido pelo scipy:
odeint : interface para o solver LSODA
Flávio Codeço Coelho
Licença: GPL
"""
from scipy. integrate import odeint
import sympy #opcional
import pylab as P #para plotar resultados
import numpy as np
def fun1(y,t):
"""
Primeiro exemplo: crescimento exponencial
dy/dt = ay
"""
a = .5
return a*y[0]
def fun2(y,t):
"""
Modelo Logístico
Crescimento populacional dependente da densidade.
EDO não-linear de primeira order
dy/dt = a(1-y/k)y
"""
a = .5 #taxa de crescimento
k = 1000.0 # capacidade e suporte
return a*(1-y[0]/k)*y[0]
def fun3(y,t):
"""
Equação não-autônoma com crescimento modulado por uma senoide
dy/dt = ay
"""
a = .5
return a*np.sin(.2*t)*y[0]
t = np.arange(0,100,.01)
y0 = 1e-6
#========= Resolvendo com odeint ============
y = odeint(fun2,y0,t)
#Visualizando a solução
P.plot(t,y)
P.grid()
P.ylabel('$y(t)$')
P.xlabel('$t$')
P.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment