Python SciPy integrate.solve_bvp用法及代码示例

本文简要介绍 python 语言中 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)#

求解 ODE 系统的边值问题。

此函数在数值上求解服从 two-point 边界条件的 ODE 的一阶系统：

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

这里 x 是一维自变量，y(x) 是 N-D vector-valued 函数，p 是与 y(x) 一起找到的未知参数的 k-D 向量。对于要确定的问题，必须有 n + k 个边界条件，即 bc 必须是 (n + k)-D 函数。

系统右侧的最后一个单数术语是可选的。它由 n-by-n 矩阵 S 定义，使得解必须满足 S y(a) = 0。此条件将在迭代过程中强制执行，因此不得与边界条件相矛盾。请参阅 [2] 了解在数值求解 BVP 时如何处理该术语的说明。

复杂领域中的问题也可以得到解决。在这种情况下，y 和 p 被认为是复数，f 和 bc 被假定为 complex-valued 函数，但 x 保持实数。请注意，f 和 bc 必须是复可微的(满足 Cauchy-Riemann 方程 [4])，否则您应该分别针对实部和虚部重写您的问题。要解决复杂域中的问题，请使用复杂数据类型传递 y 的初始猜测(见下文)。

参数 ：：

fun： 可调用的

系统的右侧。调用签名为 fun(x, y) ，如果存在参数，则为 fun(x, y, p) 。所有参数都是 ndarray：形状为 (m,) 的 x、形状为 (n, m) 的 y，这意味着 y[:, i] 对应于 x[i] ，形状为 (k,) 的 p 。返回值必须是形状为 (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 (默认)，则假定问题不依赖于任何参数。

S： 数组 形状为 (n, n) 或无

定义单数项的矩阵。如果无(默认)，则无需单数项即可解决问题。

fun_jac： 可调用或无，可选

函数计算 f 关于 y 和 p 的导数。调用签名是 fun_jac(x, y) 或 fun_jac(x, y, p) 如果参数存在。返回必须按以下顺序包含 1 或 2 个元素：

• df_dy : array_like with shape (n, n, m), where an element (i, j, q) equals to d f_i(x_q, y_q, p) / d (y_q)_j.

• df_dp : array_like with shape (n, k, m), where an element (i, j, q) equals to d f_i(x_q, y_q, p) / d p_j.

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

如果fun_jac 为无(默认)，则将通过前向有限差分来估计导数。

bc_jac： 可调用或无，可选

函数计算 bc 关于 ya、yb 和 p 的导数。调用签名是 bc_jac(ya, yb) 或 bc_jac(ya, yb, p) 如果参数存在。返回必须按以下顺序包含 2 或 3 个元素：

• dbc_dya : array_like with shape (n, n), where an element (i, j) equals to d bc_i(ya, yb, p) / d ya_j.

• dbc_dyb : array_like with shape (n, n), where an element (i, j) equals to d bc_i(ya, yb, p) / d yb_j.

• dbc_dp : array_like with shape (n, k), where an element (i, j) equals to d bc_i(ya, yb, p) / d p_j.

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

如果bc_jac 为无(默认)，则将通过前向有限差分来估计导数。

tol： 浮点数，可选

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

max_nodes： 整数，可选

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

verbose： {0, 1, 2}，可选

算法的详细程度：

• 0 (default) : work silently.

• 1 : display a termination report.

• 2 : display progress during iterations.

bc_tol： 浮点数，可选

边界条件残差所需的绝对容差：公元前值应满足abs(bc) < bc_tolcomponent-wise。等于tol默认。最多允许 10 次迭代来实现此容差。

返回 ：：

定义了以下字段的 Bunch 对象：

sol： 聚苯胺

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

p： ndarray 或无，形状(k，)

找到参数。无，如果问题中不存在参数。

x： ndarray，形状(m，)

最终网格的节点。

y： ndarray，形状(n，m)

网格节点处的解值。

yp： ndarray，形状(n，m)

网格节点处的解导数。

rms_residuals： ndarray，形状(m - 1，)

每个网格间隔上的相对残差的 RMS 值(参见 tol 参数的说明)。

niter： int

完成的迭代次数。

status： int

算法终止原因：

• 0: The algorithm converged to the desired accuracy.

• 1: The maximum number of mesh nodes is exceeded.

• 2: A singular Jacobian encountered when solving the collocation system.

message： string

终止原因的口头说明。

success： bool

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

注意：

该函数实现了一个 4 阶搭配算法，其残差控制类似于 [1]。如[3] 中所述，搭配系统通过具有affine-invariant 准则函数的阻尼牛顿法求解。

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

参考：

[1] (1,2)

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

[2]

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

[3]

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

[4]

Cauchy-Riemann equations 在维基百科上。

例子：

在第一个示例中，我们解决了 Bratu 问题：

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

对于 k = 1。

我们将方程重写为 first-order 系统并实现其右侧求值：

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，非平凡解 y = A * sin(k * x) 是可能的，其中 n 是整数。为了建立归一化常数 A = 1，我们添加了一个边界条件：

y'(0) = k

再次，我们将方程重写为 first-order 系统并实现其右侧评估：

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()