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

Question

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

U'' + a*B' = 0

B'' + b*U' = 0

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

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，但它们只是我输入的参数。我有这组方程的解析解，所以我知道存在一个不平凡的方程。任何帮助将不胜感激。

标签： pythonnumpyodedifferential-equations

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，我得到以下明显满足边界条件的解图。