集成 (scipy.integrate )¶
这个 scipy.integrate 子软件包提供了几种积分技术，包括常微分方程积分器。help命令提供了该模块的概述：
>>> 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 )¶
The function quad is provided to integrate a function of one
variable between two points. The points can be \(\pm\infty\)
(\(\pm\) inf) to indicate infinite limits. For example,
suppose you wish to integrate a bessel function jv(2.5, x) along
the interval \([0, 4.5].\)
这可以使用以下公式计算 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) 已被遗忘)。函数的功能 special.expn 可通过定义新函数进行复制 vec_expint 基于例行公事 quad ：
>>> from scipy.integrate import quad
>>> 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 )。这种情况下的积分是
>>> 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 。
广义多重积分 (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 它执行固定阶高斯正交。第二个函数是 quadrature 其执行多阶高斯求积，直到积分估计的差低于用户提供的某个容差。这两个函数都使用模块 scipy.special.orthogonal 其可以计算各种正交多项式的根和正交权重(多项式本身作为返回多项式类实例的特殊函数可用)-例如， special.legendre )。
Romberg集成¶
使用示例进行集成¶
如果样本是等间距的，并且可用的样本数量是 \(2^{{k}}+1\) 对于某个整数 \(k\) ，然后是龙伯格 romb 积分可以用来使用可用的样本获得积分的高精度估计。Romberg积分在步长与2的幂相关的步长上使用梯形规则，然后对这些估计进行Richardson外推，以更高的精度逼近积分。
在任意间距样本的情况下，这两个函数 trapezoid 和 simpson 都是有空的。他们分别使用一阶和二阶牛顿-科茨公式进行积分。梯形法则将函数近似为邻接点之间的直线，而辛普森法则将三个邻接点之间的函数近似为抛物线。
对于等间距的奇数个样本，如果函数是3阶或更低阶的多项式，辛普森法则是精确的。如果样本间距不相等，则只有当函数是2阶或更低阶多项式时，结果才是精确的。
>>> 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)
>>> print(I1)
21.0
这正好对应于
而集成第二函数
>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x)
>>> print(I2)
61.5
不对应于
因为f2中多项式的阶数大于2。
使用低级回调函数更快地集成¶
希望减少积分时间的用户可以将C函数指针传递给 scipy.LowLevelCallable 至 quad ， dblquad ， tplquad 或 nquad 它将被集成，并以Python格式返回结果。这里的性能提升源于两个因素。主要的改进是函数求值速度更快，这是通过编译函数本身提供的。此外，我们还通过删除C和Python之间的函数调用提供了加速 quad 。对于像正弦这样的普通函数，这种方法可以提供~2倍的速度改进，但对于更复杂的函数，可以产生更显著的改进(10倍以上)。因此，此功能面向具有数值密集型集成的用户，该用户愿意编写少量的C语言来显著减少计算时间。
该方法可以使用，例如，通过 ctypes 只需几个简单的步骤：
1.)用C语言编写带有函数签名的被积函数 double f(int n, double *x, void *user_data) ，在哪里 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]; /* corresponds to c + x - y * z */
}
2.)现在将该文件编译成一个共享/动态库(快速搜索会对此有所帮助，因为它依赖于操作系统)。用户必须链接所使用的任何数学库等。在Linux上，它看起来如下所示：：
$ gcc -shared -fPIC -o testlib.so testlib.c
输出库将称为 testlib.so ，但它可能具有不同的文件扩展名。现在已经创建了一个库，可以使用以下命令将其加载到Python中 ctypes 。
3.)使用将共享库加载到Python中 ctypes 并将其设置为 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函数和ctypeargtype中)。请注意，坐标作为双精度数组传入，而不是作为单独的参数传入。
4.)现在像往常一样集成库函数，这里使用 nquad ：
>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11)
Python元组按预期在较短的时间内返回。所有可选参数都可以与此方法一起使用，包括指定奇点、无限边界等。
常微分方程 (solve_ivp )¶
积分一组给定初始条件的常微分方程(ODE)是另一个有用的例子。该函数 solve_ivp 在SciPy中可用于积分一阶向量微分方程：
给定初始条件 \(\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)}}.\) 众所周知，这个微分方程在这些边界条件下的解就是艾里函数。
这给出了一种使用以下方法检查积分器的方法 special.airy 。
首先，通过设置将此颂歌转换为标准格式 \(\mathbf{{y}}=\left[\frac{{dw}}{{dz}},w\right]\) 和 \(t=z\) 。这样，微分方程就变成了
换句话说，
作为一个有趣的提醒，如果 \(\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 与其标准参数相比，与AIY函数有较大的偏差。为了使该偏差最小化，可以使用相对公差和绝对公差。
>>> 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 可以说雅可比是 带状 。对于已知刚性的大型微分方程系统，这可以显著提高性能。
作为例子，我们将使用直线法求解一维Gray-Scott偏微分方程 [MOL]. 函数的Gray-Scott方程 \(u(x, t)\) 和 \(v(x, t)\) 在区间上 \(x \in [0, L]\) 是
哪里 \(D_u\) 和 \(D_v\) 是各组分的扩散系数 \(u\) 和 \(v\) 、和 \(f\) 和 \(k\) 是常量。(有关该系统的更多信息，请参见http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/)
我们将假设Neumann(即“无通量”)边界条件：
为了应用直线方法，我们离散化 \(x\) 变量，通过定义 \(N\) 积分 \(\left\{{x_0, x_1, \ldots, x_{{N-1}}\right\}}\) ，具有 \(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}}\}}\) 是
在这两种情况下，只有五条非平凡的对角线，但是当变量交错时，带宽要小得多。也就是说，主对角线和紧靠其上的两条对角线和紧靠在主对角线下面的两条对角线是非零对角线。这一点很重要，因为输入 mu 和 ml 的 odeint 是雅可比矩阵的上带宽和下带宽。当变量交错时， 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):
"""
Differential equations for the 1-D Gray-Scott equations.
The ODEs are derived using the method of lines.
"""
# The vectors u and v are interleaved in y. We define
# views of u and v by slicing y.
u = y[::2]
v = y[1::2]
# dydt is the return value of this function.
dydt = np.empty_like(y)
# Just like u and v are views of the interleaved vectors
# in y, dudt and dvdt are views of the interleaved output
# vectors in dydt.
dudt = dydt[::2]
dvdt = dydt[1::2]
# Compute du/dt and dv/dt. The end points and the interior points
# are handled separately.
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