求解两个耦合的二阶边值问题

Question

我已经使用模块solve_bvp 求解了具有两个边界条件的单个二阶微分方程。但是，现在我正在尝试求解两个二阶微分方程组；

U'' + a*B' = 0

B'' + b*U' = 0

边界条件 U(+/-0.5) = +/-0.01 和 B(+/-0.5) = 0。我已将其拆分为一阶常微分方程系统，并尝试使用 @987654322 @ 以数字方式解决它们。但是，对于我的解决方案，我只是得到了充满零的数组。我相信我错误地执行了边界条件。我不清楚如何处理文档中的两个以上方程。我的尝试如下

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import solve_bvp
alpha = 1E-8
zeta = 8E-3
C_k = 0.05
sigma = 0.01
def fun(x, y):

    return np.vstack((y[1],-((alpha)/(C_k*sigma))*y[2],y[2], -(1/(C_k*zeta))*y[1]))

def bc(ya, yb):

    return np.array([ya[0]+0.001, yb[0]-0.001,ya[0]-0, yb[0]-0])

x = np.linspace(-0.5, 0.5, 5000)
y = np.zeros((4, x.size))
print(y)


sol = solve_bvp(fun, bc, x, y)
print(sol)

在我的问题中，我刚刚重新标记了 a 和 b，但它们只是我输入的参数。我有这组方程的解析解，所以我知道存在一个不平凡的方程。任何帮助将不胜感激。

Answer 1

如果您在评论中至少声明一次或通过分配给特定命名的变量您希望如何组成状态向量，大多数时候真的很有帮助。

通过导数返回向量的形式，我想你打算

U, U',  B, B'

这意味着U=y[0], U'=y[1] 和B=y[2],B'=y[3]，这样你的导数向量应该是正确的

return y[1], -((alpha)/(C_k*sigma))*y[3],  y[3], -(1/(C_k*zeta))*y[1]

和边界条件

return ya[0]+0.001, yb[0]-0.001, ya[2]-0, yb[2]-0

特别是你的边界条件应该在第一步中抛出算法，因为一个奇异的雅可比行列式，总是检查解决方案结构的.success 字段和.message 字段。

---

注意，默认实验solve_bvp的绝对和相对容差为1e-3，节点数量限制为500。

将初始节点数设置为 50（5000 太多了，求解器会在必要时进行细化），并将容差设置为 1-6，我得到以下解图，明显满足边界条件。