【问题标题】:Running into error while trying to solve an matrix eigen value problem using numpy that I previously did not run into尝试使用我以前没有遇到过的 numpy 解决矩阵特征值问题时遇到错误
【发布时间】:2023-02-23 05:24:38
【问题描述】:

这是我得到的错误: " phi = arctan2(-2泽塔万,万2-w2)

TypeError:输入类型不支持 ufunc 'arctan2',并且根据转换规则 ''safe'' 无法将输入安全地强制转换为任何受支持的类型” 另外我收到这条消息: " ComplexWarning: 将复数值转换为实数会丢弃虚数部分 A[:, n] = b*X 回溯(最后一次通话):”

我正在尝试使用 numpy eig、inv、transpose、arctan2 等解决三自由度弹簧阻尼器问题。我做了一个以前的问题,我能够输出一个显示强制、自由和总振动的图形模型。我最初没有遇到任何问题,现在我试图使用代码来绘制对不同问题的响应,我在 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()

【问题讨论】:

  • 确定问题数组及其数据类型

标签: python numpy matrix transpose numpy-ufunc


【解决方案1】:

复杂的参数会引发此错误

In [92]: np.arctan2(1,2j)
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''

【讨论】:

    猜你喜欢
    • 2017-06-20
    • 2020-08-01
    • 2014-10-28
    • 2021-09-27
    • 2014-10-09
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多