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)[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。此條件將在迭代期間強制執行，因此不得與邊界條件相矛盾。有關數值求解 BVP 時如何處理此項的說明，請參閱 [2]。

複數域中的問題也可以求解。在這種情況下，y 和 p 被視為複數，並且 f 和 bc 被假定為複數值函數，但 x 保持實數。請注意，f 和 bc 必須是複數可微分的（滿足柯西-黎曼方程式 [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：形狀為 (n,) 的 ya 和 yb，以及形狀為 (k,) 的 p。傳回值必須是形狀為 (n + k,) 的陣列。

xarray_like，形狀為 (m,)

初始網格。必須是嚴格遞增的實數序列，其中 x[0]=a 且 x[-1]=b。

yarray_like，形狀為 (n, m)

網格節點處函數值的初始猜測值，第 i 列對應於 x[i]。對於複數域中的問題，傳遞具有複數資料類型的 y（即使初始猜測值純粹是實數）。

parray_like，形狀為 (k,) 或 None，選用

未知參數的初始猜測值。如果為 None（預設值），則假定問題不依賴於任何參數。

Sarray_like，形狀為 (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。

verbose{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,)

每個網格區間上相對殘差的 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”.

範例

在第一個範例中，我們求解布拉圖問題

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
在本頁中