用 scipy 的 solve_bvp 求解 BVP

Posted 2023-03-12

【中文标题】用 scipy 的 solve_bvp 求解 BVP

【英文标题】：Solving a BVP with scipy's solve_bvp

【发布时间】：2017-12-05 15:16:20

【问题描述】：

我有一个具有 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)

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

BVP 的 MATLAB 解决方案

【参考方案1】：

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

您可以通过检查sol.status 来确定@​​987654327@ 是否收敛失败。如果不为 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()

【讨论】：

有没有办法知道这是否是真正的答案？就像在 MATLAB 中一样，如果我的起始猜测错误，我会得到一个错误，但这里一切都给了我 some 答案。所以我现在才确定如何验证我的解决方案在 Python 中是对还是错。我假设没有诸如溢出警告之类的错误消息意味着它有效。我刚刚尝试使用np.linspace(700, 400, n)])，它似乎也失败了。我想我需要在某种程度上包含变量的范围？

"我刚用 np.linspace(700, 400, n)])" 试过了 不是必须的，当然也不能保证收敛，但是好像传递满足边界条件的初始猜测是个好主意。

谢谢。这是验证它的好方法。您使用 y[2] 初始猜测作为从 800 到 200 的 linspace 数组，而不是反过来，因为它就是这样，对吧？你为什么不为y[1] 做一个从600 到更高的linspace，而只使用full_like 代替？ y[2] 更敏感吗？你以前是怎么知道的？

“你使用 y[2] 初始猜测作为从 800 到 200 的 linspace 数组，而不是反过来，因为它是这样的，对吧？” 是的，并且边界条件是y[2] x=1 是 200。我没有使用 linspace 代替 y[1] 因为我不需要 - 一旦我使 y[0] 增加和 y[2] 减少，求解器可靠收敛。但是，我只是尝试将np.linspace(673, 800, n) 用于y[1]，但它未能收敛！去搞清楚。较小的坡度有效，例如np.linspace(673, 700, n).

我明白了。那很有意思。我想在做这些问题时要记住保持低斜率。谢谢。