This is the error I am getting: " phi = arctan2(-2zetawn, wn2-w2)
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'' " In addition I am getting this message: " ComplexWarning: Casting complex values to real discards the imaginary part A[:, n] = b*X Traceback (most recent call last):"
I am trying to solve a three degree of freedom spring damper problem using numpy eig, inv, transpose, arctan2, etc. I did a previous problem where I was able to output a graphical model showing the forced, free, and total vibrations. I was not getting either issue originally and now that I am attempting to use the code to graph the response on a different problem, I am getting both messages on Spyder. I will post relevant code to show my process. All I changed were the initial values, initial boundary conditions, and the input function to reflect the problem after doing a FBD and putting the EOM's into state space matrix form.
-------previous code config-----------
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)
---------------New code config--------------------------------
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)
I just replaced the values and made the matrix reflect my new problem and now I am running into issues I don't know how to fix.
-------------------Last bit of code that is the same for both---------------------------
# 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()
