SciPy: solve_bvp 问题二阶差分。方程式

Question

正在尝试解决二阶差异。当量。有 2 个边界条件，但我尝试的任何方法似乎都不起作用，而且我找不到包含 all/similar 术语的教程，至少对我来说，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)

谢谢！

==========编辑========================= 有一个解析解，我已经计算过了，只是看看我是否也能在数值上找到相同的解，我认为主要问题是解有两个独立的区域，其中 r < R (in)和一个 r > R (out)。这导致两种不同的解决方案（仅在各自的域中有效），条件是 y_in(R) == y_out(R) 和 y_in'(R) == y_out'(R)。 Full 2-part solution where Radius = 1, a=99, b = 1 and constant = 1, y(inf) = constant

根据 Lutz Lehmann 的解决方案，它得到了正确的形状（至少对于内部区域而言，尽管不是在正确的尺度上）。

我只是不确定您将如何编写所有等式解决方案的代码，我想甚至首先获得他们的解决方案，尽管 Lutz 的回答在正确的方向上是一个了不起的观点。谢谢

Answer 1

问题

方程的阶数为2，因此状态向量的维数为2，取值总是y[0]，导数y[1]，没有y[2]，可能是 Matlab 翻译的残余。
同样在边界条件中，没有ya[2]，导数值为ya[1]，第二个中的函数值为yb[0]。
初始解猜测必须具有相同数量的 2 个状态分量。
为什么要用相同的数据计算两个解？
备注：不需要将return值转换为numpy类型，求解器会检查并转换。

具有奇点处理的 BVP 框架

ODE 在 r=0 处是奇异的，因此必须以特殊方式处理第一段。中值定理给出

(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²。所以如果你想要 y(r) 中的精度 1e-9，那么第一段不应该长于 1e-3.

与其尝试从 y(h) 和 y'(h) 的方程中消除 y0 以获得连接 ya[0] 和 ya[1] 的条件，不如让求解器也做这项工作并将 y0 作为参数添加到系统中。然后边界条件有3个slot对应为参数增加的虚拟维度，自然可以填入方程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 和生成的绘图

进行收敛求解器调用

SciPy: solve_bvp 问题二阶差分。方程式

SciPy: solve_bvp Problem 2nd Order Diff. Eq

python

scipy

differential-equations

问题

具有奇点处理的 BVP 框架