```python
import sympy as sm
```


```python
sm.init_printing()
```

# Define variables


```python
m, M, Ic, k, b, l, g = sm.symbols('m, M, I_c, k, b, l, g', real=True, positive=True)
```


```python
t = sm.symbols('t', real=True)
```


```python
x = sm.Function('x')(t)
v = sm.Function('v')(t)
theta = sm.Function('theta')(t)
omega = sm.Function('omega')(t)
```

# Kinetic energy


```python
vc_squared = (v - omega*l*sm.sin(theta))**2 + (omega*l*sm.cos(theta))**2
vc_squared
```


```python
eq_repl = {theta.diff(t): omega, x.diff(t): v}
eq_repl
```


```python
T = m*v**2/2 + M*vc_squared/2 + Ic*omega**2/2
T
```

# Potential energy


```python
U = k*x**2/2 -m*g*x - M*g*(x + l*sm.cos(theta))
U
```

# Damping


```python
R = b*v**2/2
R
```

# Lagrange's Equation


```python
L = T - U
L
```


```python
leq_v = L.diff(v).diff(t).subs(eq_repl) - L.diff(x) -(-R.diff(v))
leq_v
```


```python
leq_omega = L.diff(omega).diff(t).subs(eq_repl) - L.diff(theta) -(-R.diff(omega))
leq_omega
```


```python
leq_omega = sm.simplify(leq_omega)
leq_omega
```

# Put EoM in explicit first order form


```python
f = sm.Matrix([leq_v, leq_omega])
f
```


```python
g = f.subs({omega.diff(t): 0, v.diff(t): 0})
g
```


```python
I = f.jacobian([omega.diff(t), v.diff(t)])
I
```


```python
sdot = -I.inv()*g
sdot
```


```python
sdot = -I.LUsolve(g)
sdot
```


```python
sm.simplify(sdot)
```

# Simulate the non-linear system


```python
from resonance.nonlinear_systems import MultiDoFNonLinearSystem
```


```python
sys = MultiDoFNonLinearSystem()
```


```python
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
```


```python
# 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
```


```python
sys.states
```




    _StatesDict([('theta', 1.0), ('x', 0.0), ('omega', 0.0), ('v', 0.0)])




```python
eval_mass = sm.lambdify((m, M), m+M)
```


```python
eval_mass(1.0, 2.0)
```


```python
g = sm.symbols('g', real=True, positive=True)
```


```python
eval_sdot = sm.lambdify((theta, x, omega, v, m, k, b, l, M, Ic, g), [sdot[0], sdot[1]])
```


```python
eval_sdot(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0)
```


```python
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
```


```python
eval_derivatives(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0)
```


```python
sys.diff_eq_func = eval_derivatives
```


```python
trajectories = sys.free_response(5.0)
```


```python
%matplotlib widget
```


```python
trajectories[['x', 'theta']].plot(subplots=True)
```


    Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …





    array([<matplotlib.axes._subplots.AxesSubplot object at 0x7f8dab65f4a8>,
           <matplotlib.axes._subplots.AxesSubplot object at 0x7f8dab5d3be0>],
          dtype=object)



# Find the equilibrium


```python
sdot
```


```python
sdot_no_motion = sdot.subs({omega: 0, v: 0})
```


```python
sm.simplify(sdot_no_motion)
```


```python
sm.solve(sdot_no_motion, theta, x)
```

# Linearize


```python
f
```


```python
vbar = sm.Matrix([omega.diff(t), v.diff(t), omega, v, theta, x])
vbar
```


```python
vbareq = sm.Matrix([0, 0, 0, 0, 0, g*(M + m)/k])
vbareq
```


```python
vbar_repl = {k: v for k, v in zip(vbar, vbareq)}  # dictionary comphrension
vbar_repl
```


```python
flin = f.subs(vbar_repl) + f.jacobian(vbar).subs(vbar_repl)*(vbar - vbareq)
flin
```


```python
f.jacobian(vbar)
```

# Linear canonical form


```python
cbar = sm.Matrix([theta, x])
cbar
```


```python
K = flin.jacobian(cbar)
K
```


```python
sbar = sm.Matrix([omega, v])
sbar
```


```python
C = flin.jacobian(sbar)
C
```


```python
M_mat = flin.jacobian(sbar.diff(t))
M_mat
```


```python
zero_repl = {k: v for k, v in zip(vbar, [0, 0, 0, 0, 0, 0])}
zero_repl

```


```python
F = -flin.subs(zero_repl)
F
```

# Simulate the linear system


```python
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
```


```python
import numpy as np
```


```python
np.array([[1, 2], [3, 4]])
```




    array([[1, 2],
           [3, 4]])




```python
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
```


```python
sys.canonical_coeffs_func = calc_coeff_mats
```


```python
sys.canonical_coefficients()
```




    (array([[0.   , 1.5  ],
            [0.375, 0.   ]]), array([[0., 5.],
            [0., 0.]]), array([[ 0.    , 10.    ],
            [ 2.4525,  0.    ]]))




```python
sys.coordinates['x'] = 0.1
```


```python
traj = sys.free_response(5.0)
```


```python
traj[['x', 'theta']].plot(subplots=True)
```


    Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …





    array([<matplotlib.axes._subplots.AxesSubplot object at 0x7f8dab554f60>,
           <matplotlib.axes._subplots.AxesSubplot object at 0x7f8dab444ba8>],
          dtype=object)



59f4f7902ac43582971ca74b9a85821585a27543
ipynb
Jupyter Notebook
content/materials/notebooks/2020/l11_linearizing_a_bouncy_pendulum.ipynb
moorepants/eng122
7bcd502eae4ab0ec9d463389acca2290b263ba64
CC-BY-4.0
2020-05-11T05:56:54.000Z
2020-12-20T13:44:12.000Z
content/materials/notebooks/2020/l11_linearizing_a_bouncy_pendulum.ipynb
moorepants/eng122
7bcd502eae4ab0ec9d463389acca2290b263ba64
CC-BY-4.0
2016-09-20T01:26:33.000Z
2020-01-22T02:10:33.000Z
content/materials/notebooks/2020/l11_linearizing_a_bouncy_pendulum.ipynb
moorepants/eng122
7bcd502eae4ab0ec9d463389acca2290b263ba64
CC-BY-4.0
2017-06-09T01:28:45.000Z
2020-02-29T02:36:18.000Z
Qwen/Qwen-72B
1. YES
2. YES