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)[源代码][源代码]#
求解常微分方程系统的边界值问题。
此函数数值求解一个一阶ODE系统，满足两点边界条件:
dy / dx = f(x, y, p) + S * y / (x - a), a <= x <= b bc(y(a), y(b), p) = 0
这里 x 是一个 1 维自变量，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(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,) 的数组。- xarray_like, 形状 (m,)
初始网格。必须是一个严格递增的实数序列，且
x[0]=a和x[-1]=b。- yarray_like, 形状 (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) 的 array_like，其中元素 (i, j, q) 等于 d f_i(x_q, y_q, p) / d (y_q)_j。
df_dp : 形状为 (n, k, m) 的 array_like，其中元素 (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) 的 array_like，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d ya_j。
dbc_dyb : 形状为 (n, n) 的 array_like，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d yb_j。
dbc_dp : 形状为 (n, k) 的 array_like，其中元素 (i, j) 等于 d bc_i(ya, yb, p) / d p_j。
如果问题在没有未知参数的情况下得到解决，则不应返回 dbc_dp。
如果 bc_jac 为 None（默认），导数将通过前向有限差分法估计。
- tolfloat, 可选
解的期望容差。如果我们定义
r = y' - f(x, y)，其中 y 是找到的解，那么求解器在每个网格区间上尝试实现norm(r / (1 + abs(f)) < tol，其中norm是以均方根意义估计的（使用数值积分公式）。默认值为 1e-3。- max_nodesint, 可选
网格节点的最大允许数量。如果超过此数量，算法将终止。默认值为1000。
- 详细{0, 1, 2}, 可选
算法详细程度的级别：
0 (默认) : 静默工作。
1 : 显示终止报告。
2 : 在迭代过程中显示进度。
- bc_tolfloat, 可选
边界条件残差的期望绝对容差：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,)
每个网格区间上相对残差的均方根值（参见 tol 参数的描述）。
- niter整数
已完成迭代的次数。
- 状态整数
算法终止的原因：
0: 算法已收敛到所需的精度。
1: 网格节点数量超过最大限制。
2: 在求解配置系统时遇到奇异雅可比矩阵。
- 消息字符串
终止原因的口头描述。
- 成功布尔
如果算法收敛到所需的精度（
status=0），则为真。
注释
此函数实现了一个四阶的配点算法，其残差控制类似于 [1]。配点系统通过一个阻尼牛顿方法求解，该方法使用了一个仿射不变的准则函数，如 [3] 中所述。
注意，在 [1] 中，积分残差是未按区间长度归一化定义的。因此，它们的定义与这里使用的定义相差一个 h**0.5 的倍数（h 是区间长度）。
Added in version 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 和 H. Xu，“一个用户友好的 Fortran BVP 求解器”。
[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()
我们看到这两种解决方案具有相似的形状，但在规模上有显著差异。
在第二个例子中，我们解决了一个简单的 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()
