1.5.8 数值积分：

scipy.integrate.quad()是最常见的积分程序：

In [1]:

from scipy.integrate import quad
res, err = quad(np.sin, 0, np.pi/2)
np.allclose(res, 1)

Out[1]:

True

In [2]:

np.allclose(err, 1 - res)

Out[2]:

True

其他的积分程序可以在fixed_quad、 quadrature、romberg中找到。

scipy.integrate 可提供了常微分公式(ODE)的特色程序。特别的，scipy.integrate.odeint() 是使用LSODA（Livermore Solver for Ordinary Differential equations with Automatic method switching for stiff and non-stiff problems）的通用积分器，更多细节请见ODEPACK Fortran 库。

odeint解决如下形式的第一顺序ODE系统：

$dy/dt = rhs(y1, y2, .., t0,…)$

作为一个介绍，让我们解一下在初始条件下$y(t=0) = 1$，这个常微分公式$dy/dt = -2y$在$t = 0..4$时的值。首先，这个函数计算定义位置需要的导数：

In [3]:

def calc_derivative(ypos, time, counter_arr):
    counter_arr += 1
    return -2 * ypos

添加了一个额外的参数counter_arr用来说明这个函数可以在一个时间步骤被调用多次，直到收敛。计数器数组定义如下：

In [4]:

counter = np.zeros((1,), dtype=np.uint16)

现在计算轨迹线：

In [5]:

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)

因此，导数函数被调用了40多次（即时间步骤数）：

In [6]:

counter

Out[6]:

array([129], dtype=uint16)

前十个时间步骤的累积循环数，可以用如下方式获得：

In [7]:

info['nfe'][:10]

Out[7]:

array([31, 35, 43, 49, 53, 57, 59, 63, 65, 69], dtype=int32)

注意，求解器对于首个时间步骤需要更多的循环。导数答案yvec可以画出来：

阻尼弹簧重物振子（二阶振荡器）是使用scipy.integrate.odeint()的另一个例子。链接到弹簧的重物的位置服从二阶常微分方程$y'' + 2 eps wo y' + wo^2 y = 0$，其中$wo^2 = k/m$ 弹簧的常数为k, m是重物质量，$eps=c/(2 m wo)$，c是阻尼系数。例如，我们选择如下参数：

In [8]:

mass = 0.5  # kg
kspring = 4  # N/m
cviscous = 0.4  # N s/m

因此系统将是欠阻尼的，因为：

In [9]:

eps = cviscous / (2 * mass * np.sqrt(kspring/mass))
eps < 1

Out[9]:

True

对于scipy.integrate.odeint()求解器，二阶等式需要被变换为系统内向量$Y=(y, y')$的两个一阶等式。为了方便，定义$nu = 2 eps * wo = c / m$和$om = wo^2 = k/m$：

In [10]:

nu_coef = cviscous / mass
om_coef = kspring / mass

因此函数将计算速度和加速度：

In [11]:

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

如下的Matplotlib图片显示了最终的位置和速度：

在Sicpy中没有偏微分方程（PDE）求解器。存在其他求解PDE的Python包，比如fipy或SfePy。