这是我得到的错误：“φ =反正切2（-2* ζ wn，wn2-w2）
TypeError：输入类型不支持ufunc 'arctan 2'，并且根据强制转换规则“”safe "“，无法将输入安全地强制为任何支持的类型”“此外，我还收到以下消息：“复杂警告：将复数值转换为真实的会丢弃虚部A[：，n] = BX Traceback（最近的最后调用）：“
我正在尝试使用numpy eig、inv、transpose、arctan 2等解决一个三自由度Spring阻尼器问题。我之前做过一个问题，可以输出一个图形模型，显示受迫振动、自由振动和总振动。最初我没有得到任何一个问题，现在我尝试使用代码来绘制不同问题的响应。我在Spyder上收到了这两条消息。我将发布相关代码来展示我的过程。我所做的修改只是初始值、初始边界条件和输入函数，以反映在执行FBD并将EOM转换为状态空间矩阵形式后的问题。
------以前的代码配置----------

x0 = array([x10, x20, x30], dtype=float)
v0 = array([v10, v20, v30], dtype=float)
M = array([[m1, 0, 0], [0, m2, 0], [0, 0, m3]], dtype=float)
C = array([[c1, -c1, 0], [-c1, c1+c2, -c2], [0, -c2, c2]], dtype=float)
K = array([[k1+k2, -k2, 0], [-k2, k2+k3, -k3], [0, -k3, k3]], dtype=float)
F0 = array([0, 0, f0], dtype=float)
# Eigenvalue problem
D, V = eig(inv(M)@K)
wn = sqrt(D)
# Normalization of mode shapes w.r.t. the mass matrix
A = zeros((DOF, DOF), dtype=float)
for n in range(DOF):
    X = V[:, n]
    b = 1/sqrt(transpose(X)@M@X)
    A[:, n] = b*X
# Modal damping factors and damped natural angular frequenices
zeta = diag(transpose(A)@C*A)/(2*wn)
wd = wn*sqrt(1-zeta**2)
# Modal force vector
u0 = transpose(A)@F0
# Initial conditions in the modal coordinates
qx0 = transpose(A)@M@x0
qv0 = transpose(A)@M@v0
# Forced response amplitudes and phase angles
Q0 = u0/sqrt((wn**2-w**2)**2 + (2*zeta*wn)**2)
phi = arctan2(-2*zeta*wn, wn**2-w**2)

• -------------新代码配置------------------

x0 = array([x10, x20, x30], dtype=float)
v0 = array([v10, v20, v30], dtype=float)
M = array([[m1, 0, 0], [0, m2, 0], [0, 0, m3]], dtype=float)
C = array([[c1+c2, -c1, -c2], [c1, -c2, 0], [c2, 0, -c2]], dtype=float)
K = array([[k1+k2, -k1, -k2], [k1, k3-k1, 0], [k2, 0, k4-k2]], dtype=float)
F0 = array([f0, -k3*x_0, -k4*x_0], dtype=float)
# Eigenvalue problem
D, V = eig(inv(M)@K)
wn = sqrt(D)
# Normalization of mode shapes w.r.t. the mass matrix
A = zeros((DOF, DOF), dtype=float)
for n in range(DOF):
    X = V[:, n]
    b = 1/sqrt(transpose(X)@M@X)
    A[:, n] = b*X
# Modal damping factors and damped natural angular frequenices
zeta = diag(transpose(A)@C*A)/(2*wn)
wd = wn*sqrt(1-zeta**2)
# Modal force vector
u0 = transpose(A)@F0
# Initial conditions in the modal coordinates
qx0 = transpose(A)@M@x0
qv0 = transpose(A)@M@v0
# Forced response amplitudes and phase angles
Q0 = u0/sqrt((wn**2-w**2)**2 + (2*zeta*wn)**2)
phi = arctan2(-2*zeta*wn, wn**2-w**2)

我只是替换了值，使矩阵反映了我的新问题，现在我遇到了我不知道如何解决的问题。
------------------对两者都相同的最后一位代码------------------------

# Unknown coefficients in the free vibration responses
c1 = qx0 + Q0*sin(phi)
c2 = 1/wd*(qv0+zeta*wn*c1-w*Q0*sin(phi))
# Modal responses
t = linspace(0, 0.1, 1000) 
qh = zeros([DOF, 1000], dtype=float)
qp = zeros([DOF, 1000], dtype=float)
for n in range(DOF):
    qh[n, :] = exp(-zeta[n]*wn[n]*t)*(c1[n]*cos(wd[n]*t)+c2[n]*sin(wd[n]*t))
    qp[n, :] = Q0[n]*sin(w*t+phi[n])
# Responses in the physical coordinates
xh = A@qh
xp = A@qp
# Plots

for n in range(DOF):
    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()