使用scipy.solve_bvp求解ODE时遇IndexError问题求助

我已经花了几周时间尝试把一段常微分方程（ODE）转换成Python代码，但运行时不断出现各种错误。我自认对solve_bvp有基础理解，却始终找不到问题根源，反复阅读库文档也没完全搞明白。目前遇到的错误是：IndexError: too many indices for array。

原代码 #

import numpy as np
import scipy

#Constants#
P1 = 0.05
P2 = 0.1
Tau = 0.5
W = 1
Uh = 1
Uv = 1.562
Bv1 = 1.59
Bv2 = 1.59
A1 = 1.2
A2 = 0.8
Roe1 = 0.177
Roe2 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
Qv = 100
H = 1

#g is lamda this variable will change with time#
#This is variable test area#

g = 1

BH1 = 1
BH2 = 2

x1 = np.array([0, 1])
x2 = (0, 1)
t = np.linspace(0, 1, num=365)
#############################################

i = input('What is the value of i?')
i = ([int(i)])
#Functions#


def Bhinv(i):
    if i == 1:
        l = BH1*min(Nv/(Qv*int(Nh)), 1)+int(Roe1)*H*min(Nv/(Qh*int(Nh)), 1)
        return l
    if i == 2:
        l = BH2*min(Nv/(Qv*Nh), 1)+Roe2*H*min(Nv/(Qh*Nh), 1)
        return l
    else:
        return print("Null")


def sys1(x1, x2, i):
    if i.all() == 1:
        dx1dt = ((int(P1)/g)-1)*Uh*x1[0] + \
            (((int(Bhinv(1)))/g)*(int(Nh)/Nv)*x2[0])
        dx2dt = ((int(Bv1)/g)*(Nv/Nh)*x1[0])-int(Uv)*x2[0]
        return dx1dt, dx2dt

    if i.all() == 2:
        dx1dt = ((P2/g)-1)*Uh*x1(t)+(((BH2)/g)*(Nh/Nv)*x2(t))
        dx2dt = ((Bv2/g)*(Nv/Nh)*x1(t))-Uv*x2(t)
        return dx1dt, dx1dt

# BRN(Bhi) (i = 1,2)


def step1(x1, x2):
    x1 = 1
    x2 = 0
    return np.array([1, 0, 0])


y = ([0, 0],
     [0, 0])
v1s1 = scipy.integrate.solve_bvp(sys1, step1, x1, y, x2, y)

错误分析与修正方案 #

核心错误点 #

1. solve_bvp参数格式不匹配

   • 微分方程函数必须遵循fun(x, y)格式，原代码sys1额外传入i，且调用solve_bvp时传递了多余参数，导致参数解析混乱。
   • 边界条件函数要求为bc(ya, yb)，计算边界残差（即ya和yb需满足的等式，残差为0时满足条件），原step1逻辑完全错误，且返回值长度与方程数量不匹配。

2. 变量命名冲突与类型错误

   • x1被同时用作节点数组和微分方程的状态变量，导致索引错误；i被转为列表，用i.all()判断完全无效（列表无all()方法）。
   • y是嵌套列表而非numpy数组，solve_bvp要求状态数组为二维（形状(n, m)，n为方程数，m为节点数）。

3. 函数调用错误

   • sys1中i=2分支里的x1(t)、x2(t)属于错误调用，x1是数组不是函数，应改为索引访问。
   • Bhinv中不必要的int()转换会丢失精度（比如Nh=9.3转成9），影响计算结果。

修正后的代码 #

import numpy as np
from scipy.integrate import solve_bvp

# 常量定义
P1 = 0.05
P2 = 0.1
Uh = 1
Uv = 1.562
Bv1 = 1.59
Bv2 = 1.59
Roe1 = 0.177
Roe2 = 0.177
Nh = 9.3
Nv = 128
Qh = 10
Qv = 100
H = 1
g = 1
BH1 = 1
BH2 = 2

# 获取参数i，转为整数
i = int(input('What is the value of i?'))

def Bhinv(i_val):
    if i_val == 1:
        term1 = BH1 * min(Nv/(Qv * Nh), 1)
        term2 = Roe1 * H * min(Nv/(Qh * Nh), 1)
        return term1 + term2
    elif i_val == 2:
        term1 = BH2 * min(Nv/(Qv * Nh), 1)
        term2 = Roe2 * H * min(Nv/(Qh * Nh), 1)
        return term1 + term2
    else:
        raise ValueError("i must be 1 or 2")

# 定义微分方程函数，符合solve_bvp要求的格式：fun(x, y)
def ode_sys(x, y):
    # y[0]对应x1，y[1]对应x2
    x1 = y[0]
    x2 = y[1]
    if i == 1:
        b_hinv = Bhinv(1)
        dx1dt = ((P1/g) - 1) * Uh * x1
        dx1dt += (b_hinv / g) * (Nh / Nv) * x2
        dx2dt = (Bv1 / g) * (Nv / Nh) * x1 - Uv * x2
    elif i == 2:
        dx1dt = ((P2/g) - 1) * Uh * x1
        dx1dt += (BH2 / g) * (Nh / Nv) * x2
        dx2dt = (Bv2 / g) * (Nv / Nh) * x1 - Uv * x2
    else:
        raise ValueError("i must be 1 or 2")
    # 返回二维数组，形状(2, len(x))
    return np.vstack((dx1dt, dx2dt))

# 定义边界条件函数：bc(ya, yb)，返回残差数组
def boundary_cond(ya, yb):
    # 示例边界条件：初始点x=0时x1=1，x2=0；你需要根据实际模型修改这里
    res1 = ya[0] - 1  # x=0时x1=1
    res2 = ya[1] - 0  # x=0时x2=0
    # 如果有终点边界条件，添加到这里，比如res3 = yb[0] - ...
    return np.array([res1, res2])

# 设置节点数组（x轴，对应原代码的t）
x_nodes = np.linspace(0, 1, num=365)
# 设置初始猜测的状态数组，形状(2, 365)
y_guess = np.zeros((2, len(x_nodes)))

# 调用solve_bvp
solution = solve_bvp(ode_sys, boundary_cond, x_nodes, y_guess)

# 检查求解结果
if solution.success:
    print("求解成功")
    # 可以获取解：solution.x是节点，solution.y是对应的状态值
else:
    print(f"求解失败，原因：{solution.message}")

关键修改说明 #

• 重命名了冲突变量，确保x_nodes作为求解节点数组，状态变量用y[0]、y[1]表示，避免命名冲突。
• 调整了微分方程函数和边界条件函数的格式，完全符合solve_bvp的要求。
• 修正了i的类型，直接用整数判断，去掉无效的列表转换。
• 移除了不必要的int()转换，保留计算精度。
• 添加了求解结果的检查逻辑，方便排查问题。

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