Python求解器无返回值：梁的横竖固有频率参数求解故障

问题概述 #

需要计算梁的竖直固有频率和水平固有频率，使其尽可能逼近给定的目标值。计算依赖y、z方向的转动惯量与表面积，待优化参数为宽度B、高度H以及厚度t_h、t_v。

分两步使用SymPy求解：

• 第一步求解B和H，运行正常，能输出有效结果；
• 第二步求解t_h和t_v时，求解器始终返回空列表[]，尝试linsolve()、拆分方程等方法均无效（拆分方程缺少转动惯量/表面积初始值）。

原代码 #

import numpy as np
from sympy import symbols, Eq, solve

w_2_vb_g=4.341
w_2_hb_g=2.601
E=210E9
rho=7700
L=100
kL_2=7.85
kL_3=11.00

B,H,t_h,t_v=symbols('B H t_h t_v', real=True, positive=True)
w_2_vb_rec=Eq(((((kL_2)**2)/(L**2))*(((E*((1/12)*H*B**3))/(rho*(B*H)))**(1/2))), w_2_vb_g)
w_2_hb_rec=Eq(((((kL_2)**2)/(L**2))*(((E*((1/12)*B*H**3))/(rho*(H*B)))**(1/2))), w_2_hb_g)
Solution=solve((w_2_hb_rec,w_2_vb_rec), (B,H))
print(Solution)

print("Solutions:")
for var, val in Solution.items():
    print(f"{var}: {val}")

list_from_dict = list(Solution.items())
print("List from dict:", list_from_dict)

B=round(float(list_from_dict[1][1]*1.03),3)
H=round(float(list_from_dict[0][1]*1.03),3)

B_list=np.array([t_h, (H-2*t_h), B, t_v])
H_list=np.array([B, t_v, t_h, (H-2*t_h)])
D0_list=np.array([0, (0.5*B-0.5*t_v), (0.5*H-0.5*t_h),0])

w_2_vb=Eq(((((kL_2)**2)/(L**2))*(((E*(2*(((1/12)*B_list[0]*H_list[0]**3)+B_list[0]*H_list[0]*D0_list[0]**2)+2*(((1/12)*B_list[1]*H_list[1]**3)+B_list[1]*H_list[1]*D0_list[1]**2))/(rho*(((2*(B_list[0]*H_list[0])))+(2*(B_list[1]*H_list[1]))))**(1/2))))), w_2_vb_g)
w_2_hb=Eq(((((kL_2)**2)/(L**2))*(((E*(2*(((1/12)*B_list[2]*H_list[2]**3)+B_list[2]*H_list[2]*D0_list[2]**2)+2*(((1/12)*B_list[3]*H_list[3]**3)+B_list[3]*H_list[3]*D0_list[3]**2))/(rho*(((2*(H_list[2]*B_list[2])))+(2*(H_list[3]*H_list[3]))))**(1/2))))), w_2_hb_g)
Solution_beam=solve((w_2_hb,w_2_vb), (t_h,t_v))
print(Solution_beam)

控制台输出 #

{H: 0.279980237079456, B: 0.467279588297547}
Solutions:
H: 0.279980237079456
B: 0.467279588297547
List from dict: [(H, 0.279980237079456), (B, 0.467279588297547)]
[]

问题根源分析 #

1. 符号变量被覆盖：第一步求解得到符号变量B、H的解析解后，直接将其赋值为浮点数（B=round(float(...),3)），导致后续构建t_h、t_v的方程时，B、H已变为数值而非符号变量，破坏了符号方程的合法性。
2. 方程笔误：w_2_hb方程的分母部分错误使用H_list[3]*H_list[3]计算表面积，正确应为B_list[3]*H_list[3]，导致方程物理模型错误，无法得到有效解。
3. 符号求解局限性：第二步的方程包含t_h、t_v的高次项与交叉项，属于非线性方程组，SymPy的符号求解器难以处理这类复杂方程，无法找到解析解。

修正方案 #

1. 避免符号变量覆盖 #

将第一步得到的B、H数值用新变量存储，保留原符号变量的作用域。

2. 修正方程笔误 #

修正表面积计算的错误项，确保物理模型符合实际。

3. 改用数值求解 #

由于方程是非线性高次方程组，改用数值优化库（如scipy.optimize.root）进行求解。

修正后代码示例 #

import numpy as np
from sympy import symbols, Eq, solve
from scipy.optimize import root

w_2_vb_g=4.341
w_2_hb_g=2.601
E=210E9
rho=7700
L=100
kL_2=7.85
kL_3=11.00

# 第一步：求解B和H的符号解
B_sym, H_sym, t_h, t_v = symbols('B H t_h t_v', real=True, positive=True)
w_2_vb_rec = Eq((kL_2**2 / L**2) * np.sqrt((E * (1/12 * H_sym * B_sym**3)) / (rho * B_sym * H_sym)), w_2_vb_g)
w_2_hb_rec = Eq((kL_2**2 / L**2) * np.sqrt((E * (1/12 * B_sym * H_sym**3)) / (rho * H_sym * B_sym)), w_2_hb_g)
solution_bh = solve((w_2_hb_rec, w_2_vb_rec), (B_sym, H_sym))
print("第一步求解结果：", solution_bh)

# 提取B和H的数值，避免覆盖符号变量
B_val = round(float(solution_bh[B_sym] * 1.03), 3)
H_val = round(float(solution_bh[H_sym] * 1.03), 3)
print(f"调整后的B: {B_val}, H: {H_val}")

# 定义目标函数，用于数值求解
def equations(x):
    t_h, t_v = x
    # 构建各参数列表，使用数值B_val、H_val
    B_list = np.array([t_h, (H_val - 2*t_h), B_val, t_v])
    H_list = np.array([B_val, t_v, t_h, (H_val - 2*t_h)])
    D0_list = np.array([0, (0.5*B_val - 0.5*t_v), (0.5*H_val - 0.5*t_h), 0])
    
    # 计算竖直固有频率残差
    I_v = 2 * ((1/12 * B_list[0] * H_list[0]**3) + B_list[0] * H_list[0] * D0_list[0]**2) + \
          2 * ((1/12 * B_list[1] * H_list[1]**3) + B_list[1] * H_list[1] * D0_list[1]**2)
    A_v = 2*(B_list[0]*H_list[0]) + 2*(B_list[1]*H_list[1])
    w_v = (kL_2**2 / L**2) * np.sqrt((E*I_v)/(rho*A_v)) - w_2_vb_g
    
    # 计算水平固有频率残差（修正表面积计算笔误）
    I_h = 2 * ((1/12 * B_list[2] * H_list[2]**3) + B_list[2] * H_list[2] * D0_list[2]**2) + \
          2 * ((1/12 * B_list[3] * H_list[3]**3) + B_list[3] * H_list[3] * D0_list[3]**2)
    A_h = 2*(H_list[2]*B_list[2]) + 2*(B_list[3]*H_list[3])
    w_h = (kL_2**2 / L**2) * np.sqrt((E*I_h)/(rho*A_h)) - w_2_hb_g
    
    return [w_v, w_h]

# 设置初始猜测值（可根据实际情况调整）
initial_guess = [0.02, 0.02]
# 调用数值求解器
result = root(equations, initial_guess)

if result.success:
    t_h_sol, t_v_sol = result.x
    print(f"求解得到t_h: {round(t_h_sol, 5)}, t_v: {round(t_v_sol, 5)}")
else:
    print("求解失败，原因：", result.message)

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