Scipy solve_bvp迭代时数组尺寸变化报错的解决求助

问题描述 #

我正在求解弹性地基梁的四阶ODE，用于模拟埋设在多层土层中的横向受荷桩，土层参数可在代码中指定。数组c存储系统刚度参数，需要与y[0]逐元素相乘，但迭代过程中y[0]的长度减少了1个元素，出现如下错误：

ValueError: operands could not be broadcast together with shapes (40,) (39,)

我尝试在fun(x,y)的return语句前添加以下代码，但未解决问题：

if c.size < y[0].size:
        np.delete(c,-1)

完整代码 #

# ODE-Solver ---------------------------------------------------

from scipy.integrate import solve_bvp
import numpy as np


#Input ---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*

layers = [2,2,6] #m
M_Es = (10**3)*[45,20,60] #MN/m2

E_c = 33.6*10**6

F = 1000 #kN
D = 1 #m
nodes = 40

#---------------------------------------------------

#total length of the pile:
L = sum(layers) #m

#x-coordinates of the pile
x = np.linspace(0,L,nodes)

#Empty arr. for soil's compression modulus:
M_E = np.zeros(nodes)

#soil layer boundaries' depths:
depth = np.zeros(len(layers)+1)

for i in range(1,len(depth)):
    depth[i] = sum(layers[:i])


i=0

#assign M_E to the nodes of the current soil layer "j"
#while x lies between two boundaries assign the i-th value in "M_Es" array
for j in range(len(depth)-1):
    while x[i] >= sum_layers[j] and x[i] < sum_layers[j+1]:
        M_E[i] = M_Es[j]
        i += 1
    
#calculate the spring stiffness of the soil "k_sh" and the flexural rigidity of the pile "EI"
k_sh = 1.4*M_E/D #kN/m3
EI = E_c*0.25*np.pi*(D*0.5)**4 #kNm2


c = k_sh*D/(EI)

def fun(x, y):
    return np.vstack((y[1],y[2],y[3],np.multiply(-c,y[0])))

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

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

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

问题根源与解决方法 #

问题原因 #

solve_bvp在迭代过程中会自动调整节点数量（删减或新增节点），导致传入fun的x和y维度与初始定义的固定长度c不匹配。你手动删除c元素的方法无效，因为np.delete不会修改原数组（它返回新数组），且无法动态匹配每次迭代的节点数。另外原代码中sum_layers未定义，属于语法错误。

正确解决方案 #

不要预先计算固定长度的c，而是在fun内部根据当前传入的x值实时计算对应位置的刚度参数，确保参数长度始终和y[0]一致：

修正后的完整代码

# ODE-Solver ---------------------------------------------------

from scipy.integrate import solve_bvp
import numpy as np


#Input ---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*---*

layers = [2,2,6] #m
M_Es = (10**3)*np.array([45,20,60]) #MN/m2 转为numpy数组方便索引

E_c = 33.6*10**6

F = 1000 #kN
D = 1 #m
nodes = 40

#---------------------------------------------------

#total length of the pile:
L = sum(layers) #m

#x-coordinates of the pile
x = np.linspace(0,L,nodes)

#calculate the flexural rigidity of the pile "EI"
EI = E_c*0.25*np.pi*(D*0.5)**4 #kNm2

def fun(x, y):
    # 根据当前x的位置，实时计算每个点对应的M_E值
    M_E_current = np.zeros_like(x)
    depth = np.cumsum([0] + layers)  # 简洁计算土层深度边界
    for j in range(len(layers)):
        mask = (x >= depth[j]) & (x < depth[j+1])
        M_E_current[mask] = M_Es[j]
    
    # 实时计算当前x对应的刚度参数c
    k_sh_current = 1.4 * M_E_current / D
    c_current = k_sh_current * D / EI
    
    return np.vstack((y[1], y[2], y[3], -c_current * y[0]))

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

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

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

额外说明 #

• 使用np.cumsum替代手动循环求和，更简洁地生成土层深度边界
• 将M_Es转为numpy数组，确保索引和赋值操作更稳定
• 实时计算c_current完全适配solve_bvp的节点调整机制，从根源解决维度不匹配问题

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