Frequency response.ipynb

{
 "worksheets": [
  {
   "cells": [
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Exercise for the course [Python for MATLAB users](http://sese.nu/python-for-matlab-users/)",
     "level": 2
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Original exercise by Claus Führer, modified by Olivier Verdier"
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "%pylab\n%matplotlib inline\nfrom __future__ import division",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "-----"
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "System matrix",
     "level": 2
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Consider the matrix\n\\\\[\n A=\\begin{bmatrix}\n  0 & I \\\\\nK & D\n \\end{bmatrix}\n\\\\]\nwhere $0$ and $I$ are the $2 \\times 2$ zero and identity matrices and $K$ and $D$ are $2 \\times 2$\nmatrices of the following form:\n\\\\[\n K=\\begin{bmatrix}\n    -k & 0.5 \\\\ 0.5 & -k\n   \\end{bmatrix}\n\\qquad\nD=\\begin{bmatrix}\n    -d & 1.0 \\\\ 1.0 & -d\n   \\end{bmatrix}\n\\\\]\nwith $k$ and $d$ being real parameters.\n\n"
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Write a function `stiffness` which constructs the matrix $K$ above.",
     "level": 3
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "def stiffness(k):\n    return zeros([2,2]) # implement this!",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "assert(allclose(stiffness(1.), array([[-1.,.5],[.5,-1.]])))",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Write a function `damping` which constructs the matrix $D$ above.",
     "level": 3
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "def damping(d):\n    return zeros([2,2]) # implement this!",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "assert(allclose(damping(1.), array([[-1.,1.],[1.,-1.]])))",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Write a function `system_matrix` which takes $k$ and $d$ as input and which generates the matrix $A$.",
     "level": 3
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "\nHint: use the function `concatenate`. Check its documentation by running:"
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "concatenate?",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Use also `identity` (or `eye`), `zeros` (or `zeros_like`)."
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "def system_matrix(d, k):\n    return zeros([4,4]) # implement this!",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "A = system_matrix(10.,20.)\nassert(allclose(A[:2,:2], zeros([2,2])))\nassert(allclose(A[:2,2:4], identity(2)))\nassert(allclose(A[2:4,:2], stiffness(20.)))\nassert(allclose(A[2:4,2:4], damping(10.)))",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "For samples of the values $d \\in [0,100]$ and the fixed value $k=1000$, plot the four eigenvalues on the complex plane.",
     "level": 3
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "for d in linspace(0,100,200):\n    pass",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Bonus question: there is a bifurcation in the diagram above. Can you find the bifurcation point programmatically?",
     "level": 3
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "heading",
     "source": "Frequency Response Plot",
     "level": 2
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "In technical applications there occurs often linear systems of the form\n\\\\[\n\\dot x(t) = A x(t) + B u(t)\n\\\\]\nwhere $u$ is an given input signal. $x$ is called the state. From the state some quantities $y(t)$ can be\nmeasured, this is decribed by the equation\n\\\\[\n y(t)=C x(t).\n\\\\]\nWe 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. \n\nThe relationship between the in- and out-amplitude is given by the formula\n\\\\[\n \\mathrm{amplitude}(\\omega)=\\\\|(G(i\\omega))\\\\|\\quad\\text{where}\n\\quad G(i\\omega)=C \\cdot (i\\omega I -A)^{-1} \\cdot B\n\\\\]\nand $i$ is the imaginary unit.\n\n"
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "inv?",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "norm?",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "def amplitude(A, B, C, omega):\n    pass",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Plot the amplitude versus omega, for $\\omega \\in [0, 160]$, with $A$ being the system matrix above with $d=20$ and $k=500$, and\n\\\\[\n C=\\begin{bmatrix} 1 & 0 & 0 & 0  \\end{bmatrix} \\qquad B=\\begin{bmatrix}0 \\\\ 0 \\\\ 0\\\\ 1 \\end{bmatrix} .\n\\\\]\n"
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    },
    {
     "metadata": {},
     "cell_type": "markdown",
     "source": "Find out the relationship between $A$'s eigenvalues and the peak(s) in the figure.\n"
    },
    {
     "metadata": {},
     "cell_type": "code",
     "input": "",
     "outputs": [],
     "language": "python",
     "trusted": true,
     "collapsed": false
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "gist_id": "708b799a2d96af9bc96f",
  "name": "",
  "signature": "sha256:0eb21b1ee169b5ea6eefb141d2dc2c79b6c7abefaad88179ca4758e9fb6560ec"
 },
 "nbformat": 3
}