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]

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

yEy_linhaNã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.