使用scipy solve_bvp求解非齐次四阶BVP时遇维度不匹配错误

错误原因 #

scipy.integrate.solve_bvp 是自适应求解器，它不会严格固定你初始传入的100个节点，而是会根据求解过程动态调整节点数量（比如合并或拆分节点）。你预先定义的 vector 是固定长度100的数组，但求解器调用fun时传入的y[0]长度是它当前使用的节点数（这里是99），直接导致维度不匹配报错。

齐次情况能正常运行，是因为没有用到固定长度的外部数组，所有运算都基于求解器传入的y和x，维度自然匹配。

修复方法 #

不要预先定义固定长度的非齐次项数组，而是在fun函数内部根据传入的x参数生成对应长度的数组。如果你的K是和初始x对应的自定义数据，需要用插值方法（比如np.interp）根据求解器当前的x节点获取对应的值。

针对全1非齐次项的修复代码 #

from scipy.integrate import solve_bvp
import numpy as np
    
L = 10 
nodes = 100

A = 1000
B = 1500
C = 0.05

x = np.linspace(0,L,nodes)

def fun(x, y):
    # 生成与当前求解器节点x长度一致的非齐次项
    vector = np.ones_like(x)
    ans = vector - np.multiply(C, y[0])
    return np.vstack((y[1], y[2], y[3], ans))

def bc(ya, yb):
    return np.array([ya[2], yb[2], ya[3]+A/B, yb[3]])

y_a = np.zeros((4, x.size))

res_a = solve_bvp(fun, bc, x, y_a)

res1 = res_a.sol(x)[0]
res2 = res_a.sol(x)[1]
res3 = B*res_a.sol(x)[2]
res4 = B*res_a.sol(x)[3]

自定义K向量的扩展处理 #

如果你的K是和初始x对应的自定义向量（比如随机数组），需要通过插值匹配求解器的动态节点：

from scipy.integrate import solve_bvp
import numpy as np
    
L = 10 
nodes = 100

A = 1000
B = 1500
C = 0.05

x_initial = np.linspace(0,L,nodes)
# 自定义初始x对应的K向量
K = np.random.rand(nodes)

def fun(x, y):
    # 插值获取当前求解器节点对应的K值
    current_K = np.interp(x, x_initial, K)
    ans = current_K - np.multiply(C, y[0])
    return np.vstack((y[1], y[2], y[3], ans))

def bc(ya, yb):
    return np.array([ya[2], yb[2], ya[3]+A/B, yb[3]])

y_a = np.zeros((4, x_initial.size))

res_a = solve_bvp(fun, bc, x_initial, y_a)

res1 = res_a.sol(x_initial)[0]
res2 = res_a.sol(x_initial)[1]
res3 = B*res_a.sol(x_initial)[2]
res4 = B*res_a.sol(x_initial)[3]

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