使用scipy solve_bvp求解四阶ODE系统时遇IndexError问题求助

你遇到的IndexError是因为定义的状态变量数组y只有6个元素（索引范围0-5），但代码中错误地使用了y[6]，超出了数组的有效索引范围。

你的ODE系统包含：一个二阶方程（对应2个状态变量）和一个四阶方程（对应4个状态变量），总共需要6个状态变量。需要重新梳理状态变量的对应关系，修正所有索引错误：

1. 状态变量映射：

   • y[0] = f，y[1] = f'（对应二阶ODE的两个状态）
   • y[2] = g，y[3] = g'，y[4] = g''，y[5] = g'''（对应四阶ODE的四个状态）

2. 修正后的代码：

from scipy.integrate import solve_bvp
import numpy as np


def ode(r, y):
    Re = 100
    # 状态变量：y[0]=f, y[1]=f'; y[2]=g, y[3]=g', y[4]=g'', y[5]=g'''
    return np.vstack((
        y[1],  # f' = y[1]
        2 / r**y[0] - (2/3)*Re / r**2 * y[2] * y[0],  # f''的表达式
        y[3],  # g' = y[3]
        y[4],  # g'' = y[4]
        y[5],  # g''' = y[5]
        12 / r**2 * y[3] - 24 / r**3 * y[3] + (2*Re)/(5*r**3)*(
            5*y[0]*(r*y[1] - 2*y[0]) +
            2*r*y[2]*y[5] +
            (r*y[3] - 4*y[2])*y[4] -
            18/r*y[3]*y[4] +
            48/r**2*y[2]**2
        )
    ))


def bc(ya, yb):
    # 对应各状态变量的边界条件
    return np.array([
        ya[0],        # f在r=0.7处为0
        yb[0] - 1,    # f在r=1处为1
        ya[1],        # f'在r=0.7处为0
        yb[1],        # f'在r=1处为0
        ya[2],        # g在r=0.7处为0
        yb[2],        # g在r=1处为0
        ya[3],        # g'在r=0.7处为0
        yb[3],        # g'在r=1处为0
        ya[4],        # g''在r=0.7处为0
        yb[4],        # g''在r=1处为0
        ya[5],        # g'''在r=0.7处为0
        yb[5]         # g'''在r=1处为0
    ])

x_flow = np.linspace(0.7, 1, 10)
y_flow = np.ones((6, x_flow.shape[0]))

res_flow = solve_bvp(ode, bc, x_flow, y_flow, verbose=2)

• 替换了所有越界的y[6]引用为正确的y[5]
• 调整了g相关变量的索引（比如原代码中y[3]对应g，修正为y[2]）
• 同步修正了边界条件函数bc中的所有索引，确保与状态变量映射一致

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