SciPy：solve_bvp问题2阶差。式

Question

试图解决二阶差异。 eq。有两个边界条件，但我尝试的任何方法似乎都没有用，而且我找不到包含与表达式中所用术语完全/相似的术语的教程，至少对我而言，scipy文档并未真正说明如何使用显然是solve_bvp。

我有等式：y''+ 2 / r * y'= A * y + b * y ^ 3，其中y是r的函数。

我将其重写如下：

y1 = y（r）

y2 = y1'

使得y2'= -2 / r * y2 + y1（A + b * y1 ^ 2）

[当y'（0）= 0时，y（r = 10）=常数

A，b和常数是已知的。

我有以下代码，但是它似乎不起作用，并且如上所述，我对文档有些困惑，因此不胜感激！

def fun(x, y):
    return np.vstack((y[2], -2/x*y[1]+y[0]*(A+b*y[0]*y[0])))

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


x = np.linspace(1, 10, 10)
ya = np.zeros((3, x.size))
yb = np.zeros((3, x.size))

sol_1 = solve_bvp(fun, bc, x, ya)
sol_2 = solve_bvp(fun, bc, x, yb)

谢谢！

Answer 1

问题

• 等式的阶次为2，因此状态向量的维数为2，值始终为y[0]，导数y[1]，没有y[2]，可能是从Matlab。

• 在边界条件下，也没有ya[2]，微分值为ya[1]，第二个函数值为yb[0]。

• 最初的解决方案猜测必须具有相同数量的2个状态分量。

• 为什么要用相同的数据计算两个解？

• 注：不必将返回值转换为numpy类型，求解器仍会进行检查和转换。

具有奇异性处理的BVP框架

ODE在r=0处是奇异的，因此必须以特殊方式处理第一段。平均值定理为(y'(r)-y'(0))/r->y''(0)给出r->0，因此在该极限r->0中得到3*y''(0) = a*y(0) + b*y(0)^3。这允许将第一个弧定义为

y(r) = y0 + 0.5*(a*y0 + b*y0^3)*r^2
y'(r) = (a*y0 + b*y0^3)*r

最高订购r³和r²。因此，如果要在1e-9中设置精度y(r)，则第一段不应长于1e-3。

而不是尝试从y0和y(h)的方程式中消除y'(h)以获得连接ya[0]和ya[1]的条件，让求解器也进行这项工作，并将y0作为参数添加到系统。然后边界条件具有3个与参数添加的虚拟维相对应的时隙，可以自然地用等式y(h)=ya[0]和ya[1]=y'(h)以及右边界条件填充。

您总共可以将系统定义为

h = 1e-3;

def fun(r, y, p):
    return  y[1], -2/r*y[1]+y[0]*(a+b*y[0]*y[0]) 

def bc(ya, yb, p):
    y0, = p
    yh = y0 + 0.5*y0*(a+b*y0*y0)*h*h;
    dyh = y0*(a+b*y0*y0)*h
    return ya[0]-yh, ya[1]-dyh, yb[0]-c


x = np.linspace(h, 10, 10)
ya = np.zeros((2, x.size))

sol = solve_bvp(fun, bc, x, ya, p=[1])
print(sol.message);
plt.plot(sol.x, sol.y[0]);
因此，使用示例参数a, b, c = -1, 0.2, 3，您将得到带有结果图的收敛解算器调用