Frequency response.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise for the course [Python for MATLAB users](http://sese.nu/python-for-matlab-users-ht15/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Original exercise by Claus Führer, modified by Olivier Verdier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%pylab\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## System matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the matrix\n",
    "\\\\[\n",
    " A=\\begin{bmatrix}\n",
    "  0 & I \\\\\n",
    "K & D\n",
    " \\end{bmatrix}\n",
    "\\\\]\n",
    "where $0$ and $I$ are the $2 \\times 2$ zero and identity matrices and $K$ and $D$ are $2 \\times 2$\n",
    "matrices of the following form:\n",
    "\\\\[\n",
    " K=\\begin{bmatrix}\n",
    "    -k & 0.5 \\\\ 0.5 & -k\n",
    "   \\end{bmatrix}\n",
    "\\qquad\n",
    "D=\\begin{bmatrix}\n",
    "    -d & 1.0 \\\\ 1.0 & -d\n",
    "   \\end{bmatrix}\n",
    "\\\\]\n",
    "with $k$ and $d$ being real parameters.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write a function `stiffness` which constructs the matrix $K$ above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def stiffness(k):\n",
    "    return zeros([2,2]) # implement this!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert(allclose(stiffness(1.), array([[-1.,.5],[.5,-1.]])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write a function `damping` which constructs the matrix $D$ above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def damping(d):\n",
    "    return zeros([2,2]) # implement this!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "assert(allclose(damping(1.), array([[-1.,1.],[1.,-1.]])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Write a function `system_matrix` which takes $k$ and $d$ as input and which generates the matrix $A$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Hint: use the function `concatenate`. Check its documentation by running:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "concatenate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use also `identity` (or `eye`), `zeros` (or `zeros_like`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def system_matrix(d, k):\n",
    "    return zeros([4,4]) # implement this!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "A = system_matrix(10.,20.)\n",
    "assert(allclose(A[:2,:2], zeros([2,2])))\n",
    "assert(allclose(A[:2,2:4], identity(2)))\n",
    "assert(allclose(A[2:4,:2], stiffness(20.)))\n",
    "assert(allclose(A[2:4,2:4], damping(10.)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### For samples of the values $d \\in [0,100]$ and the fixed value $k=1000$, plot the four eigenvalues on the complex plane."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for d in linspace(0,100,200):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bonus question: there is a bifurcation in the diagram above. Can you find the bifurcation point programmatically?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequency Response Plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In technical applications there occurs often linear systems of the form\n",
    "\\\\[\n",
    "\\dot x(t) = A x(t) + B u(t)\n",
    "\\\\]\n",
    "where $u$ is an given input signal. $x$ is called the state. From the state some quantities $y(t)$ can be\n",
    "measured, this is decribed by the equation\n",
    "\\\\[\n",
    " y(t)=C x(t).\n",
    "\\\\]\n",
    "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. \n",
    "\n",
    "The 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",
    "\\\\]\n",
    "and $i$ is the imaginary unit.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "inv?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "norm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def amplitude(A, B, C, omega):\n",
    "    pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find out the relationship between $A$'s eigenvalues and the peak(s) in the figure.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "gist_id": "708b799a2d96af9bc96f",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.4.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}