用solve_bvp求解耦合微分方程组

Question

我想要解一个由七个耦合二阶微分方程组成的边值问题。有7个函数，y1(x),...y7(x)，每个函数都由一个形式的微分方程来描述。

d^2yi/dx^2 = -(1/x)*dyi/dx - Li(y1,...,y7) for 0 < a <= x <= b,

其中Li是一个给出y1,...,y7线性组合的函数。我们给出了一阶导数dyi/dx在x=a和函数yi在x=b的边界条件。

dyi/dx(a) = Ai,
yi(b) = Bi.

因此，我们可以把它重写为一个由14个耦合的一阶ODE组成的系统：

dyi/dx = zi,
dzi/dx = -(1/x)*zi - Li(y1,...,y7),

zi(a) = Ai,
yi(b) = Bi.

我想用scipy.integrate.solve_bvp函数来求解这个方程组。但是，我很难理解文档(https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html)中描述的函数的确切输入参数应该是什么。

该函数需要的第一个参数是可调用的fun(x,y)。据我理解，输入参数y必须是一个由yi和zi值组成的数组，并给出zi和dzi/dx的值作为输出。所以我的函数看起来如下(伪代码)：

def fun(x,y):
    y1, z1, y2, z2, ..., y7, z7 = y
    return [z1, -(1/x)*z1 - L1(y1,...,y7),
            ...,
            z7, -(1/x)*z7 - L7(y1,...,y7)]

对吗？

然后，solve_bvp的第二个参数是一个可调用的bc(ya,yb)，它应该计算边界条件的剩余值。在这里，我真的很难理解如何定义这样的函数。我也不清楚ya和yb究竟是什么数组，它们应该具有什么形状？

第三个参数是x，它是带有形状(m,)的“初始网格”。x应该只包含点a和b，这是我们知道边界条件的地方吗？还是应该是别的什么？

最后，第四个参数是y，它是“网格节点上函数值的初始猜测”，具有形状(n,m)。它的第1列对应于x[i]。我想第一行对应于y1，第二行对应于z1，第三行对应于y2等等，对吗？此外，应该把哪些值放在这里？我们可以在x=a和x=b上输入已知的边界条件，但是我们不知道函数在任何其他点上的样子。此外，这个y与bc(ya,yb)函数的关系如何？输入参数ya,yb是从这个y派生出来的吗？

如果能帮助您理解solve_bvp的语法及其在本例中的应用，我们将不胜感激。

Answer 1

对于边界条件，bc返回方程的残差。也就是说，转换方程，使右侧为零，然后返回左侧的向量。ya,yb是x=a,b点处的状态向量。

对于最初的状态，它可能是任何东西。然而，在大多数情况下，“任何”都会由于奇异雅可比或过大的网格而导致收敛失败。求解者所做的就是求解一个大型非线性方程组(参见配置法)。和任何这样的系统一样，如果你开始的时候离实际的解足够近，收敛是最有保证的。所以你可以试试yi

具有值的Bi

• linear函数与斜率Ai的
• 常数函数和x=b

上的值Bi

用这些函数的导数来表示zi的值。这不能保证这是可行的，特别是当接近于零的时候，你有接近常量解和对数的基解，a越接近于零，它就变得越奇异。