Exemplos de integração numérica de sistemas de equações diferenciais ordinárias.

ode_tut.py

# 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()

odeint_tut.py

# 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()