使用Python求解四阶BVP时solve_bvp边界条件返回值形状错误

问题描述 #

我要解一个适用于软土地基的四阶常微分方程：
$$EI \cdot \frac{d^4 y}{dx^4} = -k \cdot y$$
将其转化为一阶方程组：

• $y_1' = y_2$
• $y_2' = y_3$
• $y_3' = y_4$
• $y_4' = -\frac{k \cdot y_1}{EI}$

对应的边界条件：

• $y_3(0) = y_3(L) = 0$
• $y_4(0) = -\frac{F}{EI}$
• $y_4(L) = 0$

我的代码 #

from scipy.integrate import solve_bvp
import numpy as np

F = 300 #kN
EI = 20000 #kNm2
D = 1 #m
M_E = 45000 #kN/m2
k = 1.4*M_E/D #kN/m3
L = 10 #m

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

p = np.array([k,EI,F])

print(p)

def fun(x, y, p):
    k = p[0]
    EI = p[1]
    return np.vstack((y[1],y[2],y[3],-k*y[0]/EI))

def bc(ya, yb, p):
    F = p[2]
    EI = p[1]
    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, p)

遇到的错误 #

ValueError: 'bc' return is expected to have shape (7,), but actually has (4,)

问题原因与解决办法 #

错误原因 #

solve_bvp的p参数是用来反求未知参数的，不是传递已知常量的工具。当你传入p时，函数要求边界条件函数bc的返回值长度等于「状态变量数量 + 参数数量」：

• 这里状态变量是4个（y1到y4）
• 你传入的p包含3个元素
• 因此要求返回长度为4+3=7，但你的bc只返回了4个边界条件，触发了报错。

解决办法 #

既然k、EI、F都是你预先定义好的已知量，完全不需要放到p里。直接在fun和bc里调用这些全局变量即可，修正后的代码如下：

from scipy.integrate import solve_bvp
import numpy as np

F = 300 #kN
EI = 20000 #kNm2
D = 1 #m
M_E = 45000 #kN/m2
k = 1.4*M_E/D #kN/m3
L = 10 #m

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

# 移除p参数，直接使用全局定义的常量
def fun(x, y):
    return np.vstack((y[1], y[2], y[3], -k*y[0]/EI))

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)

# 验证求解结果并可视化
print(res_a.success)
if res_a.success:
    y_sol = res_a.sol(x)[0]
    import matplotlib.pyplot as plt
    plt.plot(x, y_sol)
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.title('Displacement on Soft Soil Foundation')
    plt.show()

如果你的场景确实需要反求p中的参数，那你需要补充额外3个边界条件（总共7个），让bc返回长度为7的数组，这样solve_bvp才能同时求解状态变量和未知参数。但从你的问题描述来看，这些参数都是已知量，用上面的方法就能解决问题。

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