<jupyter_start><jupyter_text>Define variables<jupyter_code>m, M, Ic, k, b, l, g = sm.symbols('m, M, I_c, k, b, l, g', real=True, positive=True)
t = sm.symbols('t', real=True)
x = sm.Function('x')(t)
v = sm.Function('v')(t)
theta = sm.Function('theta')(t)
omega = sm.Function('omega')(t)<jupyter_output><empty_output><jupyter_text>Kinetic energy<jupyter_code>vc_squared = (v - omega*l*sm.sin(theta))**2 + (omega*l*sm.cos(theta))**2
vc_squared
eq_repl = {theta.diff(t): omega, x.diff(t): v}
eq_repl
T = m*v**2/2 + M*vc_squared/2 + Ic*omega**2/2
T<jupyter_output><empty_output><jupyter_text>Potential energy<jupyter_code>U = k*x**2/2 -m*g*x - M*g*(x + l*sm.cos(theta))
U<jupyter_output><empty_output><jupyter_text>Damping<jupyter_code>R = b*v**2/2
R<jupyter_output><empty_output><jupyter_text>Lagrange's Equation<jupyter_code>L = T - U
L
leq_v = L.diff(v).diff(t).subs(eq_repl) - L.diff(x) -(-R.diff(v))
leq_v
leq_omega = L.diff(omega).diff(t).subs(eq_repl) - L.diff(theta) -(-R.diff(omega))
leq_omega
leq_omega = sm.simplify(leq_omega)
leq_omega<jupyter_output><empty_output><jupyter_text>Put EoM in explicit first order form<jupyter_code>f = sm.Matrix([leq_v, leq_omega])
f
g = f.subs({omega.diff(t): 0, v.diff(t): 0})
g
I = f.jacobian([omega.diff(t), v.diff(t)])
I
sdot = -I.inv()*g
sdot
sdot = -I.LUsolve(g)
sdot
sm.simplify(sdot)<jupyter_output><empty_output><jupyter_text>Simulate the non-linear system<jupyter_code>from resonance.nonlinear_systems import MultiDoFNonLinearSystem
sys = MultiDoFNonLinearSystem()
sys.constants['m'] = 1.0  # kg
sys.constants['k'] = 10.0  # N/m
sys.constants['b'] = 5.0  # Ns
sys.constants['l'] = 0.5  # m
sys.constants['M'] = 0.5  # kg
sys.constants['Ic'] = 0.5*0.5**2  # kg m**2
sys.constants['g'] = 9.81  # m/s**2
# order of entry matters!
sys.coordinates['theta'] = 1.0  # rad
sys.coordinates['x'] = 0.0  # m
sys.speeds['omega'] = 0.0 # rad/s
sys.speeds['v'] = 0.0 # m/s
sys.states
eval_mass = sm.lambdify((m, M), m+M)
eval_mass(1.0, 2.0)
g = sm.symbols('g', real=True, positive=True)
eval_sdot = sm.lambdify((theta, x, omega, v, m, k, b, l, M, Ic, g), [sdot[0], sdot[1]])
eval_sdot(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0)
def eval_derivatives(theta, x, omega, v, m, k, b, l, M, Ic, g):
    omegadot, vdot = eval_sdot(theta, x, omega, v, m, k, b, l, M, Ic, g)
    thetadot = omega
    xdot = v
    return thetadot, xdot, omegadot, vdot
eval_derivatives(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0)
sys.diff_eq_func = eval_derivatives
trajectories = sys.free_response(5.0)
%matplotlib widget
trajectories[['x', 'theta']].plot(subplots=True)<jupyter_output><empty_output><jupyter_text>Find the equilibrium<jupyter_code>sdot
sdot_no_motion = sdot.subs({omega: 0, v: 0})
sm.simplify(sdot_no_motion)
sm.solve(sdot_no_motion, theta, x)<jupyter_output><empty_output><jupyter_text>Linearize<jupyter_code>f
vbar = sm.Matrix([omega.diff(t), v.diff(t), omega, v, theta, x])
vbar
vbareq = sm.Matrix([0, 0, 0, 0, 0, g*(M + m)/k])
vbareq
vbar_repl = {k: v for k, v in zip(vbar, vbareq)}  # dictionary comphrension
vbar_repl
flin = f.subs(vbar_repl) + f.jacobian(vbar).subs(vbar_repl)*(vbar - vbareq)
flin
f.jacobian(vbar)<jupyter_output><empty_output><jupyter_text>Linear canonical form<jupyter_code>cbar = sm.Matrix([theta, x])
cbar
K = flin.jacobian(cbar)
K
sbar = sm.Matrix([omega, v])
sbar
C = flin.jacobian(sbar)
C
M_mat = flin.jacobian(sbar.diff(t))
M_mat
zero_repl = {k: v for k, v in zip(vbar, [0, 0, 0, 0, 0, 0])}
zero_repl

F = -flin.subs(zero_repl)
F<jupyter_output><empty_output><jupyter_text>Simulate the linear system<jupyter_code>from resonance.linear_systems import MultiDoFLinearSystem

sys = MultiDoFLinearSystem()

sys.constants['m'] = 1.0  # kg
sys.constants['k'] = 10.0  # N/m
sys.constants['b'] = 5.0  # Ns
sys.constants['l'] = 0.5  # m
sys.constants['M'] = 0.5  # kg
sys.constants['Ic'] = 0.5*0.5**2  # kg m**2
sys.constants['g'] = 9.81  # m/s**2
sys.coordinates['theta'] = 1.0  # rad
sys.coordinates['x'] = 0.0  # m
sys.speeds['omega'] = 0.0 # rad/s
sys.speeds['v'] = 0.0 # m/s
import numpy as np
np.array([[1, 2], [3, 4]])
def calc_coeff_mats(m, k, b, l, M, Ic, g):
    M_mat = np.array([[0, M+m],
                  [Ic+m*l**2, 0]])
    C = np.array([[0, b],
                  [0, 0]])
    K = np.array([[0, k],
                  [M*g*l, 0]])
    return M_mat, C, K
sys.canonical_coeffs_func = calc_coeff_mats
sys.canonical_coefficients()
sys.coordinates['x'] = 0.1
traj = sys.free_response(5.0)
traj[['x', 'theta']].plot(subplots=True)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19872359154929578, 'fraction_numerical': 0.03433098591549296, 'mean_word_length': 4.812020460358056, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 45, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '83175', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 2242, 'n_tokens_neox': 2128, 'n_words': 398, 'path': 'content/materials/notebooks/2020/l11_linearizing_a_bouncy_pendulum.ipynb', 'repo_name': 'moorepants/eng122', 'size': '4544'}