整合（`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) 表示无限的范围。例如，假设您需要计算一个贝塞尔函数 jv(2.5, x) 在区间 \([0, 4.5].\) 上的积分。

\[I=\int_{0}^{4.5}J_{2.5}\left(x\right)\, dx.\]

这可以通过 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- 函数作为参数。接下来的两个参数是积分的上下限。返回结果是一个组，第一个元素保存了积分的估计值，第二个元素保存了绝对积分误差的估计值。请注意，在这种情况下，这个积分的真实值为

\[I=\sqrt{\frac{2}{\pi}}\left(\frac{18}{27}\sqrt{2}\cos\left(4.5\right)-\frac{4}{27}\sqrt{2}\sin\left(4.5\right)+\sqrt{2\pi}\textrm{Si}\left(\frac{3}{\sqrt{\pi}}\right)\right),\]

其中

\[\textrm{Si}\left(x\right)=\int_{0}^{x}\sin\left(\frac{\pi}{2}t^{2}\right)\, dt.\]

是菲涅耳正弦积分。请注意，根据数字计算的积分在 \(1.04\times10^{-11}\) 精度范围内是准确的结果 — 远远低于报告的估计误差。

如果积分函数接收其他参数，它们可以在 args 参数中提供。假设要计算以下积分

\[I(a,b)=\int_{0}^{1} ax^2+b \, dx.\]

可以用以下代码评估此积分

>>> 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 作为参数之一。例如，假设需要计算指数积分

\[E_{n}\left(x\right)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt.\]

的一个数值 (special.expn(n,x) 可以计算出这个积分，但这个事实往往被遗忘)。special.expn 函数的功能可以通过根据 quad 例程定义一个新的函数 vec_expint 来复制

>>> 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获得的积分中可能包含数值误差，错误界限可能会低估错误）。这种情况下积分如下：

\[I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt\, dx=\frac{1}{n}.\]

>>> 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))  # compare against theoretical result
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)

作为非恒定限的示例，请考虑积分

\[I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y \, dx\, dy=\frac{1}{96}.\]

可以使用以下表达式对这个积分进行评估（注意，利用了非恒定 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。积分边界是一个可迭代对象：恒定边界的列表，或非恒定积分边界的函数列表。积分顺序（因此也是边界）是从最内层积分到最外层积分。

上面的积分

\[I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}\, dt\, dx=\frac{1}{n}\]

可以计算为

>>> 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]\) 中。

非恒定积分边界可以用类似的方式处理；上面的示例

\[I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y \, dx\, dy=\frac{1}{96}.\]

可以通过以下方式评估

>>> 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\)（其中 \(k\) 为某个整数），则可以利用 Romberg romb 积分获得积分的高精度估计。Romberg 积分以 2 的幂相关联的步长使用梯形规则，然后对这些估计执行 Richardson 外推，以更大精度逼近积分。

对于任意间隔的样本，可以使用这两个函数 trapezoid 和 simpson。它们分别使用 1 级和 2 级的牛顿-科茨公式来执行积分。梯形规则将相邻点之间的函数逼近为一条直线，而辛普森规则将三个相邻点之间的函数逼近为一个抛物线。

对于等间隔的奇数样本，如果函数是 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=x)
>>> print(I1)
21.0

这与以下公式完全对应

\[\int_{1}^{4} x^2 \, dx = 21,\]

然而，对第二个函数进行积分

>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x=x)
>>> print(I2)
61.5

与以下公式不对应

\[\int_{1}^{4} x^3 \, dx = 63.75\]

因为 f2 中多项式的阶数大于 2。

使用低级回调函数进行更快的积分#

希望减少积分时间的用户可以通过将 C 函数指针传递给 scipy.LowLevelCallable 来获得 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]; /* corresponds to 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 中都可以省略）。请注意，坐标以 double 数组的形式传递，而不是作为一个单独的参数传递。

4.) 现在按照普通方法集成库函数，这里使用 nquad

>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11)

Python 元组以缩减的时间量返回，如预期的那样。所有可选参数都可以与这个方法一同使用，包括指定奇点、无限界限等等。

常微分方程 (`solve_ivp`)#

对给定初始条件的一组常微分方程 (ODE) 进行积分是另一个有用的示例。SciPy 中提供了函数 solve_ivp 用于积分一阶向量微分方程

\[\frac{d\mathbf{y}}{dt}=\mathbf{f}\left(\mathbf{y},t\right),\]

初始条件为 \(\mathbf{y}\left(0\right)=y_{0}\)，其中 \(\mathbf{y}\) 是长度为 \(N\) 的向量，\(\mathbf{f}\) 是从 \(\mathcal{R}^{N}\) 到 \(\mathcal{R}^{N}.\) 的映射。将中间导数引入 \(\mathbf{y}\) 向量中，即可始终将高阶常微分方程降阶为该类型的微分方程。

例如，假设希望求解以下二阶微分方程的解

\[\frac{d^{2}w}{dz^{2}}-zw(z)=0\]

初始条件为 \(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)}.\) 已知，具有这些边界条件的微分方程的解是艾里函数

\[w=\textrm{Ai}\left(z\right),\]

这提供了使用 special.airy 检查积分器的方法。

首先，通过设定 \(\mathbf{y}=\left[\frac{dw}{dz},w\right]\) 和 \(t=z\) 将该 ODE 转换为标准形式。因此，该微分方程变为

\[\begin{split}\frac{d\mathbf{y}}{dt}=\left[\begin{array}{c} ty_{1}\\ y_{0}\end{array}\right]=\left[\begin{array}{cc} 0 & t\\ 1 & 0\end{array}\right]\left[\begin{array}{c} y_{0}\\ y_{1}\end{array}\right]=\left[\begin{array}{cc} 0 & t\\ 1 & 0\end{array}\right]\mathbf{y}.\end{split}\]

换句话说，

\[\mathbf{f}\left(\mathbf{y},t\right)=\mathbf{A}\left(t\right)\mathbf{y}.\]

作为有趣的提示，如果 \(\mathbf{A}\left(t\right)\) 在矩阵乘法下与 \(\int_{0}^{t}\mathbf{A}\left(\tau\right)\, d\tau\) 可交换，则该线性微分方程有一个使用矩阵指数的精确解

\[\mathbf{y}\left(t\right)=\exp\left(\int_{0}^{t}\mathbf{A}\left(\tau\right)d\tau\right)\mathbf{y}\left(0\right),\]

但是，在这种情况下，\(\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)

如果已知函数的 jacobian 矩阵，可以将其传递给 solve_ivp 以取得更好的效果。但请注意默认积分方法 RK45 不支持 jacobian 矩阵，因此必须选择另一种积分方法。一种支持 JacobIan 矩阵的积分方法是以下示例中的 Radau 方法。

>>> def gradient(t, y):
...     return [[0,t], [1,0]]
>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient)

求解带状 Jacobian 矩阵的系统#

odeint 可以得知 Jacobian 是分段的。对于已知僵硬的大型微分方程组，这可以显著提升性能。

作为示例，我们将使用线方法 [MOL] 求解一维 Gray-Scott 偏微分方程。区间 \(x \in [0, L]\) 上函数 \(u(x, t)\) 和 \(v(x, t)\) 的 Gray-Scott 方程为

\[\begin{split}\begin{split} \frac{\partial u}{\partial t} = D_u \frac{\partial^2 u}{\partial x^2} - uv^2 + f(1-u) \\ \frac{\partial v}{\partial t} = D_v \frac{\partial^2 v}{\partial x^2} + uv^2 - (f + k)v \\ \end{split}\end{split}\]

其中 \(D_u\) 和 \(D_v\) 分别为组分 \(u\) 和 \(v\) 的扩散系数，而 \(f\) 和 \(k\) 是常数。（有关该系统的更多信息，请参阅 http://groups.csail.mit.edu/mac/projects/amorphous/GrayScott/）

假定诺伊曼（即“无通量”）边界条件

\[\frac{\partial u}{\partial x}(0,t) = 0, \quad \frac{\partial v}{\partial x}(0,t) = 0, \quad \frac{\partial u}{\partial x}(L,t) = 0, \quad \frac{\partial v}{\partial x}(L,t) = 0\]

要应用线法，我们通过把 \(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\) 导数。即

\[\frac{\partial^2 u}{\partial x^2}(x_j, t) \rightarrow \frac{u_{j-1}(t) - 2 u_{j}(t) + u_{j+1}(t)}{(\Delta x)^2}\]

然后我们得到一个 \(2N\) 常微分方程组

(1)#\[\begin{split} \begin{split} \frac{du_j}{dt} = \frac{D_u}{(\Delta x)^2} \left(u_{j-1} - 2 u_{j} + u_{j+1}\right) -u_jv_j^2 + f(1 - u_j) \\ \frac{dv_j}{dt} = \frac{D_v}{(\Delta x)^2} \left(v_{j-1} - 2 v_{j} + v_{j+1}\right) + u_jv_j^2 - (f + k)v_j \end{split}\end{split}\]

为方便起见，忽略了 \((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)\) 按类似方式定义。

然后

(2)#\[\begin{split} \begin{split} \frac{du_0}{dt} = \frac{D_u}{(\Delta x)^2} \left(2u_{1} - 2 u_{0}\right) -u_0v_0^2 + f(1 - u_0) \\ \frac{dv_0}{dt} = \frac{D_v}{(\Delta x)^2} \left(2v_{1} - 2 v_{0}\right) + u_0v_0^2 - (f + k)v_0 \end{split}\end{split}\]

以及

(3)#\[\begin{split} \begin{split} \frac{du_{N-1}}{dt} = \frac{D_u}{(\Delta x)^2} \left(2u_{N-2} - 2 u_{N-1}\right) -u_{N-1}v_{N-1}^2 + f(1 - u_{N-1}) \\ \frac{dv_{N-1}}{dt} = \frac{D_v}{(\Delta x)^2} \left(2v_{N-2} - 2 v_{N-1}\right) + u_{N-1}v_{N-1}^2 - (f + k)v_{N-1} \end{split}\end{split}\]

我们完整的\(2N\)常微分方程组对于\(k = 1, 2, \ldots, N-2\)是(1)，加上(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}\}\)对变量进行排序时，雅可比矩阵的非零元素的模式是

\[\begin{split}\begin{smallmatrix} * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 \\ 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 \\ 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 \\ 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * \\ * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 \\ 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 \\ 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 \\ 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & * & * & * \\ \end{smallmatrix}\end{split}\]

雅可比模式，它交错变量为 \(\{u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}\}\) 是

\[\begin{split}\begin{smallmatrix} * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * \\ \end{smallmatrix}\end{split}\]

在两种情况下，只有五个非平凡对角线，但当变量交错时，带宽会小得多。也就是说，主对角线和主对角线正上方和正下方的两个对角线是非零对角线。这一点非常重要，因为 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):
    """
    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 雅可比矩阵是带状的。这允许 underlying 求解器 (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

参考资料#

[MOL]

https://zh.wikipedia.org/wiki/%E7%BA%BF%E9%97%B4%E5%88%B7%E6%96%B9%E6%B3%95

整合（scipy.integrate）#

一般积分（quad）#

一般多重积分(dblquad, tplquad, nquad)#