Section_12_Scipy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scipy\n",
    "\n",
    "- [Linear Equations](#Linear-Equations)\n",
    "- [Non-linear Equations](#Non-linear-Equations)\n",
    "- [Integration - Single](#Integration---Single)\n",
    "- [Integration - Double](#Integration---Double)\n",
    "- [Optimisation](#Optimisation)\n",
    "- [Solve ODE](#Solve-ODE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Equations\n",
    "\n",
    "### Example - Solve Linear System of Equations\n",
    "\n",
    "- Solve the following system of linear equations:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    3x_1 + 6x_2 - 4x_3 &= 4 \\\\\n",
    "    2x_2 - x_3 &= 1 \\\\\n",
    "    x_1 + x_3 &= 0.25 \\\\\n",
    "\\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Rewrite into Linear Algebra Form:\n",
    "\n",
    "$$\n",
    "\\mathbf{A} \\cdot \\mathbf{x} = \\mathbf{b}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "        3 & 6 & -4 \\\\\n",
    "        0 & 2 & -1 \\\\\n",
    "        1 & 0 & 1 \\\\\n",
    "\\end{bmatrix} \\cdot \\begin{bmatrix}\n",
    "        x_1\\\\\n",
    "        x_2\\\\\n",
    "        x_3\\\\\n",
    "\\end{bmatrix}\n",
    " = \\begin{bmatrix}\n",
    "        4\\\\\n",
    "        1\\\\\n",
    "        0.25\\\\\n",
    "\\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Module `scipy` $\\to$ module `linalg` $\\to$ function `solve`\n",
    "    - `import scipy.linalg as la` $\\to$ `la.solve`\n",
    "    - Solves system of $N\\times N$ full rank equations\n",
    "- Test is solution is correct:\n",
    "    - Matrix multiplication `numpy.dot`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import scipy as sp\n",
    "\n",
    "from scipy.linalg import solve\n",
    "from scipy import linalg as la\n",
    "import scipy.linalg as la\n",
    "\n",
    "la.solve?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "\n",
    "A = np.array([\n",
    "    [3, 6, -4],\n",
    "    [0, 2, -1],\n",
    "    [1, 0,  1]\n",
    "])\n",
    "b = np.array([4, 1,0.25])\n",
    "\n",
    "x = la.solve(A, b)\n",
    "print(x)\n",
    "print(np.dot(A, x)-b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Fitting straight line through data points $\\to$ Least Squares\n",
    "\n",
    "- Consider a straight line $y = mx + c$\n",
    "\n",
    "\n",
    "- We want to fit it through the following points:\n",
    "    $$\n",
    "    \\begin{array}{lllll}\n",
    "        x & y \\\\ \\hline\n",
    "        0 & 2.1 \\\\\n",
    "        1 & 2.5 \\\\\n",
    "        2 & 2.9 \\\\\n",
    "        3 & 3.6 \\\\\n",
    "        4 & 3.9 \\\\\n",
    "    \\end{array}\n",
    "    $$\n",
    "\n",
    "\n",
    "- Which $m$ and $c$ will give the best fit?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Consider a straight line $y = mx+c$ and substitute each point:\n",
    "    $$\n",
    "    \\begin{align}\n",
    "        m\\times 0 + c &= 2.1\\\\\n",
    "        m\\times 1 + c &= 2.5\\\\\n",
    "        m\\times 2 + c &= 2.9\\\\\n",
    "        m\\times 3 + c &= 3.6\\\\\n",
    "        m\\times 4 + c &= 3.9\\\\\n",
    "    \\end{align}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Two unknowns and five equations. Cannot satisfy each equation - minimise error.\n",
    "    $$\n",
    "    \\begin{bmatrix}\n",
    "            0 & 1 \\\\\n",
    "            1 & 1 \\\\\n",
    "            2 & 1 \\\\\n",
    "            3 & 1 \\\\\n",
    "            4 & 1 \\\\\n",
    "    \\end{bmatrix}\n",
    "    \\cdot\n",
    "    \\begin{bmatrix}\n",
    "            m\\\\\n",
    "            c\\\\\n",
    "    \\end{bmatrix}\n",
    "     = \n",
    "    \\begin{bmatrix}\n",
    "            2.1\\\\\n",
    "            2.5\\\\\n",
    "            2.9\\\\\n",
    "            3.6\\\\\n",
    "            3.9\\\\\n",
    "    \\end{bmatrix}\n",
    "    $$\n",
    "\n",
    "\n",
    "- Cannot solve\n",
    "    $$\\mathbf{A}\\mathbf{x} = \\mathbf{b}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Least squares solution is given by:\n",
    "    $$\\mathbf{A}^{\\textrm{T}}\\mathbf{A}\\mathbf{x} = \\mathbf{A}^{\\textrm{T}}\\mathbf{b}$$\n",
    "\n",
    "\n",
    "- Python: Plot data and the least squares line fitted through the data.\n",
    "\n",
    "\n",
    "- LibreOffice: Plot data and insert linear trend line through the data.\n",
    "\n",
    "\n",
    "- Compare the coefficients?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Background on where does least squares come from?\n",
    "\n",
    "\n",
    "- The error $\\mathbf{E} = \\mathbf{A}\\mathbf{x}-\\mathbf{b}$ can be written as the error squares\n",
    "    $$\n",
    "    \\begin{align}\n",
    "        \\mathbf{E}^{\\textrm{T}}\\mathbf{E} &=\n",
    "        \\left(\n",
    "            \\mathbf{A}\\mathbf{x}-\\mathbf{b}\n",
    "        \\right)^{\\textrm{T}}\n",
    "        \\left(\n",
    "            \\mathbf{A}\\mathbf{x}-\\mathbf{b}\n",
    "        \\right) \\\\ &=\n",
    "        \\mathbf{x}^{\\textrm{T}}\\mathbf{A}^{\n",
    "            \\textrm{T}}\\mathbf{A}\\mathbf{x} - \n",
    "            2\\mathbf{A}^{\\textrm{T}\n",
    "        }\\mathbf{b}\\mathbf{x}\\\\\n",
    "    \\end{align}\n",
    "    $$\n",
    "\n",
    "\n",
    "- Minimum error\n",
    "    $$\n",
    "    \\frac{\n",
    "        \\mathrm{d} \\mathbf{E}^{\\textrm{T}}\\mathbf{E}\n",
    "    }{\n",
    "        \\mathrm{d} \\mathbf{x}\n",
    "    } = \n",
    "    2\\mathbf{A}^{\\textrm{T}}\\mathbf{A}\\mathbf{x} -\n",
    "    2\\mathbf{A}^{\\textrm{T}}\\mathbf{b} = \\mathbf{0}\n",
    "    $$\n",
    "    $$\n",
    "    \\begin{array}{lll}\n",
    "        \\mathbf{A}^{\\textrm{T}}\\mathbf{A}\\mathbf{x} = \\mathbf{A}^{\\textrm{T}}\\mathbf{b}\\\\\n",
    "    \\end{array}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np \n",
    "import scipy.linalg as la \n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "A = np.array([\n",
    "    [0, 1],\n",
    "    [1, 1],\n",
    "    [2, 1],\n",
    "    [3, 1],\n",
    "    [4, 1]\n",
    "])\n",
    "b = np.array([2.1, 2.5, 2.9, 3.6, 3.9])\n",
    "\n",
    "coefficients = la.solve(\n",
    "    np.dot(A.transpose(), A),\n",
    "    np.dot(A.transpose(), b)\n",
    ")\n",
    "\n",
    "x = np.linspace(0, 5, 20)\n",
    "y = coefficients[0]*x + coefficients[1]\n",
    "plt.plot(x, y, 'r-', label='fitted line')\n",
    "\n",
    "plt.plot(\n",
    "    [0, 1, 2, 3, 4],\n",
    "    [2.1, 2.5, 2.9, 3.6, 3.9],\n",
    "    'go', label='data'\n",
    ")\n",
    "\n",
    "plt.title('Line fitted through data points')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y')\n",
    "plt.legend(loc='upper left')\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-linear Equations\n",
    "\n",
    "- We cannot write a system of non-linear equations into a similar form as with a system of linear equations $\\mathbf{A}\\mathbf{x}=\\mathbf{b}$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We therefore write it into residual form:\n",
    "    $$\n",
    "    \\mathbf{R}(\\mathbf{x}) = \\begin{bmatrix}\n",
    "            R_0(\\mathbf{x})\\\\\n",
    "            R_1(\\mathbf{x})\\\\\n",
    "            \\vdots\\\\\n",
    "            R_n(\\mathbf{x})\\\\\n",
    "        \\end{bmatrix}\n",
    "     = \n",
    "     \\begin{bmatrix}\n",
    "            0\\\\\n",
    "            0\\\\\n",
    "            \\vdots\\\\\n",
    "            0\\\\\n",
    "     \\end{bmatrix} = \\mathbf{0}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For example:\n",
    "    $$\n",
    "    \\mathbf{R}(\\mathbf{x}) = \\begin{bmatrix}\n",
    "            0.8\\cos(\\frac{\\pi}{4}) + 2\\cos(x_0) - 2\\cos(x_1) - 3\\\\\n",
    "            0.8\\sin(\\frac{\\pi}{4}) + 2\\sin(x_0) - 2\\sin(x_1)\\\\\n",
    "    \\end{bmatrix} = \n",
    "    \\begin{bmatrix}\n",
    "            0\\\\\n",
    "            0\\\\\n",
    "    \\end{bmatrix}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Need to define a function for the residual $\\mathbf{R}({x})$.\n",
    "\n",
    "\n",
    "- Function takes array (unknown variables) as input and returns an array (result of each equation).\n",
    "\n",
    "\n",
    "- `from scipy.optimize import fsolve`\n",
    "    - First argument: **Function** that defines the residual.\n",
    "    - Second argument: initial guess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "\n",
    "def residual(x):\n",
    "    r1 = 0.8*np.cos(np.pi/4) + 2*np.cos(x[0]) - 2*np.cos(x[1]) - 3\n",
    "    r2 = 0.8*np.sin(np.pi/4) + 2*np.sin(x[0]) - 2*np.sin(x[1])\n",
    "    return r1, r2\n",
    "\n",
    "\n",
    "x0 = np.array([0, 0])\n",
    "ans = opt.fsolve(residual, x0)\n",
    "print(ans)\n",
    "print(residual(ans))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext nbtutor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%nbtutor -r -f --depth 3\n",
    "def twister():\n",
    "    print(\"Red Lorry, Yellow Lorry\")\n",
    "\n",
    "def repeat(func, N):\n",
    "    for i in range(N):\n",
    "        func()\n",
    "\n",
    "repeat(twister, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Kinematics (Four-bar linkage)\n",
    "\n",
    "- Four-bar linkage (kinematics) with $L_1$, $L_2$, $L_3$ and $L_4$ the length of the links.\n",
    "\n",
    "\n",
    "- Link 1 is fixed in a horizontal position.\n",
    "\n",
    "\n",
    "- The link lengths are $L_1=3$, $L_2=0.8$, $L_3=2$ and $L_4=2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/fourbar.svg\" alt=\"Four Bar Linkage\" style=\"height: 300px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- For a given $\\theta_1$, we can solve for $\\theta_2$ and $\\theta_3$ using the following two non-linear equations:\n",
    "    $$\n",
    "    \\begin{array}{l}\n",
    "        L_2\\cos(\\theta_1) + L_3\\cos(\\theta_2) + L_4\\cos(\\theta_3) - L_1 = 0 \\\\\n",
    "        L_2\\sin(\\theta_1) + L_3\\sin(\\theta_2) - L_4\\sin(\\theta_3) = 0 \\\\\n",
    "    \\end{array}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Solve for $\\theta_2$ and $\\theta_3$ for $\\theta_1=\\pi/4$\n",
    "\n",
    "\n",
    "- As a guess to the solution of the non-linear system use $[0.707; 0.707]$.\n",
    "\n",
    "\n",
    "- Plot the three links for the solved values of $\\theta_1$ and $\\theta_2$.\n",
    "\n",
    "\n",
    "- The coordinates of link three are given by\n",
    "    $$\n",
    "    \\begin{align}\n",
    "        x_1 &= L_2\\cos(\\theta_1) \\\\\n",
    "        y_1 &= L_2\\sin(\\theta_1) \\\\\n",
    "        x_2 &= L_2\\cos(\\theta_1) + L_3\\cos(\\theta_2) \\\\\n",
    "        y_2 &= L_2\\sin(\\theta_1) + L_3\\sin(\\theta_2)\n",
    "    \\end{align}\n",
    "    $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy import pi, cos, sin\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def fourbar(theta23):\n",
    "    L = [3.0, 0.8, 2.0, 2.0]\n",
    "    theta = [pi/4, theta23[0], theta23[1]]\n",
    "    R = [\n",
    "        L[1]*cos(theta[0]) + L[2]*cos(theta[1]) + L[3]*cos(theta[2]) - L[0],\n",
    "        L[1]*sin(theta[0]) + L[2]*sin(theta[1]) - L[3]*sin(theta[2])\n",
    "    ]\n",
    "    return R\n",
    "\n",
    "\n",
    "def plot_links(theta23):\n",
    "    L = [3.0, 0.8, 2.0, 2.0]\n",
    "    theta = [pi/4, theta23[0], theta23[1]]\n",
    "    plt.plot(\n",
    "        [0, L[1]*cos(theta[0]), L[1]*cos(theta[0]) + L[2]*cos(theta[1]), L[0], 0],\n",
    "        [0, L[1]*sin(theta[0]), L[1]*sin(theta[0]) + L[2]*sin(theta[1]),    0, 0],\n",
    "        'k-'\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.optimize as opt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x0 = [0.707, 0.707]\n",
    "\n",
    "theta23 = opt.fsolve(fourbar, x0)\n",
    "\n",
    "plot_links(theta23)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Works but every time we change a length or $\\theta_1$ we need to change it multiple places (2 cells)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Improve by defining additional arguments.\n",
    "    - Modify the function to take additional arguments that are known parameters of the problem e.g. $\\theta_1$ and lengths.\n",
    "    - `fsolve` pass a third argument of additional parameters.\n",
    "    - Define additional parameters in $()$ as third argument e.g. `(theta1,lengths)`.\n",
    "    - Order of additional parameters define the order of the parameters in the tuple argument.\n",
    "    - Additional arguments are available for many `scipy` functions.\n",
    "    - `fsolve(func, x0, args=(),...)`\n",
    "        - `args`: tuple $\\to$ Any extra arguments to `func`.\n",
    "    - Python mechanism to pass known information (parameters) to a function that is given as an arguments for another function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from numpy import pi, cos, sin\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def fourbar(theta23, L, theta1):\n",
    "    theta = [theta1, theta23[0], theta23[1]]\n",
    "    R = [\n",
    "        L[1]*cos(theta[0]) + L[2]*cos(theta[1]) + L[3]*cos(theta[2]) - L[0],\n",
    "        L[1]*sin(theta[0]) + L[2]*sin(theta[1]) - L[3]*sin(theta[2])\n",
    "    ]\n",
    "    return R\n",
    "\n",
    "\n",
    "def plot_links(theta23, L, theta1):\n",
    "    theta = [theta1, theta23[0], theta23[1]]\n",
    "    plt.plot(\n",
    "        [0, L[1]*cos(theta[0]), L[1]*cos(theta[0]) + L[2]*cos(theta[1]), L[0], 0],\n",
    "        [0, L[1]*sin(theta[0]), L[1]*sin(theta[0]) + L[2]*sin(theta[1]),    0, 0],\n",
    "        'k-'\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import scipy.optimize as opt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "theta1 = pi/4\n",
    "x0 = [0.707, 0.707]\n",
    "L = [3.0, 0.8, 2.0, 2.0]\n",
    "\n",
    "theta23 = opt.fsolve(fourbar, x0, args=(L, theta1))\n",
    "\n",
    "plot_links(theta23, L, theta1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Solve for $\\theta_2$ and $\\theta_3$ for $\\theta_1$ between 0 and $\\pi$ using 10 linearly spaced intervals.\n",
    "\n",
    "- Plot then link three for the various solved values of $\\theta_1$ and $\\theta_2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import scipy.optimize as opt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x0 = [0.707, 0.707]\n",
    "L = [3.0, 0.8, 2.0, 2.0]\n",
    "angles = np.linspace(0, np.pi, 10)\n",
    "\n",
    "for theta1 in angles:\n",
    "    theta23 = opt.fsolve(fourbar, x0, args=(L, theta1))\n",
    "    plot_links(theta23, L, theta1)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration - Single\n",
    "\n",
    "### Example - Simple Integral\n",
    "\n",
    "- Compute the following definite integral in python:\n",
    "    $$\n",
    "    \\int^{9}_{2} 4x^2 + 5 \\: \\mathrm{dx} =\n",
    "    \\left. \\left(\\frac{4}{3} x^3 + 5x\\right) \\right|_2^9 =\n",
    "    996.33\n",
    "    $$\n",
    "\n",
    "### Outcomes:\n",
    "-   Importing and using quad $\\to$ `from scipy import integrate` $\\to$ `integrate.quad`\n",
    "\n",
    "\n",
    "-   Passing a function as an input to another function\n",
    "\n",
    "\n",
    "-   Integrate a function from a to b\n",
    "\n",
    "\n",
    "-   Understanding the output from the `quad` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate as it\n",
    "\n",
    "def equation(x):\n",
    "    return 4*x**2 + 5\n",
    "\n",
    "ans, err = it.quad(equation, 2, 9)\n",
    "print(ans)\n",
    "print(err)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Single Integral\n",
    "\n",
    "- Compute the following integral in python:\n",
    "    $$\\int^{b}_{a} C e^{-x} \\: dx$$\n",
    "\n",
    "- First for $a = 2$, $b = 5$ and $C = 12$\n",
    "\n",
    "- Then for $a = 2$, $b = \\infty$ and $C = 22$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate as it\n",
    "\n",
    "def equation(x, C):\n",
    "    return C * np.exp(-x)\n",
    "\n",
    "for v in [12, 90, 23, -12, 9, 24]:\n",
    "    ans, err = it.quad(\n",
    "        equation, 0, np.inf,\n",
    "        args = v\n",
    "    )\n",
    "    print(ans)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outcomes:\n",
    "\n",
    "-   Passing a function as an input to another function\n",
    "\n",
    "\n",
    "-   Integrate a function from a to b\n",
    "\n",
    "\n",
    "-   Integrate a function with an infinite interval\n",
    "\n",
    "\n",
    "-   Creating an input function that has `args`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integration - Double\n",
    "\n",
    "### Example - Double Integral\n",
    "\n",
    "- Compute the following integral in python:\n",
    "    $$\\iint_D x y^2 \\: \\mathrm{dA}$$\n",
    "\n",
    "\n",
    "- Where $D$ is the rectangular region defined by $0 \\leq x \\leq 2$ and $0 \\leq y \\leq 1$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/double_integral_over_rectangle.png\" alt=\"Double Integral Rectangle\" style=\"height: 200px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "    \\int_0^1 \\left( \\int_0^2 xy^2 \\mathrm{dx}\\right) \\mathrm{dy} =\n",
    "    \\int_0^1 \\left(\\left. \\frac{x^2}{2} y^2 \\right|_{x=0}^{x=2} \\right) \\mathrm{dy} =\n",
    "    \\int_0^1 2 y^2 \\mathrm{dy} =\n",
    "    \\left. \\frac{2 y^3}{3} \\right|_{y=0}^{y=1} =\n",
    "    2/3\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "### Outcomes\n",
    "\n",
    "- Importing and using quad $\\to$ `from scipy import integrate` $\\to$ `integrate.dblquad`\n",
    "\n",
    "\n",
    "- Passing multiple functions as inputs to another function\n",
    "\n",
    "\n",
    "- Integrate a function for $x \\in [a, b]$ and $y \\in [c, d]$\n",
    "\n",
    "\n",
    "- Understanding the output from the `dblquad` function\n",
    "\n",
    "\n",
    "- Why must the inner integral limits $\\left( \\int^{2}_{x=0} \\right)$ be given as functions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import integrate as it\n",
    "\n",
    "def equation(y, x):\n",
    "    return x * y**2\n",
    "\n",
    "def inner_bound_lower(x):\n",
    "    return 0\n",
    "\n",
    "def inner_bound_upper(x):\n",
    "    return x/2\n",
    "\n",
    "ans, err = it.dblquad(\n",
    "    equation, 0, 2,\n",
    "    inner_bound_lower,\n",
    "    inner_bound_upper\n",
    ")\n",
    "print(ans)\n",
    "print(err)\n",
    "print(32/120)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Double Integral  (x First)\n",
    "\n",
    "- Compute the following integral in python:\n",
    "    $$\\iint_D x y^2 \\: \\mathrm{dA}$$\n",
    "\n",
    "\n",
    "- Where $D$ is the triangular region defined by $0 \\leq x \\leq 2$ and $0 \\leq y \\leq x/2$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/double_integral_over_triangle_x_first.png\" alt=\"Double Integral Triangle\" style=\"height: 200px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{aligned}\n",
    "    \\int_0^2 \\left(\\int_0^{x/2} xy^2 \\mathrm{dy} \\right) \\mathrm{dx} =\n",
    "    \\int_0^2 \\left(\\left.\\frac{x}{3} y^3 \\right|_{y=0}^{y=x/2}\\right) \\mathrm{dx} =\n",
    "    \\int_0^2 \\frac{x^4}{24} \\mathrm{dx} =\n",
    "    \\left. \\frac{x^5}{120} \\right|_{x=0}^{x=2} =\n",
    "    32/120\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "### Outcomes\n",
    "\n",
    "- Passing multiple functions as inputs to another function\n",
    "\n",
    "\n",
    "- Integrate a function for $x \\in [a, b]$ and $y \\in [c, d]$, where $c$ or $d$ is a function of $x$\n",
    "\n",
    "\n",
    "- Understanding the output from the `dblquad` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import integrate as it\n",
    "\n",
    "def function(y, x): \n",
    "    return x * y**2\n",
    "\n",
    "def inner_bound_lower(x): \n",
    "    return 0\n",
    "\n",
    "def inner_bound_upper(x): \n",
    "    return x / 2.0\n",
    "\n",
    "ans, err = it.dblquad(\n",
    "    function, 0, 2,\n",
    "    inner_bound_lower,\n",
    "    inner_bound_upper\n",
    ")\n",
    "\n",
    "print(ans)\n",
    "print(err)\n",
    "print(32/120)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example - Double Integral  (y First)\n",
    "\n",
    "- Compute the following integral in python:\n",
    "    $$\\iint_D x y^2 \\: \\mathrm{dA}$$\n",
    "\n",
    "\n",
    "- Where $D$ is the triangular region defined by $0 \\leq x \\leq 2$ and $0 \\leq y \\leq x/2$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/double_integral_over_triangle_y_first.png\" alt=\"Double Integral Triangle\" style=\"height: 200px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\begin{aligned}\n",
    "    \\int_0^1 \\left( \\int_{2y}^2 xy^2 \\mathrm{dx} \\right) \\mathrm{dy} =\n",
    "    \\int_0^1\\left(\\left.\\frac{x^2y^2}{2} \\right|_{x=2y}^{x=2}\\right) \\mathrm{dy} =\n",
    "    \\int_0^1 \\left( 2y^2 - 2y^4\\right) \\mathrm{dy} =\n",
    "    \\left. \\left(\\frac{2y^3}{3} - \\frac{2y^5}{5}\\right) \\right|_{y=0}^{y=1} = \n",
    "    2/3 - 2/5 = 4/15\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "\n",
    "### Outcomes\n",
    "\n",
    "- Passing multiple functions as inputs to another function\n",
    "\n",
    "\n",
    "- Integrate a function for $x \\in [a, b]$ and $y \\in [c, d]$, where $c$ or $d$ is a function of $x$\n",
    "\n",
    "\n",
    "- Understanding the output from the `dblquad` function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def function(x, y):\n",
    "    return x * y**2\n",
    "\n",
    "def inner_bound_lower(y): \n",
    "    return 2*y\n",
    "\n",
    "def inner_bound_upper(y): \n",
    "    return 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import integrate\n",
    "\n",
    "ans, error = integrate.dblquad(\n",
    "    function, 0, 1,\n",
    "    inner_bound_lower,\n",
    "    inner_bound_upper\n",
    ")\n",
    "print(ans)\n",
    "print(error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimisation\n",
    "\n",
    "### Example - Tin Can\n",
    "\n",
    "- A cylindrical-shaped tin-can must have a capacity of $1000cm^3$.\n",
    "\n",
    "\n",
    "- Determine the dimensions that require the minimum amount of tin for the can (assume no waste material)\n",
    "\n",
    "\n",
    "- The smallest can that the market will accept has a diameter of $6cm$ and a height of $4cm$\n",
    "\n",
    "\n",
    "- Variables:\n",
    "    $$x = [x_1, x_2] = [r, h]$$\n",
    "\n",
    "\n",
    "- Surface Area (Objective Function):\n",
    "    $$A(x) = 2 \\pi \\left({x_1}^2 + {x_1} {x_2} \\right)$$\n",
    "\n",
    "\n",
    "- Constraints (Inequality):\n",
    "    $$\\begin{align}\n",
    "    g_1(x) &= \\pi {x_1}^2 x_2 - 1000 \\geq 0 \\\\\n",
    "    g_2(x) &= x_1 - 3 \\geq 0 \\\\\n",
    "    g_3(x) &= x_2 - 4 \\geq 0 \n",
    "    \\end{align} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outcomes\n",
    "\n",
    "- Importing and using `fmin_slsqp` $\\to$ `from scipy.optimize import fmin_slsqp`\n",
    "\n",
    "\n",
    "- Creating the objective function and constraint functions\n",
    "\n",
    "\n",
    "- Passing multiple functions as inputs to another function\n",
    "\n",
    "\n",
    "- Understanding the output from the `fmin_slsqp` function\n",
    "\n",
    "\n",
    "- Non-linear optimisation problem\n",
    "\n",
    "\n",
    "- Need to write $g(x)$ in the form $g(x) >= 0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def function(x): \n",
    "    return np.pi * (x[0]**2 + x[0]*x[1])\n",
    "\n",
    "\n",
    "def ineq_constraints(x): \n",
    "    return [\n",
    "        np.pi * (x[0]**2) * x[1]- 1000,\n",
    "        x[0] - 3, \n",
    "        x[1] - 4\n",
    "    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import optimize\n",
    "\n",
    "x = optimize.fmin_slsqp(\n",
    "    function,\n",
    "    [1, 1],\n",
    "    f_ieqcons=ineq_constraints\n",
    ")\n",
    "\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve ODE\n",
    "\n",
    "### Example - Spring, mass and damper system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"./figures/springdamper.svg\" alt=\"Spring Damper\" style=\"height: 220px;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The position of the mass $m$ is given by $x$\n",
    "\n",
    "\n",
    "- There is some force $f(t)$ applied to the mass that can vary with time\n",
    "\n",
    "\n",
    "- $x$ is a function of time $t$ therefore $x(t)$\n",
    "\n",
    "\n",
    "- Determine the position $x(t)$ of the mass $m$ for various values of $t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Outcomes\n",
    "\n",
    "- Understanding and using the `odeint` function\n",
    "\n",
    "\n",
    "- Understanding `tuple` unpacking"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mathematics\n",
    "\n",
    "- Spring force:\n",
    "\n",
    "$$F_s(t) = -kx(t)$$\n",
    "\n",
    "\n",
    "- Damping force:\n",
    "$$F_d(t) = -cv(t) = -c\\dot{x}(t)$$\n",
    "\n",
    "\n",
    "- Force balance (Newtons law)\n",
    "$$F = m a$$\n",
    "\n",
    "$$F_s(t) + F_d(t) + f(t) = ma(t)$$\n",
    "\n",
    "$$-k x(t) - c v(t) + f(t) = ma(t)$$\n",
    "\n",
    "$$f(t) = k x(t) + c v(t) + ma(t)$$\n",
    "\n",
    "\n",
    "- Note that\n",
    "    - $v(t) = \\dot x(t)$\n",
    "    - $a(t) = \\dot v(t) = \\ddot x(t)$\n",
    "\n",
    "$$f(t) = k x(t) + c \\dot x(t) + m \\ddot x(t)$$\n",
    "\n",
    "\n",
    "- second order differential equation !! How do we solve this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Need to write this second order DE as two first order DE's and solve simultaneously\n",
    "    - Fist equation:\n",
    "        $$v(t) = \\dot x(t)$$\n",
    "    - Second equation:\n",
    "        $$f(t) = k x(t) + c v(t) + m \\dot v(t)$$\n",
    "\n",
    "\n",
    "- Re-arranging for the first time derivatives:\n",
    "    - Fist equation:\n",
    "        $$\\dot x(t) = v(t)$$\n",
    "    - Second equation:\n",
    "        $$\\dot v(t) = \\frac{f(t)}{m} - \\frac{k x(t)}{m} - \\frac{c v(t)}{m} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "\n",
    "def force(t):\n",
    "    return 0\n",
    "\n",
    "\n",
    "def func(variables, t, c, k, m):\n",
    "    xt, vt = variables\n",
    "    ft = force(t)\n",
    "    return [\n",
    "        vt,\n",
    "        ft/m - k*xt/m - c*vt/m\n",
    "    ]\n",
    "\n",
    "\n",
    "x_init = 10  # initial position\n",
    "v_init = 0  # initial velocity\n",
    "\n",
    "times = np.linspace(1, 100, 200)\n",
    "results = odeint(func, [x_init, v_init], times, args=(0.2, 1, 1))\n",
    "\n",
    "plt.plot(times, results[:, 0], 'r')\n",
    "plt.plot(times, results[:, 1], 'b')\n",
    "plt.legend(['position', 'velocity'])\n",
    "plt.xlabel('time')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}