使用Numpy求解3自由度弹簧阻尼系统特征值时遇类型错误与复数警告

问题背景 #

使用Numpy的eig、inv、transpose、arctan2等函数求解3自由度弹簧阻尼系统，修改M、C、K矩阵及F0向量后出现两个问题：

1. TypeError错误 #

phi = arctan2(-2*zeta*wn, wn**2-w**2)
TypeError: ufunc 'arctan2' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

2. ComplexWarning警告 #

ComplexWarning: Casting complex values to real discards the imaginary part
A[:, n] = b*X

问题根源 #

1. TypeError：新K矩阵结构导致inv(M)@K的特征值D出现复数，后续wn = sqrt(D)得到复数数组，而arctan2仅支持实数输入，触发类型错误。
2. ComplexWarning：特征向量V为复数，但数组A被定义为浮点类型，赋值时虚部被强制丢弃，导致模态振型信息丢失。

修复方案 #

1. 确保刚度矩阵正定（优先方案） #

常规振动系统的刚度矩阵K必须正定，新K矩阵中存在k3-k1、k4-k2这类可能为负的对角元，极易导致矩阵非正定，产生复特征值。需检查物理参数，保证K矩阵所有顺序主子式为正，确保特征值为实数。

2. 处理复数特征值与向量（兼容非正定K） #

若物理模型确实需要非正定K，需对复数运算做如下处理：

• 提取特征值的实部，过滤负值避免平方根出错；
• 将模态振型数组A改为复数类型，保留完整的特征向量信息；
• 所有后续计算中，对复数结果取实部用于最终的振动响应计算。

3. 修正模态阻尼计算错误 #

原代码中transpose(A)@C*A是逐元素相乘，正确的模态阻尼矩阵计算应使用矩阵乘法transpose(A)@C@A，否则模态阻尼系数zeta计算错误。

修改后的核心代码 #

import numpy as np

# 假设DOF=3，参数已定义：x10,x20,x30,v10,v20,v30,m1,m2,m3,c1,c2,k1,k2,k3,k4,f0,x_0,w
x0 = np.array([x10, x20, x30], dtype=float)
v0 = np.array([v10, v20, v30], dtype=float)
M = np.array([[m1, 0, 0], [0, m2, 0], [0, 0, m3]], dtype=float)
C = np.array([[c1+c2, -c1, -c2], [c1, -c2, 0], [c2, 0, -c2]], dtype=float)
K = np.array([[k1+k2, -k1, -k2], [k1, k3-k1, 0], [k2, 0, k4-k2]], dtype=float)
F0 = np.array([f0, -k3*x_0, -k4*x_0], dtype=float)

# 特征值求解：提取实部并过滤负值
D, V = np.linalg.eig(np.linalg.inv(M)@K)
D_real = np.real(D)
D_real[D_real < 0] = 0  # 避免平方根出现虚数
wn = np.sqrt(D_real)

# 归一化模态振型：用复数数组存储完整信息
A = np.zeros((3, 3), dtype=np.complex128)
for n in range(3):
    X = V[:, n]
    modal_mass = np.real(np.transpose(X)@M@X)  # 模态质量取实部确保正数
    b = 1/np.sqrt(modal_mass)
    A[:, n] = b*X

# 修正模态阻尼计算：使用矩阵乘法@
modal_damping_matrix = np.transpose(A)@C@A
zeta = np.diag(modal_damping_matrix)/(2*wn)
zeta = np.real(zeta)  # 提取实部作为阻尼系数
wd = wn*np.sqrt(1-zeta**2)

# 模态力与初始条件计算
u0 = np.transpose(A)@F0
qx0 = np.transpose(A)@M@x0
qv0 = np.transpose(A)@M@v0

# 受迫响应计算：确保arctan2输入为实数
Q0 = u0 / np.sqrt((wn**2 - w**2)**2 + (2*zeta*wn)**2)
phi = np.arctan2(np.real(-2*zeta*wn), np.real(wn**2 - w**2))

修改后的共用绘图代码 #

# 处理复数系数，提取实部用于响应计算
c1 = np.real(qx0 + Q0*np.sin(phi))
c2 = np.real(1/wd*(qv0 + zeta*wn*c1 - w*Q0*np.sin(phi)))

t = np.linspace(0, 0.1, 1000) 
qh = np.zeros([3, 1000], dtype=float)
qp = np.zeros([3, 1000], dtype=float)

for n in range(3):
    # 模态响应取实部，得到物理可解释的振动位移
    qh[n, :] = np.real(np.exp(-zeta[n]*wn[n]*t)*(c1[n]*np.cos(wd[n]*t)+c2[n]*np.sin(wd[n]*t)))
    qp[n, :] = np.real(Q0[n]*np.sin(w*t+phi[n]))

# 物理坐标响应取实部
xh = np.real(A@qh)
xp = np.real(A@qp)

# 绘图部分保持不变
import matplotlib.pyplot as plt

for n in range(3):
    plt.subplot(311)
    plt.plot(t, xh[n, :])
    plt.subplot(312)
    plt.plot(t, xp[n, :])
    plt.subplot(313)
    plt.plot(t, xh[n, :] + xp[n, :])

plt.subplot(311)
plt.ylabel('Free Vibrations')
plt.legend(['x1', 'x2', 'x3'], loc='upper right')
plt.title('Vibration Responses [m] of 3-DOF System')
plt.grid('on')
plt.xlim([0, 0.1])

plt.subplot(312)
plt.ylabel('Forced Vibrations')
plt.legend(['x1', 'x2', 'x3'], loc='upper right')
plt.grid('on')
plt.xlim([0, 0.1])

plt.subplot(313)
plt.ylabel('Total Vibrations')
plt.xlabel('Time [s]')
plt.legend(['x1', 'x2', 'x3'], loc='upper right')
plt.grid('on')
plt.xlim([0, 0.1])

plt.show()

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