scipy.integrate.

solve_bvp#

scipy.integrate.solve_bvp(fun, bc, x, y, p=None, S=None, fun_jac=None, bc_jac=None, tol=0.001, max_nodes=1000, verbose=0, bc_tol=None)[source]#

求解常微分方程组的边值问题。

此函数数值求解一阶常微分方程组，该方程组满足两点边界条件

dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b
bc(y(a), y(b), p) = 0

这里 x 是一个一维自变量，y(x) 是一个 N 维向量值函数，p 是一个 k 维未知参数向量，它需要和 y(x) 一起求解。为了确定问题，必须有 n + k 个边界条件，即 bc 必须是一个 (n + k) 维函数。

系统右手边的最后一个奇异项是可选的。它由一个 n × n 矩阵 S 定义，使得解必须满足 S y(a) = 0。此条件将在迭代过程中强制执行，因此它不应与边界条件相矛盾。有关在数值求解边值问题时如何处理此项的说明，请参见 [2]。

也可以求解复域中的问题。在这种情况下，y 和 p 被认为是复数，f 和 bc 被认为是复数值函数，但 x 保持实数。请注意，f 和 bc 必须是复可微的（满足柯西-黎曼方程 [4]），否则您应该分别为实部和虚部重写问题。要求解复域中的问题，请传递具有复数数据类型的 y 的初始猜测（见下文）。

参数：
fun可调用对象

系统的右手边。调用签名为 fun(x, y)，或者如果存在参数，则为 fun(x, y, p)。所有参数都是 ndarray：x 形状为 (m,)，y 形状为 (n, m)，这意味着 y[:, i] 对应于 x[i]，p 形状为 (k,)。返回值必须是一个形状为 (n, m) 的数组，布局与 y 相同。

bc可调用对象

计算边界条件残差的函数。调用签名为 bc(ya, yb)，或者如果存在参数，则为 bc(ya, yb, p)。所有参数都是 ndarray：ya 和 yb 形状为 (n,)，p 形状为 (k,)。返回值必须是一个形状为 (n + k,) 的数组。

x类数组，形状为 (m,)

初始网格。必须是一个严格递增的实数序列，其中 x[0]=a 和 x[-1]=b。

y类数组，形状为 (n, m)

网格节点上函数值的初始猜测，第 i 列对应于 x[i]。对于复域中的问题，请传递具有复数数据类型的 y（即使初始猜测是纯实数）。

p形状为 (k,) 的类数组或 None，可选

未知参数的初始猜测。如果为 None（默认值），则假设问题不依赖于任何参数。

S形状为 (n, n) 的类数组或 None

定义奇异项的矩阵。如果为 None（默认值），则在没有奇异项的情况下求解问题。

fun_jac可调用对象或 None，可选

计算 f 相对于 y 和 p 的导数的函数。调用签名为 fun_jac(x, y)，或者如果存在参数，则为 fun_jac(x, y, p)。返回值必须包含以下顺序的 1 或 2 个元素

  • df_dy : 类数组，形状为 (n, n, m)，其中元素 (i, j, q) 等于 d f_i(x_q, y_q, p) / d (y_q)_j。

  • df_dp : 类数组，形状为 (n, k, m)，其中元素 (i, j, q) 等于 d f_i(x_q, y_q, p) / d p_j。

这里 q 对定义 x 和 y 的节点进行编号，而 i 和 j 对向量分量进行编号。如果在没有未知参数的情况下求解问题，则不应返回 df_dp。

如果 fun_jac 为 None（默认值），则导数将通过向前有限差分估计。

bc_jac可调用对象或 None，可选

计算 bc 相对于 ya、yb 和 p 的导数的函数。调用签名为 bc_jac(ya, yb)，或者如果存在参数，则为 bc_jac(ya, yb, p)。返回值必须包含以下顺序的 2 或 3 个元素

  • dbc_dya : 类数组，形状为 (n, n)，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d ya_j。

  • dbc_dyb : 类数组，形状为 (n, n)，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d yb_j。

  • dbc_dp : 类数组，形状为 (n, k)，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d p_j。

如果在没有未知参数的情况下求解问题，则不应返回 dbc_dp。

如果 bc_jac 为 None（默认值），则导数将通过向前有限差分估计。

tol浮点数，可选

所需的解容差。如果我们定义 r = y' - f(x, y)，其中 y 是找到的解，那么求解器尝试在每个网格间隔上实现 norm(r / (1 + abs(f)) < tol，其中 norm 是使用数值积分公式以均方根意义估计的。默认值为 1e-3。

max_nodes整数，可选

允许的最大网格节点数。如果超过，算法将终止。默认值为 1000。

verbose{0, 1, 2}，可选

算法的详细程度

  • 0（默认值）：静默工作。

  • 1：显示终止报告。

  • 2：在迭代过程中显示进度。

bc_tol浮点数，可选

所需的边界条件残差绝对容差：bc 值应满足 abs(bc) < bc_tol 按分量计算。默认情况下等于 tol。最多允许 10 次迭代来实现此容差。

返回值：
具有以下定义字段的 Bunch 对象
solPPoly

找到的 y 解，作为 scipy.interpolate.PPoly 实例，一个 C1 连续三次样条。

pndarray 或 None，形状为 (k,)

找到的参数。如果问题中不存在参数，则为 None。

xndarray，形状为 (m,)

最终网格的节点。

yndarray，形状为 (n, m)

网格节点上的解值。

ypndarray，形状为 (n, m)

网格节点上的解导数。

rms_residualsndarray，形状为 (m - 1,)

每个网格间隔上相对残差的 RMS 值（参见 tol 参数的描述）。

niter整数

已完成的迭代次数。

status整数

算法终止的原因

  • 0：算法收敛到所需的精度。

  • 1：超过最大网格节点数。

  • 2：在求解配置系统时遇到奇异雅可比矩阵。

message字符串

终止原因的文字描述。

success布尔值

如果算法收敛到所需的精度 (status=0)，则为 True。

注意

该函数实现了一个类似于[1]的四阶配置算法，并控制残差。配置系统通过一个阻尼牛顿法求解，该方法使用仿射不变准则函数，如[3]所述。

请注意，在[1]中，积分残差定义时没有按区间长度归一化。因此，它们的定义与这里使用的定义不同，相差一个 h**0.5 的乘数（h 是区间长度）。

在版本 0.18.0 中添加。

参考文献

[1] (1,2)

J. Kierzenka, L. F. Shampine, “基于残差控制和 Maltab PSE 的 BVP 求解器”，ACM Trans. Math. Softw.，第 27 卷，第 3 期，第 299-316 页，2001 年。

[2]

L.F. Shampine、P. H. Muir 和 H. Xu，“用户友好的 Fortran BVP 求解器”。

[3]

U. Ascher、R. Mattheij 和 R. Russell，“常微分方程边值问题的数值解”。

示例

在第一个示例中，我们求解 Bratu 问题

y'' + k * exp(y) = 0
y(0) = y(1) = 0

对于 k = 1。

我们将方程改写为一阶系统，并实现其右手边求值

y1' = y2
y2' = -exp(y1)

>>> import numpy as np
>>> def fun(x, y):
...     return np.vstack((y[1], -np.exp(y[0])))

实现边界条件残差的求值

>>> def bc(ya, yb):
...     return np.array([ya[0], yb[0]])

定义具有 5 个节点的初始网格

>>> x = np.linspace(0, 1, 5)

众所周知，此问题有两个解。为了获得这两个解，我们对 y 使用两个不同的初始猜测。我们用下标 a 和 b 表示它们。

>>> y_a = np.zeros((2, x.size))
>>> y_b = np.zeros((2, x.size))
>>> y_b[0] = 3

现在我们准备运行求解器。

>>> from scipy.integrate import solve_bvp
>>> res_a = solve_bvp(fun, bc, x, y_a)
>>> res_b = solve_bvp(fun, bc, x, y_b)

让我们绘制两个找到的解。我们利用将解以样条形式表示的优势来生成平滑的曲线。

>>> x_plot = np.linspace(0, 1, 100)
>>> y_plot_a = res_a.sol(x_plot)[0]
>>> y_plot_b = res_b.sol(x_plot)[0]
>>> import matplotlib.pyplot as plt
>>> plt.plot(x_plot, y_plot_a, label='y_a')
>>> plt.plot(x_plot, y_plot_b, label='y_b')
>>> plt.legend()
>>> plt.xlabel("x")
>>> plt.ylabel("y")
>>> plt.show()
../../_images/scipy-integrate-solve_bvp-1_00_00.png

我们看到两个解的形状相似，但在尺度上存在显著差异。

在第二个示例中，我们求解一个简单的斯图姆-刘维尔问题

y'' + k**2 * y = 0
y(0) = y(1) = 0

众所周知，对于 k = pi * n，其中 n 是整数，非平凡解 y = A * sin(k * x) 是可能的。为了建立归一化常数 A = 1，我们添加一个边界条件

y'(0) = k

同样，我们将方程改写为一阶系统，并实现其右手边求值

y1' = y2
y2' = -k**2 * y1

>>> def fun(x, y, p):
...     k = p[0]
...     return np.vstack((y[1], -k**2 * y[0]))

请注意，参数 p 作为向量传递（在本例中只有一个元素）。

实现边界条件

>>> def bc(ya, yb, p):
...     k = p[0]
...     return np.array([ya[0], yb[0], ya[1] - k])

设置初始网格和 y 的猜测。我们的目标是找到 k = 2 * pi 的解，为此，我们将 y 的值设置为近似遵循 sin(2 * pi * x)

>>> x = np.linspace(0, 1, 5)
>>> y = np.zeros((2, x.size))
>>> y[0, 1] = 1
>>> y[0, 3] = -1

使用 6 作为 k 的初始猜测运行求解器。

>>> sol = solve_bvp(fun, bc, x, y, p=[6])

我们看到找到的 k 大致正确

>>> sol.p[0]
6.28329460046

最后，绘制解以查看预期的正弦曲线

>>> x_plot = np.linspace(0, 1, 100)
>>> y_plot = sol.sol(x_plot)[0]
>>> plt.plot(x_plot, y_plot)
>>> plt.xlabel("x")
>>> plt.ylabel("y")
>>> plt.show()
../../_images/scipy-integrate-solve_bvp-1_01_00.png
本页