Exercise for the course Python for MATLAB users

Original exercise by Claus F├╝hrer, modified by Olivier Verdier

In [ ]:
%pylab
%matplotlib inline

System matrix

Consider the matrix $$ A=\begin{bmatrix} 0 & I \\ K & D \end{bmatrix} $$ where $0$ and $I$ are the $2 \times 2$ zero and identity matrices and $K$ and $D$ are $2 \times 2$ matrices of the following form: $$ K=\begin{bmatrix} -k & 0.5 \\ 0.5 & -k \end{bmatrix} \qquad D=\begin{bmatrix} -d & 1.0 \\ 1.0 & -d \end{bmatrix} $$ with $k$ and $d$ being real parameters.

Write a function stiffness which constructs the matrix $K$ above.

In [ ]:
def stiffness(k):
    return zeros([2,2]) # implement this!
In [ ]:
assert(allclose(stiffness(1.), array([[-1.,.5],[.5,-1.]])))

Write a function damping which constructs the matrix $D$ above.

In [ ]:
def damping(d):
    return zeros([2,2]) # implement this!
In [ ]:
assert(allclose(damping(1.), array([[-1.,1.],[1.,-1.]])))

Write a function system_matrix which takes $k$ and $d$ as input and which generates the matrix $A$.

Hint: use the function concatenate. Check its documentation by running:

In [ ]:
concatenate?

Use also identity (or eye), zeros (or zeros_like).

In [ ]:
def system_matrix(d, k):
    return zeros([4,4]) # implement this!
In [ ]:
A = system_matrix(10.,20.)
assert(allclose(A[:2,:2], zeros([2,2])))
assert(allclose(A[:2,2:4], identity(2)))
assert(allclose(A[2:4,:2], stiffness(20.)))
assert(allclose(A[2:4,2:4], damping(10.)))

For samples of the values $d \in [0,100]$ and the fixed value $k=1000$, plot the four eigenvalues on the complex plane.

In [ ]:
for d in linspace(0,100,200):
    pass

Bonus question: there is a bifurcation in the diagram above. Can you find the bifurcation point programmatically?

Frequency Response Plot

In technical applications there occurs often linear systems of the form $$ \dot x(t) = A x(t) + B u(t) $$ where $u$ is an given input signal. $x$ is called the state. From the state some quantities $y(t)$ can be measured, this is decribed by the equation $$ y(t)=C x(t). $$ We assume here that the input signal is an harmonic oscillation $u(t)=\sin(\omega t)$ with a given frequency $\omega$ and amplitude one. Then, $y(t)$ is again a harmonic oscillation with the same frequency, but another amplitude. The amplitude depends on the frequency.

The relationship between the in- and out-amplitude is given by the formula $$ \mathrm{amplitude}(\omega)=\\|(G(i\omega))\\|\quad\text{where} \quad G(i\omega)=C \cdot (i\omega I -A)^{-1} \cdot B $$ and $i$ is the imaginary unit.

In [ ]:
inv?
In [ ]:
norm?
In [ ]:
def amplitude(A, B, C, omega):
    pass

Plot the amplitude versus omega, for $\omega \in [0, 160]$, with $A$ being the system matrix above with $d=20$ and $k=500$, and $$ C=\begin{bmatrix} 1 & 0 & 0 & 0 \end{bmatrix} \qquad B=\begin{bmatrix}0 \\ 0 \\ 0\\ 1 \end{bmatrix} . $$

Find out the relationship between $A$'s eigenvalues and the peak(s) in the figure.