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

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

Answer 1

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

根据导数return向量的形式，我认为你打算

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

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

Solving two coupled second order boundary value problems

python

numpy

ode

differential-equations