基于scipy.integrate.solve_bvp求解7组耦合二阶微分边值问题

1. 自定义fun(x, y)的正确写法 #

solve_bvp要求fun(x, y)返回每个状态变量的一阶导数，其中y是形状为(14, N)的数组（N是初始网格点数量）：

• 前7行对应y1(x)到y7(x)
• 后7行对应z1(x)=dy1/dx到z7(x)=dy7/dx

根据你转化的一阶ODE，导数的计算逻辑如下：

• 对于y1到y7的导数：直接取对应的z分量，即y[7:14, :]
• 对于z1到z7的导数：按公式-(1/x)*z_i - L_i(y1,...,y7)计算，其中L_i是y1到y7的线性组合

示例代码（假设L是7×7的系数矩阵，L[i]对应L_i的系数）：

import numpy as np
from scipy.integrate import solve_bvp

def fun(x, y):
    # y的形状：(14, N)
    dydx = np.zeros_like(y)
    # 前7个方程：dyi/dx = zi
    dydx[:7, :] = y[7:, :]
    # 后7个方程：dzi/dx = -(1/x)*zi - Li(y1..y7)
    # 计算Li：用矩阵乘法，L是7x7系数矩阵，y[:7,:]是(7,N)
    L_terms = L @ y[:7, :]
    dydx[7:, :] = -(1/x) * y[7:, :] - L_terms
    return dydx

2. 边界条件残差函数bc(ya, yb)的定义 #

• ya是x=a处所有状态变量的值，形状为(14,)；yb是x=b处所有状态变量的值，形状为(14,)
• 残差函数需要返回一个长度为14的数组，每个元素对应一个边界条件的残差（残差=实际值-目标值，solve_bvp会将残差收敛到0）

根据你的边界条件：

• x=a处：z_i(a)=A_i → 残差为ya[7:] - A（A是包含A1到A7的长度7数组）
• x=b处：y_i(b)=B_i → 残差为yb[:7] - B（B是包含B1到B7的长度7数组）

示例代码：

def bc(ya, yb):
    # A = [A1, A2, ..., A7]，B = [B1, B2, ..., B7]
    res_a = ya[7:] - A
    res_b = yb[:7] - B
    return np.concatenate([res_a, res_b])

3. 初始网格x的选择 #

• 初始网格是a到b之间的一维数组，无需太密集，solve_bvp会自适应调整网格点位置
• 最简单的方式是用均匀分布的点，比如：

a = 1.0  # 你的左边界
b = 10.0 # 你的右边界
x_init = np.linspace(a, b, 10)  # 10个初始点，数量可根据问题调整

• 如果已知解在某个区间变化剧烈，可以在该区域多放几个初始点，比如在x=a附近增加点密度。

4. 初始猜测y的形状与取值规则 #

• y的形状必须和fun的输入一致：(14, N)，其中N是初始网格x_init的点数
• 取值规则：尽量贴近真实解的趋势，不需要严格满足边界条件（solve_bvp会自动修正）：
  • 前7行（y1到y7）：可以用边界条件yi(b)=Bi构造，比如设为常数Bi，或线性变化（比如从某个初始值到Bi）
  • 后7行（z1到z7）：可以用边界条件zi(a)=Ai设为常数Ai，或根据y的猜测计算导数（比如假设yi线性变化，zi=(Bi - yi_a)/(b-a)，yi_a可暂时设为0或Bi）

示例初始猜测：

# 初始化y为(14, 10)的数组
y_init = np.zeros((14, len(x_init)))
# 前7行：设为对应的Bi（x=b处的边界值）
for i in range(7):
    y_init[i, :] = B[i]
# 后7行：设为对应的Ai（x=a处的边界值）
for i in range(7):
    y_init[7+i, :] = A[i]

内容的提问来源于stack exchange，提问作者plierp