用scipy解決一個BVP solve_bvp

Question

我有一個3個微分方程組（從我認爲的代碼中可以看出）具有3個邊界條件。我設法在MATLAB中用一個循環來解決這個問題，如果它要返回一個錯誤，那麼一點一點地改變最初的猜測而不終止程序。不過，在scipy的solve_bvp，我總是得到一些的答案，雖然是錯的。所以我一直在改變我的猜測（它不斷改變答案），並且對我從實際解決方案中得到的數據給出了非常接近的數字，但它仍然不起作用。也許還有其他一些代碼問題，由於它不起作用？我剛剛編輯了他們的文檔代碼。

import numpy as np 
def fun(x, y): 
    return np.vstack((3.769911184e12*np.exp(-19846/y[1])*(1-y[0]), 0.2056315191*(y[2]-y[1])+6.511664773e14*np.exp(-19846/y[1])*(1-y[0]), 1.696460033*(y[2]-y[1]))) 
def bc(ya, yb): 
    return np.array([ya[0], ya[1]-673, yb[2]-200]) 
x = np.linspace(0, 1, 5) 
#y = np.ones((3, x.size)) 
y = np.array([[1, 1, 1, 1, 1], [670, 670, 670, 670, 670], [670, 670, 670, 670, 670] ]) 
from scipy.integrate import solve_bvp 
sol = solve_bvp(fun, bc, x, y) 

實際的解決方案如下圖所示。

MATLAB解法BVP

Answer 1

顯然，你需要一個更好的初始猜測，否則solve_bvp使用迭代方法可以y[1]，使表達exp(-19846/y[1])溢出創造價值。發生這種情況時，算法可能會失敗。該表達式中的溢出意味着y[1]中的某個值爲負值;即解算器在雜草中如此之遠以至於它很少有機會收斂到正確的解決方案。你會看到警告，有時候這個函數仍然會返回一個合理的解決方案，但是通常當溢出發生時它會返回垃圾。

您可以通過檢查sol.status來確定solve_bvp是否收斂失敗。如果它不是0，則失敗。 sol.message包含描述狀態的文本消息。

我能夠用它來創建最初的猜測得到了Matlab的解決方案：的n

n = 25 
x = np.linspace(0, 1, n) 
y = np.array([x, np.full_like(x, 673), np.linspace(800, 200, n)]) 

較小的值也行，但是當n太小，溢出的警告可能會出現。

這裏是我的腳本的修改版本，之後的情節，它產生：

import numpy as np 
from scipy.integrate import solve_bvp 
import matplotlib.pyplot as plt 


def fun(x, y): 
    t1 = np.exp(-19846/y[1])*(1 - y[0]) 
    dy21 = y[2] - y[1] 
    return np.vstack((3.769911184e12*t1, 
         0.2056315191*dy21 + 6.511664773e14*t1, 
         1.696460033*dy21)) 

def bc(ya, yb): 
    return np.array([ya[0], ya[1] - 673, yb[2] - 200]) 


n = 25 
x = np.linspace(0, 1, n) 
y = np.array([x, np.full_like(x, 673), np.linspace(800, 200, n)]) 

sol = solve_bvp(fun, bc, x, y) 

if sol.status != 0: 
    print("WARNING: sol.status is %d" % sol.status) 
print(sol.message) 

plt.subplot(2, 1, 1) 
plt.plot(sol.x, sol.y[0], color='#801010', label='$y_0(x)$') 
plt.grid(alpha=0.5) 
plt.legend(framealpha=1, shadow=True) 
plt.subplot(2, 1, 2) 
plt.plot(sol.x, sol.y[1], '-', color='C0', label='$y_1(x)$') 
plt.plot(sol.x, sol.y[2], '--', color='C0', label='$y_2(x)$') 
plt.xlabel('$x$') 
plt.grid(alpha=0.5) 
plt.legend(framealpha=1, shadow=True) 
plt.show()