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