Mass Spring Damper Example¶

The equation of motion is given as:

$m \ddot{x}(t) + c \dot{x}(t) + k x(t) = u(t)$

You are given the following physical parameters.

• c=2 N/(m/s)
• k=2 N/m
• m=1 kg

Laplace Method for Matrix Exponential¶

In [12]:
import matplotlib.pyplot as plt
%matplotlib inline
import sympy

sympy.init_printing()
t = sympy.symbols('t', real=True)
eps = 1e-10  # where to start plotting, since we assume positive t
s = sympy.symbols('s')
c = 1
k = 1

A = sympy.Matrix([[0, 1], [-k, -c]])
B = sympy.Matrix([[0], [1]])
C = sympy.Matrix([[1, 0]])
D = sympy.Matrix([[0]])
(A, B, C, D)

Out[12]:
$\displaystyle \left( \left[\begin{matrix}0 & 1\\-1 & -1\end{matrix}\right], \ \left[\begin{matrix}0\\1\end{matrix}\right], \ \left[\begin{matrix}1 & 0\end{matrix}\right], \ \left[\begin{matrix}0\end{matrix}\right]\right)$

$G(s) = C(sI - A)^{-1}B + D$

In [13]:
N = (s*sympy.eye(2) - A).inv()
G = (C*N*B + D).expand()
G.simplify()
G = G[0]
G

Out[13]:
$\displaystyle \frac{1}{s^{2} + s + 1}$

Remember from the Laplace table:

$\mathcal{L}^{-1}\left[ \dfrac{1}{s - a} \right] = e^{at}$

The same is true for the matrix version of this Laplace transform:

$\mathcal{L}^{-1}\left[ (sI- A)^{-1} \right] = e^{At}$

In [14]:
eAt = sympy.inverse_laplace_transform((s*sympy.eye(2) - A).inv(), s, t)
eAt

Out[14]:
$\displaystyle \left[\begin{matrix}\frac{\sqrt{3} \left(\sin{\left(\frac{\sqrt{3} t}{2} \right)} + \sqrt{3} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) e^{- \frac{t}{2}} \theta\left(t\right)}{3} & \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3}\\- \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3} & \frac{\sqrt{3} \left(- \sin{\left(\frac{\sqrt{3} t}{2} \right)} + \sqrt{3} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) e^{- \frac{t}{2}} \theta\left(t\right)}{3}\end{matrix}\right]$
In [15]:
sympy.plot(sympy.Heaviside(t), (t, -1, 1))

Out[15]:
<sympy.plotting.plot.Plot at 0x7f74c3dd4d00>

Spectral Method for Matrix Exponential¶

The approach here will be to use the diagonalization of A.

In [16]:
A

Out[16]:
$\displaystyle \left[\begin{matrix}0 & 1\\-1 & -1\end{matrix}\right]$
In [17]:
evects = sympy.Matrix.eigenvects(A)
evects

Out[17]:
$\displaystyle \left[ \left( - \frac{1}{2} - \frac{\sqrt{3} i}{2}, \ 1, \ \left[ \left[\begin{matrix}- \frac{2}{1 + \sqrt{3} i}\\1\end{matrix}\right]\right]\right), \ \left( - \frac{1}{2} + \frac{\sqrt{3} i}{2}, \ 1, \ \left[ \left[\begin{matrix}- \frac{2}{1 - \sqrt{3} i}\\1\end{matrix}\right]\right]\right)\right]$

$T = \left[ \xi_1 \xi_2 \right]$

Diagonolization:

$A = T diag(\lambda_1, \lambda_2) T^{-1}$

In [18]:
T = sympy.Matrix.hstack(evects[0][2][0], evects[1][2][0])
T

Out[18]:
$\displaystyle \left[\begin{matrix}- \frac{2}{1 + \sqrt{3} i} & - \frac{2}{1 - \sqrt{3} i}\\1 & 1\end{matrix}\right]$
In [19]:
T_inv = T.inv()
T_inv.simplify()
T_inv

Out[19]:
$\displaystyle \left[\begin{matrix}\frac{\frac{\sqrt{3}}{3} - i}{\sqrt{3} + i} & \frac{1}{2} - \frac{\sqrt{3} i}{6}\\\frac{\sqrt{3} i}{3} & \frac{1}{2} + \frac{\sqrt{3} i}{6}\end{matrix}\right]$
In [20]:
Lambda = sympy.diag(evects[0][0], evects[1][0])
Lambda

Out[20]:
$\displaystyle \left[\begin{matrix}- \frac{1}{2} - \frac{\sqrt{3} i}{2} & 0\\0 & - \frac{1}{2} + \frac{\sqrt{3} i}{2}\end{matrix}\right]$

Let's check that our diagonalization of A worked as expected.

In [21]:
A_check = (T*Lambda*T.inv())
sympy.expand(A_check, complex=True), A

Out[21]:
$\displaystyle \left( \left[\begin{matrix}0 & 1\\-1 & -1\end{matrix}\right], \ \left[\begin{matrix}0 & 1\\-1 & -1\end{matrix}\right]\right)$

We can calculate the exponential of the diagonal matrix containing the eigen values, it is just the exponential of each element of the diagonal.

In [22]:
LambdaExp = sympy.diag(sympy.exp(evects[0][0]*t), sympy.exp(evects[1][0]*t))
LambdaExp

Out[22]:
$\displaystyle \left[\begin{matrix}e^{t \left(- \frac{1}{2} - \frac{\sqrt{3} i}{2}\right)} & 0\\0 & e^{t \left(- \frac{1}{2} + \frac{\sqrt{3} i}{2}\right)}\end{matrix}\right]$
In [23]:
eAt_spec = (T*LambdaExp*T_inv).expand(complex=True)*sympy.Heaviside(t)
eAt_spec

Out[23]:
$\displaystyle \left[\begin{matrix}\left(\frac{\sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)}}{3} + e^{- \frac{t}{2}} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) \theta\left(t\right) & \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3}\\- \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3} & \left(- \frac{\sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)}}{3} + e^{- \frac{t}{2}} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) \theta\left(t\right)\end{matrix}\right]$
In [24]:
eAt

Out[24]:
$\displaystyle \left[\begin{matrix}\frac{\sqrt{3} \left(\sin{\left(\frac{\sqrt{3} t}{2} \right)} + \sqrt{3} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) e^{- \frac{t}{2}} \theta\left(t\right)}{3} & \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3}\\- \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3} & \frac{\sqrt{3} \left(- \sin{\left(\frac{\sqrt{3} t}{2} \right)} + \sqrt{3} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) e^{- \frac{t}{2}} \theta\left(t\right)}{3}\end{matrix}\right]$

Response given initial conditions¶

$\vec{x} = \begin{bmatrix}x \\ \dot{x} \end{bmatrix}$

$\vec{x}(t) = e^{At}\vec{x}(0) + \int_0^t e^{A(t - \tau)} B u(\tau) d\tau$

$\vec{y}(t) = C \vec{x}(t) + D\vec{u}(t)$

In [25]:
x0 = sympy.Matrix([1, 0])
x0

Out[25]:
$\displaystyle \left[\begin{matrix}1\\0\end{matrix}\right]$
In [26]:
x = eAt*x0
x

Out[26]:
$\displaystyle \left[\begin{matrix}\frac{\sqrt{3} \left(\sin{\left(\frac{\sqrt{3} t}{2} \right)} + \sqrt{3} \cos{\left(\frac{\sqrt{3} t}{2} \right)}\right) e^{- \frac{t}{2}} \theta\left(t\right)}{3}\\- \frac{2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3}\end{matrix}\right]$
In [27]:
sympy.plot(x[0], (t, eps, 10))

Out[27]:
<sympy.plotting.plot.Plot at 0x7f753055aac0>
In [28]:
sympy.plot(x[1])

Out[28]:
<sympy.plotting.plot.Plot at 0x7f7530cd0550>

Final Value Theorem¶

$G(s) = Y(s)/U(s)$

$\lim\limits_{t \rightarrow \infty} x(t) = \lim\limits_{s \rightarrow 0} s X(s)$

In [29]:
G

Out[29]:
$\displaystyle \frac{1}{s^{2} + s + 1}$

Unit step of magnitude 3¶

$X = G\cdot U$

$X = (X/U)\cdot U$

In [30]:
U1 = 3/s
X1 = G*U1
X1

Out[30]:
$\displaystyle \frac{3}{s \left(s^{2} + s + 1\right)}$
In [31]:
u1 = sympy.inverse_laplace_transform(U1, s, t)
u1

Out[31]:
$\displaystyle 3 \theta\left(t\right)$
In [32]:
x1_final = sympy.simplify(s*X1).subs(s, 0)
x1_final

Out[32]:
$\displaystyle 3$
In [33]:
x1 = sympy.inverse_laplace_transform(sympy.apart(X1, full=True), s, t)
x1

Out[33]:
$\displaystyle 3 \theta\left(t\right) - 2 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} + \frac{\pi}{3} \right)} \theta\left(t\right)$
In [34]:
sympy.plot(x1, x1_final, (t, eps, 10))

Out[34]:
<sympy.plotting.plot.Plot at 0x7f7530ff7610>

Impulse of magnitude 2¶

In [35]:
U2 = 2
X2 = G*U2
X2

Out[35]:
$\displaystyle \frac{2}{s^{2} + s + 1}$
In [36]:
u2 = sympy.inverse_laplace_transform(U2, s, t)
u2

Out[36]:
$\displaystyle 2 \delta\left(t\right)$
In [37]:
x2_final = (s*X2).subs(s, 0)
x2_final

Out[37]:
$\displaystyle 0$
In [38]:
s*X2

Out[38]:
$\displaystyle \frac{2 s}{s^{2} + s + 1}$
In [39]:
x2 = sympy.inverse_laplace_transform(X2, s, t).expand(complex=True)
x2

Out[39]:
$\displaystyle \frac{4 \sqrt{3} e^{- \frac{t}{2}} \sin{\left(\frac{\sqrt{3} t}{2} \right)} \theta\left(t\right)}{3}$
In [40]:
sympy.plot(x2, x2_final, (t, eps, 10))

Out[40]:
<sympy.plotting.plot.Plot at 0x7f75310aafa0>