使用Scipy solve_bvp求解边值问题时遭遇数组维度错误求助

嘿，我帮你定位到了这个ValueError的问题所在——根源是你第一个问题里的fun函数返回数组时，最后一行的维度和前面的不匹配！

问题拆解 #

Scipy的solve_bvp对fun(x, y)的返回值有严格要求：它必须是一个形状为(n, m)的二维数组，其中：

• n是你定义的状态变量数量（这里是5个，对应y[0]到y[4]）
• m是输入x数组的元素个数（你的代码里x.size是400）

来看你第一个问题的fun函数：

def fun(x, y):
    tmp = np.vstack((
        np.cos(y[2]),
        np.sin(y[2]),
        y[3],
        (y[4]*np.cos(y[2])),
        0  # 这里就是bug的源头！
    ))
    return tmp

前四个元素都是长度为400的一维数组（因为y的每个元素都是(400,)的数组，运算后长度不变），当用np.vstack堆叠时，它们会被转换成(1, 400)的行向量。但最后一个元素0是个标量，np.vstack会把它转成(1, 1)的数组——这就导致堆叠时，前四行的第二维度是400，最后一行是1，维度不匹配，直接触发了那个错误。

而你的测试代码能正常运行，是因为测试里的fun返回的两个元素都是长度为5的一维数组（对应x.size=5），堆叠后维度完全一致，自然没问题。

修复方法 #

把fun里的标量0替换成和其他行长度一致的全0数组就行，最简单的方式是用np.zeros_like(y[2])（它会自动匹配y[2]的长度），或者np.zeros(x.size)：

def fun(x, y):
    tmp = np.vstack((
        np.cos(y[2]),
        np.sin(y[2]),
        y[3],
        (y[4]*np.cos(y[2])),
        np.zeros_like(y[2])  # 替换原来的标量0
    ))
    return tmp

这样所有行的长度都是400，堆叠后得到(5, 400)的二维数组，完美符合solve_bvp的要求。

小技巧：验证维度 #

你可以在fun里加一行打印来确认返回值的维度是否正确：

def fun(x, y):
    tmp = np.vstack((
        np.cos(y[2]),
        np.sin(y[2]),
        y[3],
        (y[4]*np.cos(y[2])),
        np.zeros_like(y[2])
    ))
    print(tmp.shape)  # 应该输出(5, 400)
    return tmp

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