(Q3)Use sympy to determine closed-form (symbolic) expressions for the impulse, step and frequency response of Gθ and Gx

PendulumQuestion 3.py

import sympy as sym

m, ell, x1, x2, x3, x4, M, g, F = sym.symbols('m, ell, x1, x2, x3, x4, M, g, F')
# φ(F, x3, x4)
phi = 4*m*ell*x4**2*sym.sin(x3) + 4*F - 3*m*g*sym.sin(x3)*sym.cos(x3)
phi /= 4*(M+m) - 3*m*sym.cos(x3)**2
dphi_x3 = phi.diff(x3)
dphi_x4 = phi.diff(x4)
dphi_F = phi.diff(F)

psi = -3*(m*ell*x4**2*sym.sin(x3)*sym.cos(x3) + F*sym.cos(x3) - (M+m)*g*sym.sin(x3))
psi /= (4*(M+m) - 3*m*sym.cos(x3)**2)*ell
dpsi_x3 = psi.diff(x3)
dpsi_x4 = psi.diff(x4)
dpsi_F = psi.diff(F)
# Equilibrium point
Feq = 0
x3eq = 0
x4eq = 0
#------------------------------
#  adding equilibrium to x1 and x2
#------------------------------
x1eq = 0
x2eq = 0

dphi_F_eq = dphi_F.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dphi_x3_eq = dphi_x3.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dphi_x4_eq = dphi_x4.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])

dpsi_F_eq = dpsi_F.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dpsi_x3_eq = dpsi_x3.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])
dpsi_x4_eq = dpsi_x4.subs([(F, Feq), (x3, x3eq), (x4, x4eq)])

a = dphi_F_eq
b = -dphi_x3_eq
c = -dpsi_F_eq
d = dpsi_x3_eq

# -----------------------
a, b, c, d = sym.symbols('a:d', real=True, positive=True)
s, t = sym.symbols('s, t')
transfer_function_F_to_x3 = -c/(s**2 - d)
transfer_function_F_to_x1 = ((a*s**2) - (a*d) + (b*c))/(s**4 - (d * s**2))
#impulse Response(kick)
Xs_impulse = 1

#G_theta response
x3t_impulse = sym.inverse_laplace_transform(transfer_function_F_to_x3 * Xs_impulse, s, t)
sym.pprint(x3t_impulse.simplify())
#G_x response
x1t_impulse = sym.inverse_laplace_transform(transfer_function_F_to_x1 * Xs_impulse, s, t)
sym.pprint(x1t_impulse.simplify())

#Step Response(Push)
Xs_step = 1/s
#G_theta response
x3t_step = sym.inverse_laplace_transform(transfer_function_F_to_x3 * Xs_step, s, t)
sym.pprint(x3t_step.simplify())
#G_x response
x1t_step = sym.inverse_laplace_transform(transfer_function_F_to_x1 * Xs_step, s, t)
sym.pprint(x1t_step.simplify())

#Frequency Response(shake)
w = sym.symbols('w', real=True)
Xs_shake = w/(s**2 + w**2)

#G_theta response
x3t_shake = sym.inverse_laplace_transform(transfer_function_F_to_x3 * Xs_shake, s, t, w)
sym.pprint(x3t_shake.simplify())
#G_x response
x1t_shake = sym.inverse_laplace_transform(transfer_function_F_to_x1 * Xs_shake, s, t, w)
sym.pprint(x1t_shake.simplify())

print(sym.latex(x3t_impulse.simplify()))
print(sym.latex(x1t_impulse.simplify()))
print(sym.latex(x3t_step.simplify()))
print(sym.latex(x1t_step.simplify()))
print(sym.latex(x3t_shake.simplify()))
print(sym.latex(x1t_shake.simplify()))