A rotina de integração mais genérica é scipy.integrate.quad():
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
from scipy.integrate import quad
res, err = quad(np.sin, 0, np.pi/2)
np.allclose(res, 1)
np.allclose(err, 1 - res)[/pastacode]
Outros esquemas de integração estão disponíveis com fixed_quad, quadratura, romberg.
scipy.integrate também possui rotinas para integração de equações diferenciais ordinárias (ODE). Em particular, scipy.integrate.odeint() é um integrador de uso geral usando LSODA (Livermore Solver para Equações Diferenciais Ordinárias com o método automático de comutação para problemas rígidos e não-rígidos), consulte a biblioteca Fortran ODEPACK para mais detalhes.
odeint resolve sistemas de ODE de primeira ordem da forma:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
dy/dt = rhs(y1, y2, .., t0,...)[/pastacode]
Como uma introdução, vamos resolver a ODE dy / dt =-2a entre t = 0 .. 4, com a condição inicial y (t = 0) = 1. Primeiro, a função que calcula a derivada da posição precisa de ser definida:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
def calc_derivative(ypos, time, counter_arr):
counter_arr += 1
return -2 * ypos[/pastacode]
Um argumento extra counter_arr foi adicionado para ilustrar que a função pode ser chamada várias vezes por um único intervalo de tempo, até a convergência da solução. O array contador é definido como:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
counter = np.zeros((1,), dtype=np.uint16)[/pastacode]
A trajetória agora é calculada:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
from scipy.integrate import odeint
time_vec = np.linspace(0, 4, 40)
yvec, info = odeint(calc_derivative, 1, time_vec,
args=(counter,), full_output=True)[/pastacode]
Assim, a função derivada foi chamada de mais de 40 vezes (o número de intervalos de tempo):
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
counter[/pastacode]
e o número total de iterações para cada um dos 10 passos primeira vez pode ser obtido por:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
info['nfe'][:10][/pastacode]
Note que a resolução requer mais iterações para o primeiro intervalo de tempo. O solução yvec para a trajetória agora pode ser graficada:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
import numpy as np
from scipy.integrate import odeint
import pylab as pl
def calc_derivative(ypos, time):
return -2*ypos
time_vec = np.linspace(0, 4, 40)
yvec = odeint(calc_derivative, 1, time_vec)
pl.plot(time_vec, yvec)
pl.xlabel('Tempo [s]')
pl.ylabel(u'Posição y [m]')[/pastacode]
Outro exemplo com scipy.integrate.odeint() será um oscilador massa-mola amortecido (oscilador de 2ª ordem). A posição de uma massa presa a uma mola obedece à ODE de segunda ordem:[latexpage]
\begin{equation} \quicklatex{size=15} \dfrac{d^{2}y}{dx^{2}}+2.eps.w_{0}.\dfrac{dy}{dx}+w_{0}^{2}y=0 \end{equation}
com:
$\quicklatex{size=15} w_{0}^{2}=\dfrac{k}{m}$
sendo k a constante da mola e m a massa.
$\quicklatex{size=15} eps=\dfrac{c}{2 m w_{0}}$
onde c é o coeficiente de amortecimento.
Para este exemplo, nós escolhemos os parâmetros como:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
mass = 0.5 # kg
kspring = 4 # N/m
cviscous = 0.4 # N s/m[/pastacode]
de modo que o sistema será subamortecido, porque:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
eps = cviscous / (2 * mass * np.sqrt(kspring/mass))
eps < 1[/pastacode]
Para o scipy.integrate.odeint() a equação de 2ª ordem tem de ser transformada num sistema de duas equações de primeira ordem para o vetor Y = (y, y ‘). Será conveniente definir nu eps = 2 * wo = c / m e om wo = ^ 2 = k / m:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
nu_coef = cviscous / mass
om_coef = kspring / mass[/pastacode]
Assim, a função vai calcular a velocidade e aceleração por:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
def calc_deri(yvec, time, nuc, omc):
return (yvec[1], -nuc * yvec[1] - omc * yvec[0])
time_vec = np.linspace(0, 10, 100)
yarr = odeint(calc_deri, (1, 0), time_vec, args=(nu_coef, om_coef))[/pastacode]
A posição final e velocidade são mostrados na seguinte figura Matplotlib:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
import numpy as np
from scipy.integrate import odeint
import pylab as pl
mass = 0.5
kspring = 4
cviscous = 0.4
nu_coef = cviscous / mass
om_coef = kspring / mass
def calc_deri(yvec, time, nuc, omc):
return (yvec[1], -nuc * yvec[1] - omc * yvec[0])
time_vec = np.linspace(0, 10, 100)
yarr = odeint(calc_deri, (1, 0), time_vec, args=(nu_coef, om_coef))
pl.plot(time_vec, yarr[:, 0], label='y')
pl.plot(time_vec, yarr[:, 1], label="y'")
pl.legend()[/pastacode]
Não há solver de Equações Diferenciais Parciais (PDE) em Scipy. Alguns pacotes Python para a solução do PDE estão disponíveis, tais como fipy ou SfePy.