Ctrl+K
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)[源代码]#

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

此函数数值求解服从两点边值条件的一阶常微分方程组

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。此条件将在迭代期间强制执行，因此它不得与边界条件相矛盾。有关在数值求解 BVP 时如何处理此项的解释，请参见[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 参数的描述）。

niterint

已完成的迭代次数。

statusint

算法终止的原因

  • 0: 算法已收敛到期望的精度。

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

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

messagestring

终止原因的文字描述。

successbool

如果算法收敛到期望的精度（status=0），则为 True。

备注

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

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

在 0.18.0 版本中添加。

参考文献

[1] (1,2)

J. Kierzenka, L. F. Shampine, “A BVP Solver Based on Residual Control and the Maltab PSE”, ACM Trans. Math. Softw., Vol. 27, Number 3, pp. 299-316, 2001.

[2]

L.F. Shampine, P. H. Muir and H. Xu, “A User-Friendly Fortran BVP Solver”.

[3]

U. Ascher, R. Mattheij and R. Russell “Numerical Solution of Boundary Value Problems for Ordinary Differential Equations”.

示例

在第一个示例中，我们求解 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

我们看到这两个解的形状相似，但在尺度上差异很大。

在第二个示例中，我们求解一个简单的 Sturm-Liouville 问题

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
本页内容