使用 scipy 的 solve_bvp 解决 BVP

Question

我有一个包含 3 个边界条件的 3 个微分方程系统（从我相信的代码中可以很明显地看出）。我设法在 MATLAB 中用一个循环来解决它，以便在程序即将返回错误时一点一点地改变初始猜测而不终止程序。然而，在scipy's 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)

实际解决方案如下图所示。

BVP 的 MATLAB 解决方案

Answer 1

显然你需要一个更好的初始猜测，否则使用的迭代方法solve_bvp可以创建值y[1]，使表达式exp(-19846/y[1])溢出。当这种情况发生时，算法很可能会失败。该表达式中的溢出意味着某些值y[1]是负数；即求解器在杂草中是如此遥远，以至于它几乎没有机会收敛到正确的解决方案。您会看到警告，有时该函数仍会返回合理的解决方案，但通常会在发生溢出时返回垃圾。

您可以solve_bvp通过检查来确定是否收敛失败sol.status。如果它不是 0，则表示失败。 sol.message包含描述状态的文本消息。

我能够通过使用它来创建初始猜测来获得 Matlab 解决方案：

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

较小的值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()