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

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

First, collect the data:

In [12]:

deffind_bifurcation(k):eigs=[]ds=linspace(0,100,200)fordinds:M=system_matrix(d,k)eigs.append(eigvals(M))# make an array:aeigs=array(eigs)# let's plot the imaginary parts:plot(ds,aeigs.imag,'.')# compute the part where the imaginary part is close to zero:mask=aeigs.imag<1e-7# we need all eigenvalues to be real:mask_=all(mask,axis=1)# at which coefficient to we switch from False to True?i=argmax(mask_)# which d was that:returnds[i]find_bifurcation(1000)

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.