{
 "metadata": {
  "name": "",
  "signature": "sha256:d5f6b5de3f2868f5c0c7e04d3738d07a4012b4fde4e45c7ecd30527bdaaa5137"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from __future__ import division\n",
      "\n",
      "import numpy as np\n",
      "import matplotlib.pyplot as pt"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def rk4_step(y, t, h, f):\n",
      "    k1 = f(t, y)\n",
      "    k2 = f(t+h/2, y + h/2*k1)\n",
      "    k3 = f(t+h/2, y + h/2*k2)\n",
      "    k4 = f(t+h, y + h*k3)\n",
      "    return y + h/6*(k1 + 2*k2 + 2*k3 + k4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mesh = np.linspace(0, 1, 200)\n",
      "dx = mesh[1]-mesh[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 280
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "(all of the PDEs below use *periodic* boundary conditions: $u(0)=u(1)$)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Advection equation:** $u_t+u_x=0$\n",
      "\n",
      "Equivalent: $u_t=-u_x$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f_advection(t, u):\n",
      "    du = (np.roll(u, -1, axis=-1) - np.roll(u, 1, axis=-1))/(2*dx)\n",
      "    return -du"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 252
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Heat equation:** $u_t=u_{xx}$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f_heat(t, u):\n",
      "    d2u = (\n",
      "        np.roll(u, -1, axis=-1)\n",
      "        - 2*u\n",
      "        + np.roll(u, 1, axis=-1))/(dx**2)\n",
      "    return d2u"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 307
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Wave equation:** $u_{tt}=u_{xx}$\n",
      "\n",
      "NOTE: Two time derivatives $\\rightarrow$ convert to first order ODE.\n",
      "\n",
      "$$u_t=v$$\n",
      "$$v_t=u_{xx}$$\n"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def f_wave(t, w):\n",
      "    u, v = w\n",
      "    d2u = (\n",
      "        np.roll(u, -1, axis=-1)\n",
      "        - 2*u\n",
      "        + np.roll(u, 1, axis=-1))/(dx**2)\n",
      "    return np.array([v, d2u])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 434
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Initial condition**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "current_t = 0\n",
      "\n",
      "#current_u = np.sin(2*np.pi*mesh)*0.5+1\n",
      "#current_u = (mesh > 0.3) & (mesh < 0.7)\n",
      "#current_u = (mesh > 0.45) & (mesh < 0.55)\n",
      "current_u = np.exp(-(mesh-0.5)**2*150)\n",
      "#current_u = 2*np.abs(mesh-0.5)\n",
      "\n",
      "current_u = np.array([current_u], dtype=np.float64)\n",
      "current_u.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 455,
       "text": [
        "(1, 200)"
       ]
      }
     ],
     "prompt_number": 455
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Add a second component if needed (for wave equation)\n",
      "current_u = np.vstack([current_u,np.zeros(len(mesh))])\n",
      "current_u.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 456,
       "text": [
        "(2, 200)"
       ]
      }
     ],
     "prompt_number": 456
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**Time loop**"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Run this cell many times in place (using Ctrl-Enter)\n",
      "\n",
      "dt = dx # experiment with this\n",
      "\n",
      "#current_f = f_advection\n",
      "#current_f = f_heat\n",
      "current_f = f_wave\n",
      "\n",
      "for i in range(5): # takes this many time steps at a time\n",
      "    current_u = rk4_step(current_u, current_t, dt, current_f)\n",
      "    current_t += dt\n",
      "\n",
      "pt.ylim([-0.25, 1.25])\n",
      "pt.grid()\n",
      "pt.plot(mesh, current_u[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 507,
       "text": [
        "[<matplotlib.lines.Line2D at 0x34dbff50>]"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD7CAYAAABt0P8jAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XtYVPW6B/AvBp68tMVbmgyFAgmGzLAlSc3CegzU1LZZ\nG592R9DILDOsk5bu0jxHE7tpUjt2aZi3bLfdYqGzz8mYLEXxEmKiHneBAiZFQl6PcvmdP37BEkUd\nZs2sy8z38zzz6ILFWi8vM++seddv/ZafEEKAiIh8Riu9AyAiIm2x8BMR+RgWfiIiH8PCT0TkY1j4\niYh8DAs/EZGP8dc7AACw2WzYu3ev3mEQEZmK1WpFQUFBi3/OEEf8e/fuhRCCDyEwe/Zs3WMwyoO5\nYC6Yi6s/XD1gNkThJ0VJSYneIRgGc6FgLhTMhXqqC/+ECRPQrVs39O3bt9nvr1q1ClarFdHR0Rg0\naBAKCwvV7pKIiFRQXfhTUlJgt9uv+P1evXphy5YtKCwsxEsvvYTHH39c7S69WnJyst4hGAZzoWAu\nFMyFen5CCNVz9ZSUlGDkyJHYt2/fVderqqpC3759UVZW1jQIPz+4IQwiIp/iau3UtMe/dOlSDB8+\nXMtdmo7D4dA7BMNgLhTMhYK5UE+z4Zy5ublYtmwZtm7d2uz3k5OTERISAgAIDAyEzWZDfHw8AOUP\nzWXfWm5glHj0XC4oKDBUPHouNwxfNEo8Wi47HA5kZWUBQGO9dIUmrZ7CwkKMGTMGdrsdYWFhlwfB\nVg8RUYsZttVz9OhRjBkzBitXrmy26BMRkbZUF/5x48Zh4MCBOHToEIKDg7Fs2TJkZmYiMzMTADB3\n7lxUVVVh8uTJiImJQf/+/VUH7c0ubXP4MuZCwVwomAv1VPf416xZc9Xvf/DBB/jggw/U7oaIiNzE\nLT1+1UGwx09E1GKG7fETEZGxsPAbDPuXCuZCwVwomAv1WPiJiHwMe/xERCbFHj8RETmFhd9g2L9U\nMBcK5kLBXKjHwk9E5GPY4yciMin2+ImIyCks/AbD/qWCuVAwFwrmQj0WfiIiH8MePxGRSbHHT0RE\nTmHhNxj2LxXMhYK5UDAX6rHwExH5GPb4iYhMij1+IiJyCgu/wbB/qWAuFMyFgrlQj4WfiMjHqOrx\nT5gwATk5Objxxhuxb9++ZteZOnUqNm3ahLZt2yIrKwsxMTGXB8EePxFRi+nS409JSYHdbr/i9zdu\n3Ih//etfOHz4MP76179i8uTJanZHRERuoKrwDx48GB07drzi9zds2IDx48cDAOLi4lBdXY2Kigo1\nu/R67F8qmAsFc6FgLtTzaI+/vLwcwcHBjcsWiwVlZWWe3CUREV2Dx0/uXtp/8vPz8/QuTS0+Pl7v\nEAyDuVAwFwrmQj1/T248KCgIpaWljctlZWUICgpqdt3k5GSEhIQAAAIDA2Gz2Rr/wA0f7bjMZS5z\n2ZeXHQ4HsrKyAKCxXrpEqFRcXCyioqKa/V5OTo4YNmyYEEKIvLw8ERcX1+x6bgjDa+Tm5uodgmEw\nFwrmQsFcKFytnaqO+MeNG4evvvoKlZWVCA4OxiuvvIKamhoAwKRJkzB8+HBs3LgRYWFhaNeuHT78\n8EM1uyMiIjfgXD1ERCbFuXqIiMgpLPwG03Aih5iLizEXCuZCPRZ+IiIfwx4/EZFJscdPREROYeE3\nGPYvFcyFgrlQMBfqsfATEfkY9viJiEyKPX4iInIKC7/BsH+pYC4UzIWCuVCPhZ+IyMewx09EZFLs\n8RMRkVNY+A2G/UsFc6FgLhTMhXos/EREPoY9fiIik2KPn4iInMLCbzDsXyqYCwVzoWAu1GPhJyLy\nMezxExGZFHv8RETkFNWF3263IyIiAuHh4UhPT7/s+5WVlUhMTITNZkNUVBSysrLU7tKrsX+pYC4U\nzIWCuVBPVeGvq6vDlClTYLfbUVRUhDVr1uDAgQNN1snIyEBMTAwKCgrgcDjw3HPPoba2VlXQRETk\nOlWFPz8/H2FhYQgJCUFAQACSkpKQnZ3dZJ2bbroJJ0+eBACcPHkSnTt3hr+/v5rderX4+Hi9QzAM\n5kLBXCiYC/VUVeDy8nIEBwc3LlssFuzYsaPJOqmpqbjnnnvQo0cPnDp1Cp988omaXRIRkUqqjvj9\n/Pyuuc78+fNhs9lw7NgxFBQU4KmnnsKpU6fU7NarsX+pYC4UzIWCuVBP1RF/UFAQSktLG5dLS0th\nsViarLNt2zbMmjULABAaGoqePXvi0KFDiI2NbbJecnIyQkJCAACBgYGw2WyNH+ka/tBc9q3lBkaJ\nR8/lgoICQ8Wj53JBQYGh4tFy2eFwNA6QaaiXrlA1jr+2tha9e/fG5s2b0aNHD/Tv3x9r1qxBZGRk\n4zrPPvssOnTogNmzZ6OiogL9+vVDYWEhOnXqpATBcfxERC3mau1UdcTv7++PjIwMJCQkoK6uDhMn\nTkRkZCQyMzMBAJMmTcLMmTORkpICq9WK+vp6LFy4sEnRJyIibfHKXYNxOByNH/F8HXOhYC4UzIWC\nV+4SEZFTeMRPRGRSPOInIiKnsPAbzKVDGX0Zc6FgLhTMhXos/EREPoY9fiIik2KPn4iInMLCbzDs\nXyqYCwVzoWAu1GPhJyLyMezxExGZFHv8RETkFBZ+g2H/UsFcKJgLBXOhHgs/EZGPYY+fiMik2OMn\nIiKnsPAbDPuXCuZCwVwomAv1WPiJiHwMe/xERCbFHj8RETmFhd9g2L9UMBcK5kLBXKjHwk9E5GNU\n9/jtdjvS0tJQV1eHxx57DDNmzLhsHYfDgWnTpqGmpgZdunS57B2bPX4iopZztXaqKvx1dXXo3bs3\nvvjiCwQFBeH222/HmjVrEBkZ2bhOdXU1Bg0ahH/+85+wWCyorKxEly5d3BI8EZEv0+Xkbn5+PsLC\nwhASEoKAgAAkJSUhOzu7yTqrV6/Ggw8+CIvFAgCXFX1qiv1LBXOhYC4UzIV6qgp/eXk5goODG5ct\nFgvKy8ubrHP48GGcOHECQ4YMQWxsLFasWKFml0REpJK/mh/28/O75jo1NTXYs2cPNm/ejLNnz2LA\ngAG44447EB4e3mS95ORkhISEAAACAwNhs9kQHx8PQHmH94Xl+Ph4Q8XDZeMsNzBKPHotN3zNKPFo\nuexwOJCVlQUAjfXSFap6/Nu3b8ecOXNgt9sBAK+++ipatWrV5ARveno6zp07hzlz5gAAHnvsMSQm\nJmLs2LFKEOzxExG1mC49/tjYWBw+fBglJSW4cOEC1q5di1GjRjVZZ/To0fjmm29QV1eHs2fPYseO\nHejTp4+a3Xq1S4/ufBlzoWAuFMyFeqpaPf7+/sjIyEBCQgLq6uowceJEREZGIjMzEwAwadIkRERE\nIDExEdHR0WjVqhVSU1NZ+ImIdMS5eoiITIpz9RARkVNY+A2G/UsFc6FgLhTMhXos/EREPoY9fiIi\nk2KPn4iInMLCbzDsXyqYCwVzoWAu1GPhJyLyMezxExGZFHv8RETkFBZ+g7la//LoUSAvD9izB6ir\n0y4mvbCXq/CFXJw+DXzzDbB7N1BVdeX1fCEXnsbCbwLr1wMDBwKxsUBaGvCnPwFdugDPPw+cPKl3\ndETqHDwIjBwJdO8OPPcc8NhjQM+ewLhxQFGR3tF5JxZ+g7l4zvEzZ4CJE4EZM4Dp04Fjx4AdO+SL\nobAQqKwEIiOBnTv1i9eTLs6Fr/PWXCxZAgweDAwZAvz0k3x+f/stUFwM9OsH3H038MEHTX/GW3Oh\nJZ7cNaizZ4Fhw4AePYD33wfat29+vexseYS0ahVw333axkjkKiGAmTOBf/wDsNuBK91T5MABYMwY\n4JFHgD//WdMQTYEnd72Ew+FATQ3whz8At9wiC/qVij4AjB4tXzyPPCKPlLwJe7kKb8vFm28CGzfK\nnv7VbiQVGQnk5gIffQS89Zb8mrflQg+q5uMnz5gxA/D3B5YtA1o58dZ8553Au+/KN4EdO4CbbvJ8\njESuysmRhT8vT56rupbu3YHNm4G4OCA6GrjuOs/H6O3Y6jGYdevkCa7du4FOnVr2sy+9BOzaJY+k\nnLgdMpHmfvwRsNnkgIUBA1r2s19+KT/Z7t4tW6Dkeu1k4TeQn34C+vYFPvsM6N+/5T9fUyNfTJMm\nAamp7o+PSA0hgFGjgJgYYO5c17YxZ44czPD55zy4Adjj9wppacCQIQ6Xij4ABATIXujMmUBZmXtj\n0wN7uQpvyMWqVfJ5qeYk7cyZQFGRAx9/7L64fBELv0HY7bI/n5ysbjt9+gBPPCHPExAZxalT8jn5\n3ntA69aub6d1a3n9yrPPAidOuC8+X8NWjwHU1MiTVq+/DowYoX57Z87I0RCrV8sTv0R6mzEDqKgA\nsrLcs70nngDatpUniX0Ze/wm9s478mTXf/+3+/qWH38s30h27mQvlPRVXCyvOv/uO/eNOKuoAG67\nDdi+HQgLc882zUi3Hr/dbkdERATCw8ORnp5+xfV27twJf39/rFu3Tu0uvcqvv8oTXW+8IQu0u3q5\nDz8M1NfLMf5m5Q19bXcxcy5eeQV46in3FX2Hw4Fu3eToN7Y0XaOq8NfV1WHKlCmw2+0oKirCmjVr\ncODAgWbXmzFjBhITE336yL45ixYBiYmy1eNOrVoB//VfcoinL0zoRsZ04IAcXvzcc+7fdlqavBag\noMD92/Z2qgp/fn4+wsLCEBISgoCAACQlJSE7O/uy9ZYsWYKxY8eia9euanbndU6ckHOVvPyy8jV3\nzkMybBjQsSNMOwKCc7IozJqLOXNk0e/QwX3bbMhFmzbAf/yHPMChllFV+MvLyxEcHNy4bLFYUF5e\nftk62dnZmDx5MgDZkyLpjTeABx4AQkM9s30/P3nEv2CBHENNpKWDB+V0C0895bl9TJokp3347jvP\n7cMbqZqywZkinpaWhgULFjSehLhSqyc5ORkhv03aERgYCJvN1vjO3tDf9KblM2eA996Lx65dTb9/\ncS/XHfu77z7g3DkHFiwAXnzROL+/M8sNXzNKPHouFxQUIC0tzTDxOLOclRWPp58Gdu1y7/YXLVrU\nWB/atQNGj3Zg6lTgyy+N9ft7YtnhcCDrt6FRIVeb5OhahAp5eXkiISGhcXn+/PliwYIFTdbp2bOn\nCAkJESEhIaJ9+/bixhtvFNnZ2U3WURmGKb35phAPP3z513Nzc92+r9WrhbjrLrdv1uM8kQuzMlsu\niouF6NRJiBMn3L/tS3Nx6pQQXbsKceCA+/dldK7WTlXDOWtra9G7d29s3rwZPXr0QP/+/bFmzRpE\nRkY2u35KSgpGjhyJMWPGNPm6rw3nrK2VQ9DWrpUTT2mxv/BwYM0a4I47PL8/orQ0ebHVwoXa7G/e\nPNlaWrFCm/0Zhau1U1Wrx9/fHxkZGUhISEBdXR0mTpyIyMhIZGZmAgAmTZqkZvNea906wGLRpugD\ncqbP554D0tPNPbyTzOHUKTl1iJajbaZMkQdT33/vuXNm3oQXcGlMCHnUPWOGvMHEpRwOR2Nvz53O\nnpXznm/ZAkREuH3zHuGpXJiRmXKxZIl8nv3tb57Z/pVyMXOmvG/v2297Zr9GxEnaTCIvT94ycfRo\nbffbtq0cXfHaa9rul3xLfb0s/FOnar/vp54CVq4Eqqu137fZ8IhfYw8+CMTHA08/rf2+Kytlr//Q\nIeDGG7XfP3m/jRuBWbOAPXv0mSrkT38CrFY5kZsv4Fw9JvD997KvX1Jy9dspelJqKnDzzXJ8P5G7\nJSYCSUnqZ5l11e7d8ralP/wgz215O7Z6TGDxYnlj9KsV/YvHsHvC1KnAX/4CXLjg0d24hadzYSZm\nyMXBg/K+z0lJnt3P1XLRrx/Qsyfw9797NgazY+HXSFWV7D9OmaJvHH37yimbPXXijXzXkiXA448D\n11+vbxzTpsnpmn2gieAytno0snAhsG+fMcYZb9gg5zfZsYNTNpN7VFfLI+39+/W/H25dHXDrrfK1\nNnCgvrF4Gls9Bnbhghxi9uyzekcijRgB/PKLnMucyB2WLZOTAupd9AHguuvkBWRvvaV3JMbFwq+B\nv/1NHoHExFx7XS16udddJ0cVLV7s8V2pYoa+tlaMnIu6OiAjQ7shnM7kIiVFThBXXOz5eMyIhd/D\nhJCzcBrlaL9BSoq845c33JSd9PX550DXrsaaDqR9e/kcX7JE70iMiT1+D3M45P1Bi4rkzVGMZOpU\n+QKZP1/vSMjM7r0XmDABeOQRvSNp6uhR+Sm7uBj43e/0jsYzOI7foEaNkj11I05bdPgwMGgQcOSI\nvKkFUUvt2wckJMhrU1q31juayyUlyWtnpk3TOxLP4MldAzp0SJ5A/fd/d/5ntOzlhocD/fsDq1Zp\ntssWMXJfW2tGzcWSJfITrZZFvyW5ePZZeS6rttZz8ZgRC78HLVokXxRGPpp+5hn5wvDSD1zkQb/8\nIgcuGPHTbIP+/YHgYM5Keym2ejykYV6cAweA7t31jubKhACiouSR2z336B0NmUl6ujx3tXy53pFc\n3T/+Ia+jycvTOxL3Y6vHYN57T067bOSiD8gLuKZONf7QTjKW2lrgnXf0mYWzpUaNAn76Cdi2Te9I\njIOF3wPOn5cvCldOKOnRy330Ufmi+P57zXd9VUbta+vBaLlYv15O9tevn/b7bmkuGi7oevNNz8Rj\nRiz8HrB6tZwaNipK70ic07atHI6XkaF3JGQWixeb42i/QUqKHFr9ww96R2IM7PG7mRBAdLS8aOu+\n+/SOxnkNY55LSoAbbtA7GjKynTuBsWPlJ0QzTX08Y4b8NL5okd6RuA97/AaxaZP8aDl0qN6RtMzN\nN8uTu1lZekdCRvfmm3I0mJmKPiCnKfnoI96hC2Dhd7sFC+SRhauzXurZy33mGTm6p75etxCaMFpf\nW09GycXRo3Kqj4kT9YvB1VxYLMDw4cD777s3HjNSXfjtdjsiIiIQHh6O9PT0y76/atUqWK1WREdH\nY9CgQSgsLFS7S8Patk3OffPQQ3pH4ppBg2SbZ9MmvSMho3r7bXl3rQ4d9I7ENdOmyd+hpkbvSHQm\nVKitrRWhoaGiuLhYXLhwQVitVlFUVNRknW3btonq6mohhBCbNm0ScXFxl21HZRiGMWqUEO++q3cU\n6ixfLsTQoXpHQUb0669CdOokREmJ3pGoEx8vxEcf6R2Fe7haO1Ud8efn5yMsLAwhISEICAhAUlIS\nsrOzm6wzYMAAdPjt8CAuLg5lXjod5P798sYmet1r1F3++Ec5/0pRkd6RkNEsXSrPXd1yi96RqDNr\nFjBvnpxO2lepKvzl5eUIDg5uXLZYLCgvL7/i+kuXLsXw4cPV7NKwXntNDm9TOz2D3r3cf/s3eQn+\n22/rGgYA/XNhJHrnorZWDuE0wvTianNx771Ax47Ap5+6Jx4zUnVe3q8FZzBzc3OxbNkybN26Vc0u\nDenoUeCzz4x3AZSrnngC6NMHeOUVoFs3vaMhI1i7Vo786t9f70jU8/MDXn4ZeP55eT7OaNOla0FV\n4Q8KCkJpaWnjcmlpKSwWy2XrFRYWIjU1FXa7HR07dmx2W8nJyQgJCQEABAYGwmazIT4+HoDyDm/U\n5SefdCAhAQgMVL+9+Ph43X+fgwcduPtu4LXX4vH66/rnl8tyuYHW+9+82YEXXwSWLtX3929Ybvia\nmu1dfz3Qpk081q8HOnXS9/dpybLD4UDWb2OuG+qlS9ScWKipqRG9evUSxcXF4vz5882e3D1y5IgI\nDQ0VeXl5V9yOyjB0dfiwEJ07C/HLL3pH4l6lpUJ07ChERYXekZDeVq8WYsAAIerr9Y7EvbKzhbBa\nzf17uVo7VX3I8ff3R0ZGBhISEtCnTx/88Y9/RGRkJDIzM5GZmQkAmDt3LqqqqjB58mTExMSgvzd8\nVrzI7NlyHpBOndyzvUuP7vRisQDjxskrkPVilFwYgV65qKsD/vM/5fPc1WtT3M1duRg5Uv772Wdu\n2Zy5uPkNyCUGCaPF9u0Tols3IU6edN82c3Nz3bcxlY4elUf9P/+sz/6NlAu96ZWLtWuFiIsz1lGx\nO3Oxbp0Q0dFC1Na6bZOacrV2cq4eFf7wB2DwYGOMdPCUJ56QIyBefVXvSEhr9fVy3qmFC+UVr95I\nCPkanjBBPsyG99zVWH6+nG//8GFj32FLrSNHgN//Xt5GsksXvaMhLX36qbzZSn6+cdo8nrBjh3wt\nHzoEtG+vdzQtw0naNCSEnI/npZfcX/SN1te+5Rbg4YeB+fO137fRcqEnrXNRUwP8+c/A3LnGK/ru\nzkVcHHDXXcDrr7t1s4bGwu+Cdevk/Ub1nKhKS7Nny9vrcS5z3/HBB0BQEJCYqHck2nj1VTlB4bFj\nekeiDbZ6WujcOXlx09KlvnWP2rlz5TQOH3+sdyTkaadOAbfeCuTkyDafr5gxQ94re+lSvSNxHnv8\nGnn5ZVkAfe1y7zNngIgIYNUq+bGYvNf06UBFhfFvou5uv/4qn+Pr18v2jxmwx6+BoiLgL3/x7I3J\njdrXbtdO3oDjqae0m9LWqLnQg1a5KCoCPvxQjuQxKk/lokMH+RxPTfX+aZtZ+J1UXw88/ricvyYo\nSO9o9DF2LNC9O+/N662EAKZMkYMWfHWOpqQkoEcPfS9c1AJbPU564w0gO1vesLmVD79dHj4MDBgg\nh8CFhuodDbnT++8DmZnA9u3mu62iO5WUALffDmzeLK9jMDL2+D1o3z55InfnTkDNvEje4o03gM8/\nly8MX34T9CalpfJEbm4uEBWldzT6y8qSbZ+dO+VU5UbFHr+HnDkj56x57TVtir4Z+tppacD5856f\ns98MudCKJ3NRXy+vWn3mGXMUfS2eF+PHA+Hhcupmb8TCfxVCyL7+7bfLJwJJ110HrFwpL+ravVvv\naEitBQvkG/kLL+gdiXH4+clhnTk53jmEma2eq1i0SI5wyMsD2rbVOxrj+eQTYOZM+XH4CrdZIIP7\n6it5u81du+SMrNRUQYG83eQXXwBWq97RXI6tHjdbv162dzZsYNG/kocfllPbPvyw9w9/80Y//CCL\n/ooVLPpXYrMB774L3H8/4E23C2fhb8aWLXIsb3a29jeWNltf+7XX5AiQp5+WrTF3MlsuPMnduThx\nAhg1Ss7HM3SoWzftcVo/Lx56SD6/hw2TV/Z6Axb+S2zdKserr10LxMbqHY3x+fvLXO3aBcyapXc0\n5IxTp+Q0y4mJ8oI8urbnn5dH/UOHyjdNs2OP/yIbN8qTuKtWAffdp3c05lJZCcTHAw88IO/YZLQZ\nHUn69VfZnouIkGP2+XdynhDyBHhODrBpExAcrHdE7PGrIgTw3ntySNtnn7Hou6JLFzkG3G4HnnwS\nqK3VOyK61PHj8s3ZapXPdxb9lvHzk/cnSEkBBg2Sn3LNyucL/+nTQHIy8M47wNdfA3fcoW88Zu5r\nd+0KfPmlPGmYmKi+H2rmXLib2lzk5cnW5dix8voLM194p/fz4rnngLfekj3/zEz3n9vSgon//Orl\n5spLslu1kpeph4frHZH5/e53smUWGwvExAD//KfeEfm22lpg3jxg9Gh5lD9rFo/03eHBB+WB4nvv\nASNGyCufzcQne/z/+hfw4ovylnLvvCNP2pD7bd4sb1YzeLC80QWHDGpryxZ5lXXXrvLGKkboSXub\nmhr53F68GJg6Vd5/+4YbtNu/bj1+u92OiIgIhIeHIz09vdl1pk6divDwcFitVnz77bdqd+kSIYBt\n24BHH5XtnOho4OBBFn1Puvde4Lvv5JBYqxWYNs27xkIbUcPz/P775XN9+nR53oVF3zMCAuQ9Onbv\nlvfsDQ0F5syR96o2NKFCbW2tCA0NFcXFxeLChQvCarWKoqKiJuvk5OSIYcOGCSGE2L59u4iLi7ts\nOyrDuKL6eiEOHRJi3jwhwsOF6N1biPR0IaqrPbI7t8jNzdU7BI8oKxPi2WeFCAwUYvRoIT755Np/\nB2/NhSuulYvyciHeflsIm02Inj2FePddIc6d0yY2rRn5eXHwoBBPPilE585C3HOPEMuXC1FZ6bn9\nuVo7VU2+mp+fj7CwMIT8NntZUlISsrOzERkZ2bjOhg0bMP63iW7i4uJQXV2NiooKdPPAhN+nT8sb\nSRQWyo+5X34pv37//cBHH8m76rC/qY+gIDmr55w5ctz/smWyDRQbCyQkAP36yXMCnTvrHanxCSHv\nDZuXJ/vMX38tpxIeMULeQOXee8198tbMeveW7eM33pAjBFeskPc46NVLzvA7cCBw223yfKKeU1+r\n2nV5eTmCL/oMabFYsGPHjmuuU1ZWdlnhr6uTJ6IaHg3L58/LGTJPn1YeVVXAjz82fRw5Avz0k0x8\nVBRw553yhhJhYeYq9vHx8XqH4FE33AA89ph8nDkj35y/+EKeEC4okCeHw8KAm28GbrklHkVFQGCg\nnAsoMFA+brgBaN266eO668z1d25QWwv83//Jx7lzyv9PnwZ+/vniRzz++lfgf/9XPtq0Afr3l8/z\njAz5Btq6td6/jTbM8Bq5/np5xe9DD8nzALt2yef6ypWy/VleLp/nwcHy3JfFIg96AgPlncAaHm3b\nyr9rQIDyXA8IkA81b+6qCr+fk680ccnJh+Z+rnVr+Q7o7y9fxA3/b90aaN9eebRrJxNy000yWbff\nLv8fHAz07Cl/lsyhXTt5MdHIkXK5vl4euRYXA0ePyjfzffuA6uqmj1OngAsXmj7q65u+ATQ8xRr+\n78wDkEfTDU/X5v5/te+5ul6bNrJQXH+98v+2beVJ2YZH9+5y3phnnpE3QuekeOYRECBvXjRggPK1\ns2flTY1KS+V5r7IyYP9++fz+9Vfl33Pn5BvHhQvKvw3/VzMeRlXhDwoKQulF45hKS0thuWToxqXr\nlJWVIaiZexc++mhyY8soMDAQNput8Z29Ydzu1ZbLyoCwMOfXN+ryxWOUjRCP1su9egFHjzrQs6d8\nI784J1f7+bo6YODAeNTVAVu2yO/feWc8hAC+/toBIYBBg+TyN9/I5YED5fLWrcr3AWDbNgf8/OSy\nn59cbtien59c389PWf7mG7k8eLBc/vpruf7ddze/3BBffHw8/P2dy09BQQEefTRN97+PEZYXLVrU\n4vpg1GVjrjyqAAAH9klEQVSrVS5HRDhfH7KysgAAt9wSgrlz4Ro1JxZqampEr169RHFxsTh//vw1\nT+7m5eVpenLXjIx84kprzIWCuVAwFwpXa6fqcfybNm1CWloa6urqMHHiRLz44ovIzMwEAEyaNAkA\nMGXKFNjtdrRr1w4ffvghfv/73zfZhlHm6iEiMhPec5eIyMdwkjYvcXGP39cxFwrmQsFcqMfCT0Tk\nY9jqISIyKbZ6iIjIKSz8BsP+pYK5UDAXCuZCPRZ+IiIfwx4/EZFJscdPREROYeE3GPYvFcyFgrlQ\nMBfqsfATEfkY9viJiEyKPX4iInIKC7/BsH+pYC4UzIWCuVCPhZ+IyMewx09EZFLs8RMRkVNY+A2G\n/UsFc6FgLhTMhXos/EREPoY9fiIik2KPn4iInMLCbzDsXyqYCwVzoWAu1HO58J84cQJDhw7Frbfe\nivvuuw/V1dWXrVNaWoohQ4bgtttuQ1RUFN5++21VwRIRkXou9/inT5+OLl26YPr06UhPT0dVVRUW\nLFjQZJ3jx4/j+PHjsNlsOH36NPr164f169cjMjKyaRDs8RMRtZjmPf4NGzZg/PjxAIDx48dj/fr1\nl63TvXt32Gw2AED79u0RGRmJY8eOubpLIiJyA5cLf0VFBbp16wYA6NatGyoqKq66fklJCb799lvE\nxcW5ukufwP6lgrlQMBcK5kI9/6t9c+jQoTh+/PhlX583b16TZT8/P/j5+V1xO6dPn8bYsWOxePFi\ntG/fvtl1kpOTERISAgAIDAyEzWZDfHw8AOUPzWXfWm5glHj0XC4oKDBUPHouFxQUGCoeLZcdDgey\nsrIAoLFeusLlHn9ERAQcDge6d++OH3/8EUOGDMHBgwcvW6+mpgb3338/hg0bhrS0tOaDYI+fiKjF\nNO/xjxo1CsuXLwcALF++HA888MBl6wghMHHiRPTp0+eKRZ+IiLTlcuF/4YUX8D//8z+49dZb8eWX\nX+KFF14AABw7dgwjRowAAGzduhUrV65Ebm4uYmJiEBMTA7vd7p7IvdSlbQ5fxlwomAsFc6HeVXv8\nV9OpUyd88cUXl329R48eyMnJAQDceeedqK+vdz06H3RxL9fXMRcK5kLBXKjHK3cNprkL4XwVc6Fg\nLhTMhXos/EREPoaF32BKSkr0DsEwmAsFc6FgLtQzxLTMNpsNe/fu1TsMIiJTsVqtjdc1tIQhCj8R\nEWmHrR4iIh/Dwk9E5GM0Lfx2ux0REREIDw9Henp6s+tMnToV4eHhsFqt+Pbbb7UMT1PXysWqVatg\ntVoRHR2NQYMGobCwUIcoteHM8wIAdu7cCX9/f6xbt07D6LTlTC4cDgdiYmIQFRXl1ePZr5WLyspK\nJCYmwmazISoqqnEOG28zYcIEdOvWDX379r3iOi2um0IjtbW1IjQ0VBQXF4sLFy4Iq9UqioqKmqyT\nk5Mjhg0bJoQQYvv27SIuLk6r8DTlTC62bdsmqqurhRBCbNq0yadz0bDekCFDxIgRI8Snn36qQ6Se\n50wuqqqqRJ8+fURpaakQQoiff/5Zj1A9zplczJ49W7zwwgtCCJmHTp06iZqaGj3C9agtW7aIPXv2\niKioqGa/70rd1OyIPz8/H2FhYQgJCUFAQACSkpKQnZ3dZJ2L5/iPi4tDdXX1Nad7NiNncjFgwAB0\n6NABgMxFWVmZHqF6nDO5AIAlS5Zg7Nix6Nq1qw5RasOZXKxevRoPPvggLBYLAKBLly56hOpxzuTi\npptuwsmTJwEAJ0+eROfOneHv7/JkBIY1ePBgdOzY8Yrfd6Vualb4y8vLERwc3LhssVhQXl5+zXW8\nseA5k4uLLV26FMOHD9ciNM05+7zIzs7G5MmTAeCqU4CbmTO5OHz4ME6cOIEhQ4YgNjYWK1as0DpM\nTTiTi9TUVOzfvx89evSA1WrF4sWLtQ7TEFypm5q9PTr7YhWXjC71xhd5S36n3NxcLFu2DFu3bvVg\nRPpxJhdpaWlYsGBB4xS0lz5HvIUzuaipqcGePXuwefNmnD17FgMGDMAdd9yB8PBwDSLUjjO5mD9/\nPmw2GxwOB77//nsMHToUe/fuxQ033KBBhMbS0rqpWeEPCgpCaWlp43JpaWnjx9UrrVNWVoagoCCt\nQtSMM7kAgMLCQqSmpsJut1/1o56ZOZOL3bt3IykpCYA8obdp0yYEBARg1KhRmsbqac7kIjg4GF26\ndEGbNm3Qpk0b3HXXXdi7d6/XFX5ncrFt2zbMmjULABAaGoqePXvi0KFDiI2N1TRWvblUN912BuIa\nampqRK9evURxcbE4f/78NU/u5uXlee0JTWdyceTIEREaGiry8vJ0ilIbzuTiYsnJyeLvf/+7hhFq\nx5lcHDhwQNx7772itrZWnDlzRkRFRYn9+/frFLHnOJOLadOmiTlz5gghhDh+/LgICgoSv/zyix7h\nelxxcbFTJ3edrZuaHfH7+/sjIyMDCQkJqKurw8SJExEZGYnMzEwAwKRJkzB8+HBs3LgRYWFhaNeu\nHT788EOtwtOUM7mYO3cuqqqqGvvaAQEByM/P1zNsj3AmF77CmVxEREQgMTER0dHRaNWqFVJTU9Gn\nTx+dI3c/Z3Ixc+ZMpKSkwGq1or6+HgsXLkSnTp10jtz9xo0bh6+++gqVlZUIDg7GK6+8gpqaGgCu\n101O2UBE5GN45S4RkY9h4Sci8jEs/EREPoaFn4jIx7DwExH5GBZ+IiIfw8JPRORjWPiJiHzM/wNj\nxL+4sahrEgAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x34dbffd0>"
       ]
      }
     ],
     "prompt_number": 507
    }
   ],
   "metadata": {}
  }
 ]
}