用scipy的resolve_bvp解决BVP [英] Solving a BVP with scipy's solve_bvp

问题描述

我有一个包含3个边界条件的3个微分方程组(从我相信的代码中将显而易见).我设法在MATLAB中使用循环来解决该问题，以逐个更改初始猜测，而不会在程序即将返回错误时终止程序.但是，在scipy的solve_bvp上，我总是得到 some 答案，尽管这是错误的.因此，我不断改变自己的猜测(不断改变答案)，并给出了与实际解决方案相当接近的数字，但仍然无法正常工作.也许由于代码不起作用，代码是否还有其他问题?我刚刚编辑了他们文档的代码.

I have a system of 3 differential equations (will be obvious from the code I believe) with 3 boundary conditions. I managed to solve it in MATLAB with a loop to change the initial guess bit by bit without terminating the program if it is about to return an error. However, on scipy's solve_bvp, I always get some answer, although it is wrong. So I kept changing my guesses (which kept changing the answer) and am giving pretty close numbers to what I have from the actual solution and it's still not working. Is there some other problem with the code perhaps, due to which it's not working? I just edited their documentation's code.

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)

实际解决方案在下面的图中给出.

The actual solution is given below in the figure.

MATLAB BVP解决方案

MATLAB Solution to the BVP

推荐答案

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

Apparently you need a better initial guess, otherwise the iterative method used by solve_bvp can create values in y[1] that make the expression exp(-19846/y[1]) overflow. When that happens, the algorithm is likely to fail. An overflow in that expression means that some value in y[1] is negative; i.e. the solver is so far out in the weeds that it has little chance of converging to a correct solution. You'll see warnings, and sometimes the function still returns a reasonable solution, but usually it returns garbage when the overflow occurs.

您可以通过检查sol.status来确定solve_bvp是否收敛失败.如果不为0，则失败. sol.message包含描述状态的短信.

You can determine whether or not solve_bvp has failed to converge by checking sol.status. If it is not 0, something failed. sol.message contains a text message describing the status.

我可以使用此方法创建初始猜测，从而获得Matlab解决方案:

I was able to get the Matlab solution by using this to create the initial guess:

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

n的较小值也可以使用，但是当n太小时，会出现溢出警告.

Smaller values of n also work, but when n is too small, an overflow warning can appear.

这是我的脚本修改版，后面是脚本生成的图:

Here's my modified version of your script, followed by the plot that it generates:

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

这篇关于用scipy的resolve_bvp解决BVP的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！