Mathematica转Python：振动系统特征向量求解得全零解，与原结果不符

问题背景 #

需将Mathematica振动系统求解代码转为Python，原Mathematica代码求解广义特征值问题：

Solve[(K - w*M) . a == 0]
其中K、M为7阶方阵，a = {x1, x2, x3, x4, x5, x6, x7}为位移向量，w对应系统固有频率，每个w对应一组x_i的相对位移关系。

用户现有Python代码 #

from sympy.solvers import solve
from sympy import *
import numpy as np

m=3.0
k=1.50

w,x1,x2,x3,x4,x5,x6,x7 = symbols("w,x1,x2,x3,x4,x5,x6,x7")

M = Matrix([[m,0,0,0,0,0,0], [0,m,0,0,0,0,0],[0,0,m,0,0,0,0],[0,0,0,m,0,0,0],[0,0,0,0,m,0,0],[0,0,0,0,0,m,0],[0,0,0,0,0,0,m]])
K = Matrix([[2*k,-k,0,0,0,0,0], [-k,2*k,-k,0,0,0,0],[0,-k,2*k,-k,0,0,0],[0,0,-k,2*k,-k,0,0],[0,0,0,-k,2*k,-k,0],[0,0,0,0,-k,2*k,-k],[0,0,0,0,0,-k,2*k]])
xn = Matrix([x1,x2,x3,x4,x5,x6,x7])

D1=K-w*M

print("omega squared")
A=solve(D1.det(), w)
print(A)

res =(np.array(A))**0.5

for i in range(7):
  omegan=float(res[i])
  wn=round(omegan**2,5)
  D1=K-wn*M
  print(D1*xn)
  print("\n** Positions x_i ", i+1, "for omega= ", np.sqrt(wn), " son; ",solve(D1*xn,xn))

现存问题 #

• 部分K、M组合下计算耗时过长甚至无响应
• 求解x_i相对位移时始终得到全零解{x1: 0.0, x2: 0.0, ..., x7: 0.0}，无法得到Mathematica输出的非零相对位移关系

解决办法 #

核心问题解析 #

原代码的低效与全零解问题，本质是用了错误的求解路径：

1. 手动计算行列式再解方程的方式，对高阶矩阵符号计算效率极低
2. solve(D1*xn, xn)会返回齐次方程的所有解，全零解是合法的平凡解，但我们需要的是非零的特征向量（振型）

方案1：SymPy广义特征值直接求解（符号/半符号） #

利用SymPy内置的广义特征值求解方法，直接处理K*a = w²*M*a问题，同时得到特征值（固有频率平方）与特征向量（振型）：

from sympy import Matrix, N

m = 3.0
k = 1.50

# 定义质量矩阵与刚度矩阵
M = Matrix([[m,0,0,0,0,0,0], 
            [0,m,0,0,0,0,0],
            [0,0,m,0,0,0,0],
            [0,0,0,m,0,0,0],
            [0,0,0,0,m,0,0],
            [0,0,0,0,0,m,0],
            [0,0,0,0,0,0,m]])
K = Matrix([[2*k,-k,0,0,0,0,0], 
            [-k,2*k,-k,0,0,0,0],
            [0,-k,2*k,-k,0,0,0],
            [0,0,-k,2*k,-k,0,0],
            [0,0,0,-k,2*k,-k,0],
            [0,0,0,0,-k,2*k,-k],
            [0,0,0,0,0,-k,2*k]])

# 求解广义特征值问题：K*v = w_sq*M*v
eigen_pairs = K.eigenvects(M)

print("固有频率及对应振型：")
for idx, (w_sq, _, vecs) in enumerate(eigen_pairs, 1):
    # 计算固有频率
    omega = N(w_sq)**0.5
    print(f"\n第{idx}阶固有频率：{N(omega)}")
    # 取特征向量并归一化，得到相对位移关系
    mode = vecs[0]
    # 用第一个非零元素归一化
    norm_mode = mode / mode[mode.nonzero()[0][0]]
    for i, xi in enumerate(norm_mode, 1):
        print(f"x{i} = {N(xi)}")

方案2：SciPy数值广义特征值求解（高效，适合高阶矩阵） #

若矩阵规模较大，符号计算效率不足，推荐用SciPy的对称矩阵广义特征值专用函数，数值计算效率极高：

import numpy as np
from scipy.linalg import eigh

m = 3.0
k = 1.50

# 转为numpy数组
M = np.diag([m]*7)
K = np.array([[2*k,-k,0,0,0,0,0], 
              [-k,2*k,-k,0,0,0,0],
              [0,-k,2*k,-k,0,0,0],
              [0,0,-k,2*k,-k,0,0],
              [0,0,0,-k,2*k,-k,0],
              [0,0,0,0,-k,2*k,-k],
              [0,0,0,0,0,-k,2*k]])

# 求解广义特征值问题，得到特征值（w²）与特征向量（振型）
w_sq, modes = eigh(K, M)

print("固有频率及对应振型：")
for idx in range(7):
    omega = np.sqrt(w_sq[idx])
    print(f"\n第{idx+1}阶固有频率：{omega}")
    # 取对应振型列向量，归一化
    mode = modes[:, idx]
    norm_mode = mode / mode[0] if mode[0] != 0 else mode / mode[np.nonzero(mode)[0][0]]
    print("相对位移：", norm_mode)

原代码修复（仅解决全零解问题） #

若需保留原代码框架，可将solve(D1*xn,xn)替换为求解矩阵零空间的方法，获取非零基础解系：

# 原循环内修改
null_space = D1.nullspace()
if null_space:
    mode = null_space[0]
    norm_mode = mode / mode[mode.nonzero()[0][0]]
    print("相对位移关系：")
    for i, xi in enumerate(norm_mode, 1):
        print(f"x{i} = {N(xi)}")

内容的提问来源于stack exchange，提问作者Elizabeth Hernández Marín