积分 (scipy.integrate)#
scipy.integrate 子包提供了几种积分技术，包括一个常微分方程积分器。模块的概述可以通过帮助命令获得：
>>> help(integrate)
Methods for Integrating Functions given function object.
quad -- General purpose integration.
dblquad -- General purpose double integration.
tplquad -- General purpose triple integration.
fixed_quad -- Integrate func(x) using Gaussian quadrature of order n.
quadrature -- Integrate with given tolerance using Gaussian quadrature.
romberg -- Integrate func using Romberg integration.
Methods for Integrating Functions given fixed samples.
trapezoid -- Use trapezoidal rule to compute integral.
cumulative_trapezoid -- Use trapezoidal rule to cumulatively compute integral.
simpson -- Use Simpson's rule to compute integral from samples.
romb -- Use Romberg Integration to compute integral from
-- (2**k + 1) evenly-spaced samples.
See the special module's orthogonal polynomials (special) for Gaussian
quadrature roots and weights for other weighting factors and regions.
Interface to numerical integrators of ODE systems.
odeint -- General integration of ordinary differential equations.
ode -- Integrate ODE using VODE and ZVODE routines.
一般积分 (quad)#
函数 quad 用于在两个点之间对一个变量的函数进行积分。这些点可以是 \(\pm\infty\) (\(\pm\) inf) 以表示无限极限。例如，假设你希望在区间 \([0, 4.5]\) 上对贝塞尔函数 jv(2.5, x) 进行积分。
这可以通过使用 quad 来计算：
>>> import scipy.integrate as integrate
>>> import scipy.special as special
>>> result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5)
>>> result
(1.1178179380783249, 7.8663172481899801e-09)
>>> from numpy import sqrt, sin, cos, pi
>>> I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5) - 4.0/27*sqrt(2)*sin(4.5) +
... sqrt(2*pi) * special.fresnel(3/sqrt(pi))[0])
>>> I
1.117817938088701
>>> print(abs(result[0]-I))
1.03761443881e-11
quad 的第一个参数是一个“可调用”的 Python 对象（即，函数、方法或类实例）。注意在这个例子中使用了 lambda 函数作为参数。接下来的两个参数是积分的极限。返回值是一个元组，第一个元素是积分的估计值，第二个元素是积分误差的估计值。注意，在这种情况下，这个积分的真实值是
其中
是菲涅尔正弦积分。注意，数值计算的积分结果与精确结果的误差在 \(1.04\times10^{-11}\) 以内——远低于报告的误差估计。
如果被积函数需要额外的参数，可以在 args 参数中提供。假设需要计算以下积分：
可以通过以下代码进行积分计算：
>>> from scipy.integrate import quad
>>> def integrand(x, a, b):
... return a*x**2 + b
...
>>> a = 2
>>> b = 1
>>> I = quad(integrand, 0, 1, args=(a,b))
>>> I
(1.6666666666666667, 1.8503717077085944e-14)
quad 函数也允许输入为无穷大，通过使用 \(\pm\) inf 作为参数之一。例如，假设需要数值计算指数积分：
（并且忘记了该积分可以通过 special.expn(n,x) 计算）。可以通过基于 quad 的例程定义一个新函数 vec_expint 来复制 special.expn 的功能：
>>> from scipy.integrate import quad
>>> import numpy as np
>>> def integrand(t, n, x):
... return np.exp(-x*t) / t**n
...
>>> def expint(n, x):
... return quad(integrand, 1, np.inf, args=(n, x))[0]
...
>>> vec_expint = np.vectorize(expint)
>>> vec_expint(3, np.arange(1.0, 4.0, 0.5))
array([ 0.1097, 0.0567, 0.0301, 0.0163, 0.0089, 0.0049])
>>> import scipy.special as special
>>> special.expn(3, np.arange(1.0,4.0,0.5))
array([ 0.1097, 0.0567, 0.0301, 0.0163, 0.0089, 0.0049])
被积函数甚至可以使用 quad 参数（尽管由于使用 quad 可能导致的数值误差，误差界限可能会低估误差）。在这种情况下，积分是
>>> result = quad(lambda x: expint(3, x), 0, np.inf)
>>> print(result)
(0.33333333324560266, 2.8548934485373678e-09)
>>> I3 = 1.0/3.0
>>> print(I3)
0.333333333333
>>> print(I3 - result[0])
8.77306560731e-11
最后一个例子表明，可以通过重复调用 quad 来处理多重积分。
警告
数值积分算法在有限数量的点上对被积函数进行采样。因此，它们无法保证任意被积函数和积分限的准确结果（或准确度估计）。例如，考虑高斯积分：
>>> def gaussian(x):
... return np.exp(-x**2)
>>> res = integrate.quad(gaussian, -np.inf, np.inf)
>>> res
(1.7724538509055159, 1.4202636756659625e-08)
>>> np.allclose(res[0], np.sqrt(np.pi)) # 与理论结果比较
True
由于被积函数在原点附近几乎为零，我们期望有限但较大的积分限也能得到相同的结果。然而：
>>> integrate.quad(gaussian, -10000, 10000)
(1.975190562208035e-203, 0.0)
这是因为 quad 中实现的适应性积分程序虽然按设计工作，但在如此大的有限区间内没有注意到函数的小而重要的部分。为了获得最佳结果，请考虑使用紧密包围被积函数重要部分的积分限。
>>> integrate.quad(gaussian, -15, 15)
(1.772453850905516, 8.476526631214648e-11)
具有多个重要区域的被积函数可以根据需要分解为若干部分。
一般多重积分（dblquad、tplquad、nquad）#
双重积分和三重积分的机制已经被封装在函数 dblquad 和 tplquad 中。这些函数接受要积分的函数以及四个或六个参数，分别对应不同的积分限。所有内层积分的限需要定义为函数。
以下是使用双重积分计算多个 \(I_{n}\) 值的示例：
>>> from scipy.integrate import quad, dblquad
>>> def I(n):
... return dblquad(lambda t, x: np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf)
...
>>> print(I(4))
(0.2500000000043577, 1.29830334693681e-08)
>>> print(I(3))
(0.33333333325010883, 1.3888461883425516e-08)
>>> print(I(2))
(0.4999999999985751, 1.3894083651858995e-08)
作为非恒定积分限的示例，考虑以下积分：
该积分可以通过以下表达式进行计算（注意内层积分上限使用了非恒定的 lambda 函数）：
>>> from scipy.integrate import dblquad
>>> area = dblquad(lambda x, y: x*y, 0, 0.5, lambda x: 0, lambda x: 1-2*x)
>>> area
(0.010416666666666668, 1.1564823173178715e-16)
对于 n 重积分，scipy 提供了函数 nquad。积分限是一个可迭代对象：可以是常数限的列表，也可以是非恒定积分限的函数列表。积分的顺序（以及因此的积分限）是从最内层积分到最外层积分。
上述积分
可以计算为
>>> from scipy import integrate
>>> N = 5
>>> def f(t, x):
... return np.exp(-x*t) / t**N
...
>>> integrate.nquad(f, [[1, np.inf],[0, np.inf]])
(0.20000000000002294, 1.2239614263187945e-08)
注意，f 的参数顺序必须与积分区间的顺序匹配；即，关于 \(t\) 的内积分在区间 \([1, \infty]\) 上，而关于 \(x\) 的外积分在区间 \([0, \infty]\) 上。
非恒定的积分区间可以以类似的方式处理；上述示例
可以通过以下方式进行评估
>>> from scipy import integrate
>>> def f(x, y):
... return x*y
...
>>> def bounds_y():
... return [0, 0.5]
...
>>> def bounds_x(y):
... return [0, 1-2*y]
...
>>> integrate.nquad(f, [bounds_x, bounds_y])
(0.010416666666666668, 4.101620128472366e-16)
这与之前的结果相同。
高斯求积#
fixed_quad 在固定区间上执行固定阶数的高斯求积。此函数使用 scipy.special 提供的正交多项式集合，可以计算大量正交多项式的根和求积权重（多项式本身作为返回多项式类实例的特殊函数可用 — 例如，special.legendre）。
使用样本进行积分#
如果样本是等间隔的，并且可用样本的数量为 \(2^{k}+1`（对于某个整数 :math:`k\)），则可以使用 Romberg romb 积分来获得高精度的积分估计值。Romberg 积分使用与 2 的幂次相关的步长进行梯形法则，然后对这些估计值执行 Richardson 外推法，以更高精度逼近积分。
对于任意间隔的样本，两个函数 trapezoid 和 simpson 是梯形法则和辛普森法则的直接应用。
和 simpson 是可用的。它们分别使用一阶和二阶的 Newton-Coates 公式来进行积分。梯形法则将函数近似为相邻点之间的直线，而辛普森法则将函数在三个相邻点之间近似为抛物线。
对于奇数个等间距的样本，如果函数是三次或更低次的多项式，辛普森法则可以精确积分。如果样本不等间距，则结果仅在函数是二次或更低次的多项式时精确。
>>> import numpy as np
>>> def f1(x):
... return x**2
...
>>> def f2(x):
... return x**3
...
>>> x = np.array([1,3,4])
>>> y1 = f1(x)
>>> from scipy import integrate
>>> I1 = integrate.simpson(y1, x=x)
>>> print(I1)
21.0
这完全对应于
而积分第二个函数
>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x=x)
>>> print(I2)
61.5
并不对应于
因为 f2 中的多项式阶数大于二。
使用低级回调函数进行更快的积分#
希望减少积分时间的用户可以通过 scipy.LowLevelCallable 传递一个 C 函数指针给 quad、dblquad、tplquad 或 nquad，它将在 Python 中进行积分并返回结果。这里的性能提升来自两个因素。主要改进是更快的函数求值，这是由函数本身的编译提供的。此外，我们在 quad 中移除了 C 和 Python 之间的函数调用，从而提供了加速。对于像正弦这样的简单函数，这种方法可能会提供约 2 倍的速度提升，但对于更复杂的函数，可以产生更显著的改进（10 倍以上）。因此，这个功能是为那些愿意进行大量数值积分的用户设计的。
编写少量C代码可以显著减少计算时间。
例如，可以通过`ctypes`在几个简单的步骤中使用这种方法：
1.) 用函数签名``double f(int n, double *x, void *user_data)``编写一个C语言的被积函数，其中``x``是一个数组，包含函数f被评估的点，而``user_data``用于提供你想要提供的任意额外数据。
/* testlib.c */
double f(int n, double *x, void *user_data) {
double c = *(double *)user_data;
return c + x[0] - x[1] * x[2]; /* 对应于 c + x - y * z */
}
2.) 现在将这个文件编译成共享/动态库（快速搜索将有助于此操作，因为它依赖于操作系统）。用户必须链接任何使用的数学库等。在Linux上，这看起来像:
$ gcc -shared -fPIC -o testlib.so testlib.c
输出库将被引用为``testlib.so``，但它可能有不同的文件扩展名。现在已经创建了一个库，可以使用`ctypes`加载到Python中。
3.) 使用`ctypes`将共享库加载到Python中，并设置``restypes``和``argtypes`` - 这允许SciPy正确解释函数：
import os, ctypes
from scipy import integrate, LowLevelCallable
lib = ctypes.CDLL(os.path.abspath('testlib.so'))
lib.f.restype = ctypes.c_double
lib.f.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)
c = ctypes.c_double(1.0)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)
func = LowLevelCallable(lib.f, user_data)
函数中的最后一个``void *user_data``是可选的，如果不需要，可以在C函数和ctypes argtypes中省略。请注意，坐标作为双精度数组传递，而不是作为单独的参数。
4.) 现在像往常一样对库函数进行积分，这里使用`nquad`：
>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11)
Python 元组在减少的时间内按预期返回。所有可选参数都可以与该方法一起使用，包括指定奇点、无限边界等。
常微分方程 (solve_ivp)#
在给定初始条件的情况下，积分一组常微分方程 (ODE) 是另一个有用的例子。SciPy 中的函数 solve_ivp 可用于积分一阶向量微分方程：
给定初始条件 \(\mathbf{y}\left(0\right)=y_{0}\)，其中 \(\mathbf{y}\) 是一个长度为 \(N\) 的向量，而 \(\mathbf{f}\) 是从 \(\mathcal{R}^{N}\) 到 \(\mathcal{R}^{N}\) 的映射。通过将中间导数引入 \(\mathbf{y}\) 向量，总是可以将高阶常微分方程简化为这种类型的微分方程。
例如，假设希望求解以下二阶微分方程的解：
初始条件为 \(w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}\) 和 \(\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}\)。已知该微分方程在这些边界条件下的解是 Airy 函数
这提供了一种使用 special.airy 检查积分器的方法。
首先，通过设置 \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) 和 \(t=z\) 将此 ODE 转换为标准形式。因此，微分方程变为
换句话说，
作为一个有趣的提醒，如果 \(\mathbf{A}\left(t\right)\) 与 \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) 在矩阵乘法下可交换，那么这个线性微分方程 有一个精确解，使用矩阵指数：
然而，在这种情况下，\(\mathbf{A}\left(t\right)\) 和它的积分不可交换。
这个微分方程可以使用函数 solve_ivp 来求解。
它需要导数，fprime，时间跨度 [t_start, t_end]
以及初始条件向量，y0，作为输入参数，并返回
一个对象，其 y 字段是一个数组，连续的解值作为
列。因此，初始条件在第一列输出中给出。
>>> from scipy.integrate import solve_ivp
>>> from scipy.special import gamma, airy
>>> y1_0 = +1 / 3**(2/3) / gamma(2/3)
>>> y0_0 = -1 / 3**(1/3) / gamma(1/3)
>>> y0 = [y0_0, y1_0]
>>> def func(t, y):
... return [t*y[1],y[0]]
...
>>> t_span = [0, 4]
>>> sol1 = solve_ivp(func, t_span, y0)
>>> print("sol1.t: {}".format(sol1.t))
sol1.t: [0. 0.10097672 1.04643602 1.91060117 2.49872472 3.08684827
3.62692846 4. ]
如上所示，solve_ivp 如果没有指定其他时间步长，它会自动确定时间步长。
为了将 solve_ivp 的解与 airy 函数进行比较，solve_ivp 创建的时间向量被传递给 airy 函数。
>>> print("sol1.y[1]: {}".format(sol1.y[1]))
sol1.y[1]: [0.35502805 0.328952 0.12801343 0.04008508 0.01601291 0.00623879
0.00356316 0.00405982]
>>> print("airy(sol.t)[0]: {}".format(airy(sol1.t)[0]))
airy(sol.t)[0]: [0.35502805 0.328952 0.12804768 0.03995804 0.01575943 0.00562799
0.00201689 0.00095156]
使用 solve_ivp 的标准参数求解结果与 airy 函数存在较大偏差。为了最小化这种偏差，可以使用相对和绝对容差。
>>> rtol, atol = (1e-8, 1e-8)
>>> sol2 = solve_ivp(func, t_span, y0, rtol=rtol, atol=atol)
>>> print("sol2.y[1][::6]: {}".format(sol2.y[1][0::6]))
sol2.y[1][::6]: [0.35502805 0.19145234 0.06368989 0.0205917 0.00554734 0.00106409]
>>> print("airy(sol2.t)[0][::6]: {}".format(airy(sol2.t)[0][::6]))
airy(sol2.t)[0][::6]: [0.35502805 0.19145234 0.06368989 0.0205917 0.00554733 0.00106406]
为了指定 solve_ivp 求解的用户定义时间点，solve_ivp 提供了两种可以互补使用的可能性。通过将 t_eval 选项传递给函数调用，solve_ivp 会在其输出中返回这些 t_eval 时间点的解。
>>> import numpy as np
>>> t = np.linspace(0, 4, 100)
>>> sol3 = solve_ivp(func, t_span, y0, t_eval=t)
如果已知函数的雅可比矩阵，可以将其传递给 solve_ivp 以获得更好的结果。但请注意，默认的积分方法 RK45 不支持雅可比矩阵，因此必须选择另一种积分方法。支持雅可比矩阵的积分方法之一是例如以下示例中的 Radau 方法。
>>> def gradient(t, y):
... return [[0,t], [1,0]]
>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient)
求解具有带状雅可比矩阵的系统#
odeint 可以被告知雅可比矩阵是*带状的*。对于已知为刚性的大型微分方程系统，这可以显著提高性能。
作为一个例子，我们将使用线法[MOL]_来求解一维Gray-Scott偏微分方程。Gray-Scott方程在区间 \(x \in [0, L]\) 上描述了函数 \(u(x, t)\) 和 \(v(x, t)\) 的演化，形式如下：
其中 \(D_u\) 和 \(D_v\) 分别是组分 \(u\) 和 \(v\) 的扩散系数，而 \(f\) 和 \(k\) 是常数。（有关该系统的更多信息，请参见 http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/）
我们假设Neumann（即“无通量”）边界条件：
为了应用线法，我们通过定义均匀分布的 \(N\) 个点 \(\left\{x_0, x_1, \ldots, x_{N-1}\right\}\) 来离散化 \(x\) 变量，其中 \(x_0 = 0\) 且 \(x_{N-1} = L\)。我们定义 \(u_j(t) \equiv u(x_k, t)\) 和 \(v_j(t) \equiv v(x_k, t)\)，并用有限差分替换 \(x\) 的导数。即，
这样我们就得到了一个由 \(2N\) 个常微分方程组成的系统：
为方便起见，省略了 \((t)\) 参数。 为了强制边界条件，我们引入“幽灵”点 \(x_{-1}\) 和 \(x_N\)，并定义 \(u_{-1}(t) \equiv u_1(t)\)，\(u_N(t) \equiv u_{N-2}(t)\)；\(v_{-1}(t)\) 和 \(v_N(t)\) 的定义类似。
然后
以及
我们完整的 \(2N\) 个常微分方程系统是 (1) 对于 \(k = 1, 2, \ldots, N-2\)，以及 (2) 和 (3)。
我们现在可以开始在代码中实现这个系统。我们必须将 \(\{u_k\}\) 和 \(\{v_k\}\) 合并成一个长度为 \(2N\) 的向量。两个明显的选择是 \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) 和 \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\)。从数学上讲，这无关紧要，但选择会影响 odeint 解决系统的效率。原因是顺序如何影响雅可比矩阵中非零元素的模式。
当变量按 \(\{u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}\}\) 排序时，雅可比矩阵中非零元素的模式是
雅可比模式的变量交错排列为 \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\) 时为
在这两种情况下，都只有五个非平凡的对角线，但当变量交错排列时，带宽要小得多。
即，主对角线及其上方和下方紧邻的两条对角线是非零对角线。
这一点很重要，因为 odeint 的输入参数 mu 和 ml 是雅可比矩阵的上带宽和下带宽。
当变量交错排列时，mu 和 ml 为 2。
当变量堆叠排列，且 \(\{v_k\}\) 跟随 \(\{u_k\}\) 时，上带宽和下带宽为 \(N\)。
做出这一决定后，我们可以编写实现微分方程组的函数。
首先，我们定义系统源项和反应项的函数:
def G(u, v, f, k):
return f * (1 - u) - u*v**2
def H(u, v, f, k):
return -(f + k) * v + u*v**2
接下来，我们定义计算微分方程组右侧的函数:
def grayscott1d(y, t, f, k, Du, Dv, dx):
"""
一维 Gray-Scott 方程的微分方程。
使用线法推导 ODE。
"""
# 向量 u 和 v 在 y 中交错排列。我们通过切片 y 来定义 u 和 v 的视图。
u = y[::2]
v = y[1::2]
# dydt 是此函数的返回值。
dydt = np.empty_like(y)
# 就像 u 和 v 是 y 中交错向量的视图一样，dudt 和 dvdt 是 dydt 中交错输出向量的视图。
dudt = dydt[::2]
dvdt = dydt[1::2]
# 计算 du/dt 和 dv/dt。端点和内部点分别处理。
dudt[0] = G(u[0], v[0], f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2
dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2
dudt[-1] = G(u[-1], v[-1], f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2
dvdt[0] = H(u[0], v[0], f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2
dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2
dvdt[-1] = H(u[-1], v[-1], f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2
return dydt
我们不会实现一个计算雅可比矩阵的函数，但我们会告诉`odeint`，雅可比矩阵是带状的。这使得底层的求解器（LSODA）可以避免计算它知道为零的值。对于一个大型系统，这显著提高了性能，如下面的ipython会话所示。
首先，我们定义所需的输入:
In [30]: rng = np.random.default_rng()
In [31]: y0 = rng.standard_normal(5000)
In [32]: t = np.linspace(0, 50, 11)
In [33]: f = 0.024
In [34]: k = 0.055
In [35]: Du = 0.01
In [36]: Dv = 0.005
In [37]: dx = 0.025
在不利用雅可比矩阵的带状结构的情况下计时计算:
In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
1 loop, best of 3: 25.2 s per loop
现在设置``ml=2``和``mu=2``，以便`odeint`知道雅可比矩阵是带状的:
In [39]: %timeit solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2)
10 loops, best of 3: 191 ms per loop
这快了很多！
让我们确保它们计算的结果是相同的:
In [41]: np.allclose(sola, solb)
Out[41]: True