SciPy: solve_bvp问题二阶差分。情商

Question

试图解决一个二阶差异。情商。在两个边界条件下，我似乎什么都没有尝试，而且我找不到一个教程，它包含了所有/类似于我表达式中的所有/类似的术语，而且至少对我来说，参与文档并不能真正清楚地解释如何使用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)

谢谢!

==========EDIT=========================有一个解析解，我已经计算过了，它只是看我是否也能在数值上找到同样的解，我认为主要的问题是，解有两个独立的区域，一个是rR (out)。这导致了两个不同的解(仅在各自域中有效)，条件是y_in(R) == y_out(R)和y_in'(R) == y_out'(R)。全2部分解，其中半径= 1，a=99，b= 1，常数=1，y(inf) =常数

从Lehmann的解决方案，它得到了正确的形状(至少在内部区域，虽然不是在正确的尺度)。

我只是不确定你将如何编码所有的等价物解决方案，我想，甚至从一开始就得到了它们的解决方案，尽管Lutz的答案是正确方向上的一个惊人的点。谢谢

Answer 1

问题

• 方程的阶数为2，状态向量的维数为2，其值总是y[0]，导数y[1]，没有y[2]，很可能是由Matlab的平移而来的。
• 在边界条件下，不存在ya[2]，导数值为ya[1]，第二个条件中的函数值为yb[0]。
• 初始解猜测必须有相同数目的状态分量2。
• 为什么要用相同的数据计算两种解决方案？
• 备注:不需要将返回值转换为numpy类型，求解器将检查和转换。

带奇异处理的BVP框架

在r=0中，ODE是单数的，所以必须以一种特殊的方式来处理第一个段。中值定理

(y'(r)-y'(0))/r->y''(0)  for  r->0,

所以在这个限制下，你可以得到r->0

3*y''(0) = a*y(0) + b*y(0)^3`. 

这允许将第一个弧定义为

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

按照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 + y0*(a+b*y0*y0)*h*h/6;
    dyh = y0*(a+b*y0*y0)*h/3
    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,f"y(0)={sol.p[0]}");
plt.plot(sol.x, sol.y[0]);

这样，使用示例参数a, b, c = -1, 0.2, 3，您就可以得到一个具有y(0)=2.236081087849196和结果的图的收敛求解器调用。

​

​