{
 "metadata": {
  "name": "",
  "signature": "sha256:3540b80f7cce90e2f8d18309f853ba11a54fde232f02b31c17444c4d7c54fba8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, C.D. Cooper, G.F. Forsyth."
     ]
    },
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Spreading out"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Welcome to the fifth, and last, notebook of Module 4 \"_Spreading out: diffusion problems,\"_ of our fabulous course **\"Practical Numerical Methods with Python.\"**\n",
      "\n",
      "In this course module, we have learned about explicit and implicit methods for parabolic equations in 1 and 2 dimensions. So far, all schemes have been first-order in time and second-order in space. _Can we do any better?_ We certainly can: this notebook presents the Crank-Nicolson scheme, which is a second-order method in both time and space! We will continue to use the heat equation to guide the discussion, as we've done throughout this module. "
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Crank-Nicolson scheme"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The [Crank Nicolson scheme](http://en.wikipedia.org/wiki/Crank\u2013Nicolson_method) is a popular second-order, implicit method used with parabolic PDEs in particular. It was developed by John Crank and [Phyllis Nicolson](http://en.wikipedia.org/wiki/Phyllis_Nicolson). The main idea is to take the average between the solutions at $t^n$ and $t^{n+1}$ in the evaluation of the spatial derivative. Why bother doing that? Because the time derivative will then be discretized with a centered scheme, giving second-order accuracy!\n",
      "\n",
      "Remember the 1D heat equation from the [first notebook](https://github.com/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/01_Heat_Equation_1D_Explicit.ipynb)? Just to refresh your memory, here it is:\n",
      "\n",
      "\\begin{equation}\n",
      "\\frac{\\partial T}{\\partial t} = \\alpha \\frac{\\partial^2T}{\\partial x^2}.\n",
      "\\end{equation}\n",
      "\n",
      "In this case, the Crank-Nicolson scheme leads to the following discretized equation:\n",
      "\n",
      "\\begin{eqnarray}\n",
      "\\frac{T^{n+1}_i - T^n_i}{\\Delta t} = & \\nonumber \\\\\n",
      "\\alpha \\cdot \\frac{1}{2} &\\left( \n",
      "\\frac{T^{n+1}_{i+1} - 2T^{n+1}_i + T^{n+1}_{i-1}}{\\Delta x^2} \n",
      "+ \\frac{T^n_{i+1}-2T^n_i + T^n_{i-1}}{\\Delta x^2}\n",
      "\\right)\n",
      "\\end{eqnarray}\n",
      "\n",
      "Notice how the both time indices $n$ and $n+1$ appear on the right-hand side. You know we'll have to rearrange this equation, right? Now look at the stencil and notice that we are using more information than before in the update."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "![stencil-cranknicolson](./figures/stencil-cranknicolson.png)\n",
      "\n",
      "#### Figure 2. Stencil of the Crank-Nicolson scheme."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Rearraning terms so that everything that we don't know is on the left side and what we do know on the right side, we get\n",
      "\n",
      "\\begin{eqnarray} \n",
      "-T^{n+1}_{i-1} + 2\\left(\\frac{\\Delta x^2}{\\alpha\\Delta t}+1\\right)T^{n+1}_i - T^{n+1}_{i+1} = & \\nonumber\\\\\n",
      "T^{n}_{i-1} + & 2\\left(\\frac{\\Delta x^2}{\\alpha\\Delta t}-1\\right)T^{n}_i + T^{n}_{i+1}\n",
      "\\end{eqnarray}\n",
      "\n",
      "Again, we are left with a linear system of equations. Check out the left side of that equation: it looks a lot like the matrix from [notebook 2](https://github.com/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/02_Heat_Equation_1D_Implicit.ipynb), doesn't it? Apart from the slight modification in the $T_i^{n+1}$ term, the left side of the equation is pretty much the same. What about the right-hand side? Sure, it looks quite different, but that is not a problem, we know all those terms!\n",
      "\n",
      "Things don't change much for boundary conditions, either. We've seen all the cases already. Say $T_0^{n+1}$ is a Dirichlet boundary. Then the equation for $i=1$ becomes\n",
      "\n",
      "\\begin{eqnarray} \n",
      " 2\\left(\\frac{\\Delta x^2}{\\alpha\\Delta t}+1\\right)T^{n+1}_1 - T^{n+1}_{2} = & \\nonumber\\\\ \n",
      " T^{n}_{0} + & 2\\left(\\frac{\\Delta x^2}{\\alpha\\Delta t}-1\\right)T^{n}_1 + T^{n}_{2} + T^{n+1}_{0}\n",
      "\\end{eqnarray}\n",
      "\n",
      "And if we have a Neumann boundary $\\left(\\left.\\frac{\\partial T}{\\partial x}\\right|_{x=L} = q\\right)$ at $T_{n_x-1}^{n+1}$? We know this stuff, right? For $i=n_x-2$ we get\n",
      "\n",
      "\\begin{eqnarray} \n",
      "-T^{n+1}_{n_x-3} + \\left(2\\frac{\\Delta x^2}{\\alpha\\Delta t}+1\\right)T^{n+1}_{n_x-2} = & \\nonumber\\\\\n",
      "T^{n}_{n_x-3} + & 2\\left(\\frac{\\Delta x^2}{\\alpha\\Delta t}-1\\right)T^{n}_{n_x-2} + T^{n}_{n_x-1} + q\\Delta x\n",
      "\\end{eqnarray}\n",
      "\n",
      "The code will look a lot like the implicit method from the [second notebook](https://github.com/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/02_Heat_Equation_1D_Implicit.ipynb). Only some terms of the matrix and right-hand-side vector will be different, which changes some of our custom functions."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "The linear system"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Just like in [notebook 2](https://github.com/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/02_Heat_Equation_1D_Implicit.ipynb), we need to solve a linear system on every time step of the form:\n",
      "\n",
      "$$[A][T^{n+1}_\\text{int}] = [b]+[b]_{b.c.}$$\n",
      "\n",
      "The coefficient matrix is very similar to the previous case, but the right-hand side changes a lot:\n",
      "\n",
      "\\begin{align}\\left[ \\begin{array}{cccccc}\n",
      " 2\\left(\\frac{1}{\\sigma}+1\\right) & -1 & 0 & \\cdots & & 0 \\\\\n",
      " -1 & 2\\left(\\frac{1}{\\sigma}+1\\right) & -1 & 0 & \\cdots & 0 \\\\\n",
      " 0 & & \\ddots& & & \\vdots \\\\\n",
      " \\vdots & & & & 2\\left(\\frac{1}{\\sigma}+1\\right)& \\\\\n",
      " 0 & \\cdots & & & -1 & \\left(2\\frac{1}{\\sigma}+1\\right) \\end{array} \\right]\n",
      " \\cdot \n",
      " \\left[ \\begin{array}{c} \n",
      " T_1^{n+1} \\\\ T_2^{n+1} \\\\ \\vdots \\\\ \\\\ T_{N-2}^{n+1} \\end{array} \\right]\n",
      " =\n",
      " \\left[ \\begin{array}{c} \n",
      " T_0^n + 2\\left(\\frac{1}{\\sigma}-1\\right)T_1^n + T_2^n \\\\  T_1^n + 2\\left(\\frac{1}{\\sigma}-1\\right)T_2^n + T_3^n \\\\ \\vdots \\\\ \\\\  T_{n_x-3}^n + 2\\left(\\frac{1}{\\sigma}-1\\right)T_{n_x-2}^n + T_{n_x-1}^n \\end{array} \\right]\n",
      " +\n",
      " \\begin{bmatrix}\n",
      "  T_0^{n+1}\\\\\n",
      "  0\\\\\\\\\n",
      "  \\vdots\\\\\\\\\n",
      "  0\\\\\n",
      "  q\\Delta x\n",
      " \\end{bmatrix}\n",
      " \\end{align} \n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's write a function that will create the coefficient matrix and right-hand-side vectors for the heat conduction problem from [notebook 2](https://github.com/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/02_Heat_Equation_1D_Implicit.ipynb): with Dirichlet boundary at $x=0$ and zero-flux boundary $(q=0)$ at $x=L$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy\n",
      "from scipy.linalg import solve"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def generateMatrix(N, sigma):\n",
      "    \"\"\" Computes the matrix for the diffusion equation with Crank-Nicolson\n",
      "        Dirichlet condition at i=0, Neumann at i=-1\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    N: int\n",
      "        Number of discretization points\n",
      "    sigma: float \n",
      "        alpha*dt/dx^2\n",
      "    \n",
      "    Returns:\n",
      "    -------\n",
      "    A: 2D numpy array of float\n",
      "        Matrix for diffusion equation\n",
      "    \"\"\"\n",
      "   \n",
      "    # Setup the diagonal\n",
      "    d = 2*numpy.diag(numpy.ones(N-2)*(1+1./sigma))\n",
      "    \n",
      "    # Consider Neumann BC\n",
      "    d[-1,-1] = 1+2./sigma\n",
      "    \n",
      "    # Setup upper diagonal\n",
      "    ud = numpy.diag(numpy.ones(N-3)*-1, 1)\n",
      "    \n",
      "    # Setup lower diagonal\n",
      "    ld = numpy.diag(numpy.ones(N-3)*-1, -1)\n",
      "    \n",
      "    A = d + ud + ld\n",
      "    \n",
      "    return A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def generateRHS(T, sigma):\n",
      "    \"\"\" Computes right-hand side of linear system for diffusion equation\n",
      "        with backward Euler\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    T: array of float\n",
      "        Temperature at current time step\n",
      "    sigma: float\n",
      "        alpha*dt/dx^2\n",
      "    \n",
      "    Returns:\n",
      "    -------\n",
      "    b: array of float\n",
      "        Right-hand side of diffusion equation with backward Euler\n",
      "    \"\"\"\n",
      "    \n",
      "    b = T[1:-1]*2*(1./sigma-1) + T[:-2] + T[2:]\n",
      "    # Consider Dirichlet BC\n",
      "    b[0] += T[0]\n",
      "    \n",
      "    return b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will solve the linear system at every time step. Let's define a function to step in time:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def CrankNicolson(T, A, nt, sigma):\n",
      "    \"\"\" Advances diffusion equation in time with Crank-Nicolson\n",
      "   \n",
      "    Parameters:\n",
      "    ----------\n",
      "    T: array of float\n",
      "        initial temperature profile\n",
      "    A: 2D array of float\n",
      "        Matrix with discretized diffusion equation\n",
      "    nt: int\n",
      "        number of time steps\n",
      "    sigma: float\n",
      "        alpha*td/dx^2\n",
      "        \n",
      "    Returns:\n",
      "    -------\n",
      "    T: array of floats\n",
      "        temperature profile after nt time steps\n",
      "    \"\"\"\n",
      "    \n",
      "    for t in range(nt):\n",
      "        Tn = T.copy()\n",
      "        b = generateRHS(Tn, sigma)\n",
      "        # Use numpy.linalg.solve\n",
      "        T_interior = solve(A,b)\n",
      "        T[1:-1] = T_interior\n",
      "        # Enforce Neumann BC (Dirichlet is enforced automatically)\n",
      "        T[-1] = T[-2]\n",
      "\n",
      "    return T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And we are good to go! First, let's setup our initial conditions, and the matrix"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L = 1.\n",
      "nx = 21.\n",
      "alpha = 1.22e-3\n",
      "\n",
      "dx = L/(nx-1)\n",
      "\n",
      "Ti = numpy.zeros(nx)\n",
      "Ti[0] = 100\n",
      "\n",
      "sigma = 0.5\n",
      "dt = sigma * dx*dx/alpha \n",
      "nt = 10\n",
      "\n",
      "A = generateMatrix(nx, sigma)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Check the matrix..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [-1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.  0.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6. -1.\n",
        "   0.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.  6.\n",
        "  -1.]\n",
        " [ 0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0. -1.\n",
        "   5.]]\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Looks okay! Now, step in time"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T = Ti.copy()\n",
      "T = CrankNicolson(T, A, nt, sigma)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And plot,"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "from matplotlib import rcParams\n",
      "rcParams['font.family'] = 'serif'\n",
      "rcParams['font.size'] = 16\n",
      "\n",
      "x = numpy.linspace(0,L,nx)\n",
      "\n",
      "plt.plot(x, T, color='#003366', ls='-', lw=3);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAELCAYAAAAybErdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FPX9x/HXJwnhFFAO5ZLDAwUFFEsVUeOF2opVgZ5W\na71atT9tf1p/YlUs1qv21HrUu9W2gnjgiUcNCBZEkRu8CqICKiD3kYR8fn/sZNlsFkiyszvZzfv5\neOSxO7Mzk88Oy74z35nvd8zdERERSVQQdQEiItLwKBxERKQGhYOIiNSgcBARkRoUDiIiUoPCQURE\nathlOJhZJzN7ycwqs1GQiIhEb6fhYGZnAm8CvYAddogws1ZmdqeZLTKz+WY20cz6pFiuiZmNMbOF\nZjbXzKaa2ZFpvwsREQnVro4cfgmcQCwgbCfLjQP6AQPcvS8wHSg1s85Jy90BjASGuPvBwIPAy2bW\nvz7Fi4hIZtjOekibWYG7V5rZw8DZ7l4jTMzsRGAicJy7lwbzmgArgH+6+6XBvN7AAuA8d384Yf15\nwBJ3PzWsNyUiIunZ6ZGDu9fmPMNwoAyYkrBeOTA1eK3KGcSOPl5PWv91YKiZtahNwSIiknlhXK3U\nD1jm7hVJ85cAe5pZ+4TltgFLk5ZbDBQBNc5RiIhINMIIh/bA+hTz1wWP7RKW2+Q127GSlxMRkYip\nn4OIiNQQRjisBHZLMb918LgqYbmWZpZ81VPyciIiErGiELYxBxhoZkVJ5x16AivcfWUwPRv4LtCN\n6ucdegLlxK5kqsbMdLMJEZE6cveddT2olbocOezoi3o80ASId2Yzs+JgenzCck8F2zg2af1jgZfd\nfVOqjX/36vtw90b9c/3110deQ0P50b7QvtC+2PlPWOoSDimTyN1fIdbPYYyZNQ9mX0PsaOCmhOXe\nB/4KXG1m7QDM7FxiRw7X7OiXvvjmPMorttWhTBERSdeuhs+4zczeBYYBbmbvmtnMoJNbohHEmpdm\nmdkC4HCgxN2XJy33M2K9qaea2VzgfGCou8/ZUQ1rN2xmyqwP6/auREQkLTs95+Duv6zNRtx9I3Bp\nLZarAK4Nfmrt2cmzOfaw3nVZJa+UlJREXUKDoX2xnfbFdtoX4dvp8BlRMzNn4IXs07UDHzw1hpoX\nOomISCIzw7N8QjoyH336JYuWrIi6DBGRRiMnwgHg2ck7PC0hIiIhy51weEPhICKSLQ0+HKrOM7w5\n5yNWrdkQcTUiIo1Dgw+HIw7uBUBlpfPC1HkRVyMi0jg0+HAYdnS/+PMJk2dHWImISOPR8MPhqO3h\nMHHaAsrKk28bISIiYWvw4dCnVyd6dondL2j9xi1Meuf9iCsSEcl/DT4czKza0YOuWhIRybwGHw4A\npx1dPRwacq9uEZF8kBPhcNQh+9G6ZTMAlixbxbyPlkVckYhIfsuJcChuUsTJg/vGp5/VVUsiIhmV\nE+EA6LyDiEgW5Uw4nDL4IAoKYr2lp89bwher10VckYhI/sqZcGjXthVH9t8HAHfn+SlzI65IRCR/\n5Uw4AJx2dP/48wkapVVEJGNyKhwSzzu8PG0BW7aWR1iNiEj+yqlw6N1jL/bbuyMAm7aU8frb70Vc\nkYhIfsqpcABdtSQikg05Hw7qLS0iEr6cC4cjB+zL7q1bAPDp518x671PIq5IRCT/5Fw4NCkq5JTB\nB8Wn1bQkIhK+nAsH0HkHEZFMy8lwOHlwX4oKY6W/veBjln25JuKKRETyS06GQ9vdWnDUIfvFp9Vb\nWkQkXDkZDlD9Hg+6t7SISLhyNhyGJQyl8epbi9i0pSzCakRE8kvOhsM+XTtwYM9OAGzZWs5rby2M\nuCIRkfyRs+EAMOyog+PPddWSiEh4cjscEpqWnntjLpWVlRFWIyKSP3I6HI44uBft2rQEYPnKtbyz\ncGnEFYmI5IdQwsHMDjOzF81sgZnNMbPpZjYiaZlWZnanmS0ys/lmNtHM+qTzewsLC/jmEDUtiYiE\nLe1wMLMewGvAF8BB7t4PeBAYa2anJiw6DugHDHD3vsB0oNTMOqfz+6v1ltYNgEREQhHGkcM3gN2A\n37t7JYC73wusA74PYGYnAicB17n7lmC9MUAhMCqdXz708D40KSoEYNb7n/DJitXpbE5ERAgnHCqC\nxyZVM8zMiH3xV21/OFAGTKlaxt3LganBa/XWulVzSgbuH59+Tr2lRUTSFkY4/BNYBPzKzFqaWQGx\no4EmwD3BMv2AZe5ekbTuEmBPM2ufTgHV7y2t3tIiIulKOxzcfT1wPNAcWAl8DpwLDHX30mCx9sD6\nFKuvCx7bpVPDsIShNP494z02bNqyk6VFRGRXwjgh3RuYASwGdnf3DsA1wJNmdnK626+N7p3acfC+\nXQAoK6/glenqLS0iko6iELYxBmgNXObuWwHc/XEz+w7wSHA10kpgrxTrtg4eV+1o46NHj44/Lykp\noaSkJOVyw47qx9wPPwNiVy2dcewhdX0fIiI5p7S0lNLS0tC3a+neg9nMFgLu7n2S5t8KXAnsD1xB\nrKmpZeJ5BzN7Fhjo7ikvZzUzr2190+b+lyPOvRWADrvvxvKXbqOwMKf7+ImI1JmZ4e6W7nbC+Pb8\nHOhsZoVJ87sDlcBqYDyxE9RHVr1oZsXB9PgQamBQ3x503GM3AL78aj1vzV8cxmZFRBqlMMLhDmLN\nQ7+ummFmxwJnAI+7+2p3fwWYCIwxs+bBYtcA5cBNIdRAQUEBpw7R7UNFRMIQxtVK44GTgSOC4TPm\nAn8kdjnrjxIWHQHMAWaZ2QLgcKDE3ZenW0OVxKuW1FtaRKT+0j7nkEl1OecAsHHzVtod/wu2lsVO\nayye8Bt6dE6rC4WISE5pSOccGoyWzZty3GEHxKd19CAiUj95FQ6QfG9phYOISH3kXTicmjBK66SZ\n77Nuw+YIqxERyU15Fw5d99ydQ3p3A6C8YhsTpy2IuCIRkdyTd+EA1a9amjBJA/GJiNRVXoZD4iit\nz74xh61l5RFWIyKSe/IyHA49YG96dI4N9Lp2w2YNxCciUkd5GQ5mxrdPGBifHvfqOxFWIyKSe/Iy\nHABGJoTDM5Nmq2lJRKQO8jYcBh7YnZ5dYr2j1bQkIlI3eRsOZsbI4w+NT6tpSUSk9vI2HKB609LT\npbPUtCQiUkt5HQ6JTUvrNm5R05KISC3ldTgkNy2NfUVNSyIitZHX4QDw7RMPiz9/ZpKalkREaiPv\nw+HQA/ZW05KISB3lfTgkd4hT05KIyK7lfThAcoc4NS2JiOxKowiHQw/Ym14JTUsvaxhvEZGdahTh\nYGbVjh7GvTozwmpERBq+RhEOoKYlEZG6aDThoKYlEZHaazThoKYlEZHaazThAOoQJyJSW40qHA7p\n3U1NSyIitdCowsHMqh09qEOciEhqjSocoPpVSxMm6w5xIiKpNLpwUNOSiMiuNbpwUNOSiMiuNbpw\ngJpNS1u2qmlJRCRRowyHQ3p3Y5+uHQA1LYmIpBJaOJjZcDObZGZvm9lHZjbDzM5KeL2Vmd1pZovM\nbL6ZTTSzPmH9/jrWmtQhTk1LIiKJQgkHM/s5MAr4nrsfBvQG3geOS1hsHNAPGODufYHpQKmZdQ6j\nhrpKvMfDM2paEhGpJu1wMLMewM3Ahe6+DMDdK4ArgDuDZU4ETgKuc/ctwapjgEJioZJ1AxKaltar\naUlEpJowjhx+CHzl7tXaZtx9ubtXDWA0HCgDpiS8Xg5MDV7LOjUtiYjsWBjhMBj4ODjnMNnMFprZ\nVDM7N2GZfsCy4Igi0RJgTzNrH0IddaamJRGR1MIIh25AX+AXwAh3PxD4A/BXM6tqMmoPrE+x7rrg\nsV0IddSZmpZERFILIxyaAS2BK939CwB3fwJ4BhhlZs1D+B0ZEesQt/3oYeyrb0dYjYhIwxFGOKwH\nHJiVNH8W0ALoA6wEWqdYt2reqhDqqJeRxyd2iJujpiUREaAohG0sAvpTM2i2BY8GzAYGmllR0nmH\nnsAKd1+5o42PHj06/rykpISSkpIQSt5uQO9u7NutIx9+8kW8aem0Y/qH+jtERDKltLSU0tLS0Ldr\n7p7eBsy+BzwGHOXuUxPmPwZ8C+gADAEmAse6+6Tg9WJgBfCYu/9sB9v2dOurjVF/eYqbH3oJgB+c\nMohHx5yX8d8pIpIJZoa7W7rbCaNZ6XFgBnCjmbUEMLOjiF2i+ht33+zurxALhzEJ5yCuAcqBm0Ko\nIS1qWhIRqS7tcHD3SuBk4ANgvpktItb57RJ3vzlh0RHAHGCWmS0ADgdK3H15ujWkq6ppCWJXLU38\nz/yIKxIRiVYow2e4+1fufqG793D3A9y9v7s/kLTMRne/1N17u3sfdz/J3ReG8fvTFesQd2h8etxr\n6hAnIo1boxyVNZVvn7D9Hg9qWhKRxk7hEOi/f1c1LYmIBBQOATOrNpyGmpZEpDFTOCSofoc4NS2J\nSOOlcEigpiURkRiFQ4LkpqWxGsZbRBophUOS6k1Ls9m8pSzCakREoqFwSNJ//67st3esaWnDpq1M\n1DDeItIIKRySmFm14TR0hzgRaYwUDikk3uNBTUsi0hgpHFLot1/1pqWXdNWSiDQyCocUkq9aenDC\n1J0sLSKSfxQOO/CjYYPjz5+fMo+Pl0d2szoRkaxTOOzAvt06MvTwPgC4O3998o2IKxIRyR6Fw078\ndMQx8ecPTJhKWXnFTpYWEckfCoedOHXIwXTp2BaAz1et4+nSWRFXJCKSHQqHnSgqKuSC04fEp+9+\nYlKE1YiIZI/CYRfOP30IhYWx3VT6zvssWrIi4opERDJP4bALXTruzmlH94tP3zNeRw8ikv8UDrXw\nk+HbT0w/8tw0NqnHtIjkOYVDLZww6AD26doBgDXrN/H4yzMirkhEJLMUDrVQUFDAT4YfHZ+++4nJ\nEVYjIpJ5Coda+tGwwTQtLgJgxoIlvLPw44grEhHJHIVDLbVv26raUN73jNfRg4jkL4VDHSQ2Lf3j\npbdYu2FzhNWIiGSOwqEOBvffh4P37QLApi1l/P35aRFXJCKSGQqHOjCzauMt3T1+Eu4eYUUiIpmh\ncKijH5w8iJbNmwKw4L/LmTLrw4grEhEJn8Khjlq3as5Zp3w9Pq3xlkQkHykc6iHxxPQTr83ki9Xr\nIqxGRCR8Cod6GNC7G4cf3BOA8optPDThzYgrEhEJV0bCwczeMLNKM9s7E9tvCH6aMN7SvU+9QWVl\nZYTViIiEK/RwMLPhwJFAjct4zKyVmd1pZovMbL6ZTTSzPmHXkA0jTxjI7q1bALD4s5W8PG1BxBWJ\niIQn1HAws2LgFuAFwFIsMg7oBwxw977AdKDUzDqHWUc2NG9WzLnDBsenNd6SiOSTsI8cLiH2hV9j\n2FIzOxE4CbjO3bcEs8cAhcCokOvIiovO3H5i+rkpc1i6YnWE1YiIhCe0cDCzPYArgKtJfdQwHCgD\nplTNcPdyYGrwWs7Zv/ueHD/oAAAqK537n56yizVERHJDmEcO1wF/d/dPdvB6P2CZu1ckzV8C7Glm\n7UOsJWsST0zf//QUyiu2RViNiEg4QgkHM9sPGAn8ZieLtQfWp5hf1UmgXRi1ZNtpx/SnU/s2ACxf\nuZZnSmdFXJGISPrCOnK4FbjZ3VN9+ee1JkWFnH/6kPj0PU/qxLSI5L60w8HMjgL6Avekejnh+Upg\ntxTLtA4eV6VbS1QuOH0IBQWxt/raW4t4/+PPI65IRCQ9RSFs4wRiVxzNMItnwV7B4wtmVkbsaqTZ\nwEAzK0o679ATWOHuK1NtfPTo0fHnJSUllJSUhFByuLrttQenDunHhMmzAbj3ycn87ucjI65KRBqD\n0tJSSktLQ9+uZWLIaTO7Hrge6OHuS4N5JwITgWPdfVIwrxhYATzm7j9LsR3PlSGxX3pzHqf8zx0A\n7N66BZ+9cCvNmxVHXJWINDZmhrunumK0TjI1tpIlPeLurxALhzFm1jyYfQ1QDtyUoTqyZujhfejZ\nJXbB1VfrNjH21XcirkhEpP7C7iF9ipm9C1xEbPiMF8xsZsIiI4A5wCwzWwAcDpS4+/Iw64hCQUEB\nF51xVHz6nvEayltEcldGmpXCkkvNSgBfrF5H12/8X7yvw7uP/YoBvbtFXJWINCYNvVmpUeq4R2tG\nHH9ofPqe8bqsVURyk8IhZIn3mH70xems27A5wmpEROpH4RCyIQP2pW+v2CCzGzdv5dEXp0dckYhI\n3SkcQmZm1W4jes/4yeTSeRMREVA4ZMQPv3k4LYI+DnM//Iw3Z38UcUUiInWjcMiANq2a8/2TB8Wn\n79ZlrSKSYxQOGZLYtDTu1Zl8vmrdTpYWEWlYFA4ZMvDA7nytTw8AysoruOG+56ItSESkDhQOGXTd\nBd+MP//rU2+waMmKCKsREak9hUMGfXPIwRx7WG8Atm2r5Ko/PxlxRSIitaNwyCAz47eXbb899oTJ\ns5k88/0IKxIRqR2FQ4YNPLA7Z53y9fj0FX8cT2VlZYQViYjsmsIhC268+Fs0LY7dV2nGgiU8/vLb\nEVckIrJzCocs6N6pHZd997j49NV/eZotW8sjrEhEZOcUDlly9bmn0K5NSwA+Xr6KO8e+HnFFIiI7\npnDIkra7teC6C06NT//mwRdZvXZjhBWJiOyYwiGLfjL8aPbt1hGANes3ceMDz0dckYhIagqHLCpu\nUsQtl54Rn75zbCkfffplhBWJiKSmcMiyM487hMH99gGgvGIbo/7yVMQViYjUpHDIMjPj9su3d4wb\n+8o7TJv73wgrEhGpSeEQgSP67VPtXtNX/PEJ3RBIRBoUhUNEbr70DJoUFQIwdfZHPF06K+KKRES2\nUzhEZN9uHbl45DHx6avueJLyim0RViQisp3CIULXnvdN2rRqDsAHS7/g3vGTI65IRCRG4RChdm1b\ncc2PT4lP33Dfc6zdsDnCikREYhQOEfvZd46je6d2AKxcs4FbHn4p4opERBQOkWvWtAk3XXJ6fPqP\n/3yNpStWR1iRiIjCoUH47tDDGHjg3gBs2VrOr+56JuKKRKSxUzg0AAUFBdx+2Yj49KMvTufdRUsj\nrEhEGjuFQwNRclhvhh3VDwB354o/qWOciERH4dCA3Po/Z1JYGPsn+feM93hx6ryIKxKRxkrh0IAc\n2LMTF5w+JD79yz8/SYU6xolIBNIOBzMbYGb3mdkCM5tjZvPN7E9m1j5puVZmdqeZLQqWmWhmfdL9\n/flm9IXDaNWiKQDz/7uMh5/7T8QViUhjFMaRw7+AtsBAd+8HnAgMBaaaWbOE5cYB/YAB7t4XmA6U\nmlnnEGrIG3u2a81VZ58Un7727mfYsGlLhBWJSGMURjhUAle5+2YAd18G/BbYD/gGgJmdCJwEXOfu\nVd90Y4BCYFQINeSVX5x1Ip07tAVgxap1/O7RVyKuSEQamzDCoZ+7J9+QYHnw2DZ4HA6UAVOqFnD3\ncmBq8JokaNGsmBt/elp8+ra/vczylWsjrEhEGpu0w8HdK1LM3h9woGokuX7AshTLLgH2TD4/IXD2\nN4/g4H27ALBpS5k6xolIVoV+tZKZFQLnAfe7+4fB7PbA+hSLrwse24VdR64rLCzg9su3d4x7cMJU\nHtHJaRHJkkxcynotsBW4PAPbblSGHt6H4cdtv2Pchb95lP/M+SjCikSksSgKc2Nmdi4wAiipOkEd\nWAnslWKV1sHjqh1tc/To0fHnJSUllJSUpF1nLnno+nN47+MVzPtoGWXlFZxx5T28/bdRdN1z96hL\nE5EGoLS0lNLS0tC3a2EN0WBmPwR+CRzv7l8kvXYPcC7QMvG8g5k9S+wS2JSXs5qZawgJWPzZSr52\n9k2sWrsRgEMP2Js37r+SFs2KI65MRBoaM8PdLd3thNKsZGZnkRQMZnaqmV0QLDIeaAIcmbBOcTA9\nPowa8lnPLu154taLKAqG1pi5aCk/vuERjb0kIhkTRg/pHwD3AQ8DQ83srCAshgGdAdz9FWAiMMbM\nmgerXgOUAzelW0NjUHJYb+648rvx6cdfeZubHnwxwopEJJ+l3axkZquI9WdIPoxx4AZ3/3WwXEvg\nVmI9qLcBnwCXu/vCnWxbzUpJLrn1H9w1blJ8+qnbf8rpJQMirEhEGpKwmpVCO+eQCQqHmsortnHS\npX/i9bffA6Bl86b856Gr4n0iRKRxa1DnHCR7mhQVMu6WC+nVJdZvcOPmrZz2i7/w5VepupGIiNSP\nwiEHtWvbigm/vyQ+euuSZasYcdW9lJWn6qwuIlJ3Cocc1XefzvzjxvMwix09Tp75AT+77V+6gklE\nQqFwyGHDju7PTZecHp/+61NvcNe40ugKEpG8oXDIcVedcxLfP3lQfPqy343l3zMWRViRiOQDXa2U\nBzZvKeOYC3/HjAVLANijTUveeuRq9unaIdrCRCTrdLWSxDVvVsxTt/+ETu3bALB67UaG/fwvrNuw\neRdrioikpnDIE1067s7Tt/+UpsWxsRQXLl7O93/1ANu2VUZcmYjkIoVDHhl0UE8euPbs+PTzU+Zy\nzV1PR1iRiOQqhUOe+cEpX+eqc06KT9/6yEQefWFahBWJSC7SCek8tG1bJd/637t4fspcAJoWFzH5\nr1cw6KCeEVcmIpmmE9KyQ4WFBfzjxvPo06sTAFvLKjj9irv5ZMXqiCsTkVyhcMhTrVs1Z8LvL2H3\n1i0AWL5yLQN/eBMvT1sQcWUikgsUDnlsn64dqt0k6Muv1nPSpX9i1F+eoqJiW8TViUhDpnMOjcDr\nb7/H96+5nxWr1sXnHdl/H/75m/PpttceEVYmImHT/RykTr5YvY4fXvdQtWalPdq05OHrz2HY0f0j\nrExEwqRwkDqrrKzktr+9zK/ufqZa57hf/OAEbr70DIqbFEVYnYiEQeEg9TZl1od875r7+fTzr+Lz\nvtanB4/ffAE9g5sIiUhuUjhIWlat2cCPbniY596YG5/XplVzHrj2bIYff2iElYlIOhQOkjZ35w+P\nvcpVdzxJRUIz0yUjS7j98hE0a9okwupEpD4UDhKat+Yt5juj7mPJslXxeQP278bYWy5gv733jLAy\nEakrhYOEas36TZz367/x5Ovvxue1atGUe0edVe1mQiLSsCkcJHTuzt1PTOLnvx9HWXlFfP553zqS\nP1/5XVo0K46wOhGpDYWDZMy7i5by7avv48NPvojP69urM2NvuYA+vTpHWJmI7IrCQTJq/cYtXHTT\no/xz4oz4vOImRYw84VAuHlHCEf16YZb2509EQqZwkIxzdx58ZiqX/vZfbNlaXu21Aft34+KRx/D9\nkwfRsnnTiCoUkWQKB8maeR9+xkU3Pcabcz6q8VqbVs0559QjuHjEMfTusVcE1YlIIoWDZN3MRUu5\na1wp/3jpLTYnHUkAHD/oAC4ecQynHd2foqLCCCoUEYWDROardRt55Ln/cNcTk/hg6Rc1Xu/SsS0X\nnXk0558+hE7t20RQoUjjpXCQyFVWVvLaW4u464lJTJg8m8rK6v9WRYUFnHncIVwysoSjDtlPJ7BF\nskDhIA3KJytWc++Tk7nv6Sl8sXp9jdf79urMxSOP4TsnHka7tq0iqFCkccjJcDCzjsAfgIHBrLnA\n5e7+2Q6WVzjkmLLyCsa/NpO7npjElFkfplxm324dGdS3B18/qCdfP6gnA/bvStNijeMkEoacCwcz\nKwZmAIuA7wWzHwQGA4e4+8YU6ygccticDz7l7icm8fcXprNx89YdLtekqJAB+3erFhj7dutAQYHu\nYitSV7kYDhcA9wK93H1JMG9P4DPg/9z99hTrKBzywNoNm/n789N49MXpzFy0lPJa3L+67W4t4mFR\n9dhh992yUK1IbsvFcHgJ6O3uPZPmzwE2uvsRKdZROOSZrWXlzHr/U96at5jp8xfz1vwlKa94SqVn\nl/Z8rU93enRqR+cObWM/7dvQuUNbOrVvoyHGRcjNcFgGLHL345LmTwCOc/caZykVDjGlpaWUlJRE\nXUbGrF67kbeCoJg+bzHT5y1m1doarYwx65fBbqnHd9q9dQs6t29L5w5tgvBoU226U/s2dGrfJm9u\nh5rvn4u60L7YLqxwyOb/kvZAzctYYB3QwsyauvuOG6YbsXz/4O/RpiUnDz6IkwcfBMSG7Vj82Uqm\nzwsCY/5iZi5aytayip2Gw1frNvHVuk3M/++ynf6+tru1oFWLprRq3pSWzYtp1bwZLZsX07J5wrwW\nzWjZrJhWLZomzN/+WrPiIooKC2lSFPspKizY/ryo+vxMXcKb75+LutC+CF82w0GHAFIrZkavrh3o\n1bUD3wvuJVFWXsHcDz/j17++gUEnfItlK9ew7Mu1LPtyDctWrmX5yrVsS7ib3c6sWb+JNes3ZfIt\nVFMYBEe1ACkspLDAKCgooKDAKDDDLPZYUGDxeQUFBZixfb4VxF//5N2pvPrJbfHwMWP7c1LMS3q9\nKrOSwys5y2q8zs6Xr490A/T9aW/xzto70y9E4rIZDiuBVGcUWxM756CjBtmh4iZFDDywO4f03ptr\nzvtGjdcrKytZuWZDtcBY9mVigKxh+cp1rFi1tkZnvUzbtq2y1sFVJ5+v5tPZNce7apSWfc4HFXN3\nvZzUWjbPObwIHJDihPRcYL27D06xjo42RETqKNfOOTwJ3Gtm3d39Y4hfynoA8H+pVgjjDYqISN1l\n88ihCfA2sBD4AbFzEA+wvRNc9hqBRURkp7LWBdXdy4ETgW3AguCnFbHLWBUMIiINSIMeeE+krszs\nRmAUcK67PxJ1PSJhMbNOwEPAUHfP+B/2WR+8xsw6mtljZrYo+BlnZl1quW4TMxtjZgvNbK6ZTTWz\nIzNdc6bUd1+YWSczu8HMZgf7YaGZjTezg7JRdyak87lI2EZX4BfEmixz9q+edPeFmfU3s2fM7J3g\ns7HIzG7NZM2Zkub3RSczuz/YB7PNbJ6ZXW1mOdcL0szOBN4EelHHz3a9vzfdPWs/QDEwG3icWDAV\nAA8D7wMta7H+PcQG7msXTJ8HbAT6Z/N9RL0vgv3wHtAlmG4KjA32xUFRv7dsfy4StvM34FmgEjg7\n6vcVxb4gdg7vM+CIhHkXA/+N+r1lc18Ey74LzAF2D+YNADYBv436vdVjX0wD9gnef2Ud163X92a2\n3+AFwX/ycnGnAAAEkUlEQVTcHgnz9gQqgCt2sW5vYucrfpQ0fx7wXNT/eFneF3cDP06a1yvY3p+j\nfm/Z3BcJyw8EPgSG5ng4pPO5MGIXfPxv0vwi4KSo31uW90WfYN3LkuY/DSyL+r3VY18UBI91Cod0\nvjez3aw0HPjYg1FZAdz9c2Inp4fvYt0ziH34X0+a/zow1MxahFhnNqSzLy4l1vaYaHnw2DasArMo\nnX1R5XfEzjWUhV5ddqWzL4YQ+zJ4LnGmu1e4+8SQ68yGdPZFRfCYPBpjEyJoTk+Xu9e3F2W9vzez\nvZP6AYtTzF8CHFyLdbcBS5PmLyb2l1GfdIvLsnrvC3ff5kH8J9g/eCxNu7LsS+dzgZmdDjR197Eh\n1xWFdPZFVUfStsE5h3lBW/sYM2sWZpFZks7/kfeBfwAXmVl3ADM7DjgBuCPcMhu0en9vZvvETDqD\n77UHNqX4UlwXPLYLqcZsCXsgwguJHSr+PYzisqze+yLoP3ML8OMM1pdN6XwuugWP/wS+7e5vm1k/\n4HlgEHBS6NVmVrr/R84hdufJD8zsS2KXzl/m7veEX2qDVe/vzWwfOeTsFSQZENq+MLPjgW8T+0Io\nD2u7WZTOvvgpMM/d3wyrmIilsy+qjg7ud/e3Adx9DnArcKKZHZ1ucVlW730RHCmVAocB3d29C1AC\njDKzUaFUl+eyHQ7pDL63EmhpNYdvbB08rgqhvmwKZSBCM+tP7CTVMHdfFF55WVWvfWFmbYkNvZJq\n+JVcHXolnc9F1V/Zs5LmV00flmZt2ZbOvvgxcCRwpbsvB3D3d4HbgTHB/5vGoN7fm9kOhzlAzxTz\newK7GlJxNrF6uyXN7wmUEztJlUvS2RcABE0GTwHfcfdpIdaWbfXdF4cTO/E4zszeNbN3gfuC134d\nzPtVuKVmXDqfi4XBY/L/6207mN/QpbMvqs5JfJA0/wNifzjkWlDWV72/N7P9YXkS6F51ggiqDb43\nPnFBM9szKe2eInaYeWzSNo8FXvbcG4IjnX1RFQxPA2dVNakEnX5ysT21XvvC3V9y973d/ZCqH+D8\nYNFrg3k3Zuk9hCWdz8ULxIIg+a/iqs6RM8IvN6PS2RefB4/dqa5qOtdaGhLtsLkt1O/NLF+r24RY\nkv0LKCQWTg8R69DVImG5I4l9yO9KWv9uqnfmOJdYZ45+2XwfUe8LYn8VfRnsj7MSfi4HXo/6vWX7\nc5G0rRJi17efE/X7imJfELukdxmwbzDdhdhfyy9F/d6yuS+AHsBaYCLQKpi3N7G+MB8Qu7ot8vdY\nj33yMDvo5xD292ZWjxy89oPvrQe+IvYhT/QzYBww1WL3gTif2DgjczJde9jS3BejgT2Ai4j1Cq76\n+T05eNI/hM9F1TALVc1KzvZmpYGZrj9MIeyLK4ldqvmCmS0EJhH7K/tbGS49dOnsC4/1jRgErAZm\nmNls4CXgReBIz7Gbi5nZbcHnexjgwWd7ZnC1XpVQvzc18J6IiNSQayeoREQkCxQOIiJSg8JBRERq\nUDiIiEgNCgcREalB4SAiIjUoHEREpAaFg4iI1KBwEBGRGhQOIiJSw/8DJnyMyKjlUacAAAAASUVO\nRK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f207f01ec90>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Works nicely. But wait! This method has elements of explicit and implicit discretizations. Is it *conditionally stable* like forward Euler, or *unconditionally stable* like backward Euler? Try out different values of `sigma`. You'll see Crank-Nicolson is an *unconditionally stable scheme* for the diffusion equation!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "Accuracy & convergence"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Using some techniques you might have learned in your PDE class, such as separation of variables, you can get a closed expression for the rod problem. It looks like this:\n",
      "\n",
      "\\begin{eqnarray}\n",
      "T(x,t) = & \\nonumber \\\\\n",
      "100 - \\sum_{n=1}^{\\infty} & \\frac{400}{(2n-1)\\pi}\\sin\\left(\\frac{(2n-1)\\pi}{2L}x\\right) \\exp\\left[-\\alpha\\left(\\frac{(2n-1)\\pi}{2L}\\right)^2t\\right]\n",
      "\\end{eqnarray}\n",
      "\n",
      "Unfortunately, the analytical solution is a bit messy, but at least it gives a good approximation if we evaluate it for large $n$. Let's define a function that will calculate this for us:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from math import pi\n",
      "def T_analytical(x, t, n_max, alpha, L):\n",
      "    \"\"\"Computes the exact solution for 1D diffusion with T=100 at x=0 and dT/dx=0 at x=L\n",
      "        \n",
      "    Paramters:\n",
      "    ---------\n",
      "    x    : array of float\n",
      "        Spatial position\n",
      "    t    : float\n",
      "        Evaluation time\n",
      "    n_max: int \n",
      "        Number of terms to evaluate expression\n",
      "    alpha: float\n",
      "        diffusion coefficient\n",
      "    L    : float\n",
      "        Size of rod\n",
      "        \n",
      "    Returns:\n",
      "    -------\n",
      "    T    : array of float\n",
      "        Temperature at each location x\n",
      "    \"\"\"\n",
      "    T = 100\n",
      "    for n in range(1,n_max+1):\n",
      "        k = (2*n-1)*pi/(2*L)\n",
      "        \n",
      "        summation = 400/((2*n-1)*pi) * numpy.sin(k*x) * numpy.exp(-alpha*k*k*t)\n",
      "        T -= summation\n",
      "\n",
      "    return T    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And let's see how that expression looks for the time where we left the numerical solution"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T_exact = T_analytical(x, dt*nt, 100, alpha, L)\n",
      "plt.plot(x, T_exact, color='#003366', ls='-', lw=3);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAELCAYAAAAybErdAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmYFOW5/vHvMws7iIKCCAIqIqAjBvAQXDIu4BKNC2ii\nWVyiURNzsmnikiiGaKLZficaxS1qckyMBFBUBDXHYVUBlR0Eo4jKJoiyw8A8vz+qp+npaWBmurpr\nuvv+XNdc3VVdVfN00fQ99b5Vb5m7IyIikqgo6gJERKTxUTiIiEgtCgcREalF4SAiIrUoHEREpBaF\ng4iI1LLPcDCzg81sgplVZaMgERGJ3l7DwcwuBKYDhwF7vCDCzFqZ2X1mttjMFpjZRDPrnWK5UjMb\nYWaLzGyemU0zsxPSfhciIhKqfR05/BQ4nSAgbC/LjQLKgL7u3gd4A6gws05Jy90LXASc6O7HAH8B\nXjKzYxtSvIiIZIbt7QppMyty9yozexz4lrvXChMzGwxMBE5194rYvFJgFfAPd78+Nq8nsBD4trs/\nnrD+fGCZu58T1psSEZH07PXIwd3r0s8wFNgBTE1YrxKYFnut2gUERx+vJq3/KjDEzFrUpWAREcm8\nMM5WKgNWuPvOpPnLgA5m1j5huV3A8qTl3gdKgFp9FCIiEo0wwqE9sDHF/A2xx3YJy23x2u1YycuJ\niEjEdJ2DiIjUEkY4rAVap5jfJva4LmG5lmaWfNZT8nIiIhKxkhC2MRfoZ2YlSf0O3YFV7r42Nj0H\n+BrQhZr9Dt2BSoIzmWowM91sQkSkntx9b5ce1El9jhz29EU9GigF4hezmVmT2PTohOXGxrZxStL6\npwAvufuWVBv/2s0P4+4F/XP77bdHXkNj+dG+0L7Qvtj7T1jqEw4pk8jdXya4zmGEmTWPzb6V4Gjg\nroTllgAPATebWTsAM7uC4Mjh1j390henz6dy5656lCkiIuna1/AZ95jZ28C5gJvZ22b2Vuwit0TD\nCJqXZpvZQmAgUO7uK5OW+z7B1dTTzGwecBUwxN3n7qmGzzdtZcrbS+v3rkREJC177XNw95/WZSPu\nvhm4vg7L7QR+Efups+cmz+XUAUfVZ5W8Ul5eHnUJjYb2xW7aF7tpX4Rvr8NnRM3MnH7f4fDOB7J0\n7Ahqn+gkIiKJzAzPcod0ZP7z0Scsej+5hUpERDIlJ8IBgqYlERHJjpwJh3EKBxGRrGn04VBUFDSd\nvTbvPT5Zn2oIJxERCVujD4dBZYcD4O6MnzY/4mpERApDow+Hc08qiz8fN2lOhJWIiBSOxh8OJ+8O\nh5feWMj2HZURViMiUhgafTgc1a0jR3Q5CIBNW7ZT8eaSiCsSEcl/jT4czIxzTzomPj1uspqWREQy\nrdGHA8C5Jx8bf/7c5LmhjjwoIiK15UQ4nNj3CNq2bgHAh6vXM2fJRxFXJCKS33IiHEpLijlrUJ/4\n9HNTdEGciEgm5UQ4AHwlqWlJREQyJ2fC4cxBfSgpDsqduXAZKz75LOKKRETyV86EQ9vWLTjpuB7x\n6RemzouwGhGR/JYz4QDwlYQL4tS0JCKSOTkVDomntL48YxFbtu2IsBoRkfyVU+FweOcD6dX9YAC2\nba/k3zMWRVyRiEh+yqlwgJpNS7rHg4hIZuRcOCSO0vr8lLlUVVVFWI2ISH7KuXAYeMxhtG/bCoBV\n6zbw5qLlEVckIpJ/ci4ciouL+PKJGohPRCSTci4coGbTkobSEBEJX06Gw5CBvWlSWgLAnCUfsXzV\npxFXJCKSX3IyHFq3bMYp/Y+MTz+npiURkVDlZDiAmpZERDIpd8Mh4WrpV2ctYePmbRFWIyKSX3I2\nHA7teADHHtkZgB2VO3np9YURVyQikj9yNhwg6R4PaloSEQlNTodDYr/DC1PnsWuXrpYWEQlDTodD\nv16H0rFdGwDWfraJ1+e9F3FFIiL5IZRwMLP+ZvaimS00s7lm9oaZDUtappWZ3Wdmi81sgZlNNLPe\n6fzeoqIizj1ZZy2JiIQt7XAws27Av4E1wNHuXgb8BXjazM5JWHQUUAb0dfc+wBtAhZl1Suf3JzYt\naZRWEZFwhHHkcDbQGviDu1cBuPuDwAbgUgAzGwycAdzm7tXnnI4AioFb0vnlpx3fi2ZNSwFY9P5K\n3v1wTTqbExERwgmHnbHH0uoZZmYEX/zV2x8K7ACmVi/j7pXAtNhrDdaiWRMGH98rPq3bh4qIpC+M\ncPgHsBj4uZm1NLMigqOBUmBkbJkyYIW770xadxnQwczap1OA+h1ERMKVdji4+0bgNKA5sBZYDVwB\nDHH3ithi7YGNKVbfEHtsl04N5yT0O0x+eynrN2xOZ3MiIgUvjA7pnsBM4H1gf3c/ELgVGGNmZ6a7\n/bo4uP1+DOjdDYBdu6qYMH1BNn6tiEjeKglhGyOANsAP3H07gLv/08y+CjwROxtpLdAxxbptYo/r\n9rTx4cOHx5+Xl5dTXl6ecrlzTy5j5sJlQNC0dMmZx9fzbYiI5J6KigoqKipC3665e3obMFsEuLv3\nTpp/N3AjcCRwA0FTU8vEfgczew7o5+4pT2c1M69rfXOWfEjfS38FQNvWLVjz8u8oLSluwDsSEcld\nZoa7W7rbCaNDejXQycySv4m7AlXAp8Bogg7qE6pfNLMmsenRIdRAWY/OdOmwPwCfbdzC1NnvhrFZ\nEZGCFEY43EvQPPTL6hlmdgpwAfBPd//U3V8GJgIjzKx5bLFbgUrgrhBqwMxqnrWkGwCJiDRYGGcr\njQbOBL4YGz5jHvD/CE5nvTxh0WHAXGC2mS0EBgLl7r4y3RqqJY7SOm7yXNJtMhMRKVRp9zlkUn36\nHAC276ik/ek/YdOW7QAsHDWcXt0PzlR5IiKNTmPqc2g0mjYp5YyBfeLTulpaRKRh8iocoObV0uPU\n7yAi0iB5Fw5nn3A0wdBO8Nq891j72aaIKxIRyT15Fw4H7t+aQWWHAVBV5YyfOi/iikREck/ehQPU\nvMfDs5PUtCQiUl95GQ7nlfeNP39x+nw2b90eYTUiIrknL8PhqG4d6XNYMCLH1u2VvKCmJRGResnL\ncAC4eHC/+POnX34zwkpERHJP3obDRafvDofx0+axacu2vSwtIiKJ8jYcenU/mKMPV9OSiEhD5G04\nQM2jh1GvvBVhJSIiuaVgwkFNSyIidZfX4aCmJRGRhsnrcAA1LYmINERBhcMLaloSEamTvA+HxKal\nbWpaEhGpk7wPB4CLB/ePP1fTkojIvhVEOKhpSUSkfgoiHI7q1lFNSyIi9VAQ4QA1m5Y01pKIyN4V\nTDjUuCBu+nw1LYmI7EXBhMNR3TpyzBGHAEHT0vNT1LQkIrInBRMOkHxBnJqWRET2pGDDQU1LIiJ7\nVlDhoKYlEZG6KahwADUtiYjURUGHg5qWRERSK7hwUNOSiMi+FVw4AFyccPTw9CuzIqxERKRxKshw\nSGxaenH6AjUtiYgkKchw6NmtI2U9OgNqWhIRSSW0cDCzoWY2ycxmmdl/zGymmX0j4fVWZnafmS02\nswVmNtHMeof1++vrotO+EH+upiURkZpCCQcz+xFwC3CJu/cHegJLgFMTFhsFlAF93b0P8AZQYWad\nwqihvpKbljZuVtOSiEi1tMPBzLoBvwa+4+4rANx9J3ADcF9smcHAGcBt7l79LTwCKCYIlayr3bQ0\nN4oyREQapTCOHL4JrHf3GleUuftKd6++7dpQYAcwNeH1SmBa7LVIJDYtjfq3LogTEakWRjgMAj6I\n9TlMNrNFZjbNzK5IWKYMWBE7oki0DOhgZu1DqKPe1LQkIpJaGOHQBegD/BgY5u69gD8CD5lZdZNR\ne2BjinU3xB7bhVBHvalpSUQktTDCoRnQErjR3dcAuPu/gGeBW8yseQi/I2MSL4hT05KISCCMcNgI\nODA7af5soAXQG1gLtEmxbvW8dSHU0SA1xlqaNl9NSyIiQEkI21gMHEvtoNkVezRgDtDPzEqS+h26\nA6vcfe2eNj58+PD48/LycsrLy0Moebcju3agrEdn5i79iO07dvL8lLlccubxof4OEZFMqaiooKKi\nIvTtmruntwGzS4AngZPcfVrC/CeB84ADgROBicAp7j4p9noTYBXwpLt/fw/b9nTrq4s7Hx3Pzx94\nFoDzy/sy9nfXZfx3iohkgpnh7pbudsJoVvonMBP4lZm1BDCzkwhOUb3T3be6+8sE4TAioQ/iVqAS\nuCuEGtJS86wlNS2JiKQdDu5eBZwJLAUWmNligovfvufuv05YdBgwF5htZguBgUC5u69Mt4Z0Hdm1\nA8ceGZy1VN20JCJSyMLoc8Dd1wPf2ccym4Hrw/h9mXDRaf2Ys+QjAJ5+5U31O4hIQSvIUVlTUdOS\niMhuCoeY5Kal59S0JCIFTOGQ4KLTEi6Ie0UXxIlI4VI4JEhuWtqwaWuE1YiIREfhkKDWWUtTdYc4\nESlMCockNcZaUtOSiBQohUMSNS2JiCgcaulxqJqWREQUDikkNi09/fKsCCsREYmGwiGFxKalCa8t\nUNOSiBQchUMKPQ7tQN8juwBB09K4yXMirkhEJLsUDntw8eDdRw+PPDNtL0uKiOQfhcMeXHbOFyku\nDnbPpLeWsOj9yAePFRHJGoXDHnQ6sC3nnXxsfPrBMZMjrEZEJLsUDntx3bAvxZ8/8fzrbNm2I8Jq\nRESyR+GwF6cO6MkRXQ4C4LONW/jnSzMjrkhEJDsUDntRVFTENReeFJ8eOVpNSyJSGBQO+3D5uYNo\nUhrcMG/GgmW8tXh5xBWJiGSewmEf2rdtxUWnfyE+/aCOHkSkACgc6uC6obs7pp+cMENXTItI3lM4\n1MGgYw/n6MM7AbB563b+98U3Iq5IRCSzFA51YGZcO/Tk+PTI0ZNx9wgrEhHJLIVDHX3j7IG0aNYE\ngHnvfsxrc9+LuCIRkcxRONTRfq2ac+mZx8endVqriOQzhUM9JDYtPf3KLNZ9tinCakREMkfhUA/9\nenVlQO9uQDCU9+PPvxZtQSIiGaJwqKfEo4cHx0ymqqoqwmpERDJD4VBPXx3Sn/1aNQdg6fI1vDrr\nnYgrEhEJn8Khnlo2b8q3vjwwPq2OaRHJRwqHBrjmwt1NS89UzGbl2s8jrEZEJHwKhwboc3gnTjru\nCAB27qri0WemRlyRiEi4MhIOZjbFzKrM7NBMbL8xSBxv6aGxU9i1Sx3TIpI/Qg8HMxsKnADUGl/C\nzFqZ2X1mttjMFpjZRDPrHXYN2XDhqcfRvm0rAD5cvZ4Xp8+PuCIRkfCEGg5m1gT4DTAesBSLjALK\ngL7u3gd4A6gws05h1pENTZuUcuVXBsWn1TEtIvkk7COH7xF84de6n6aZDQbOAG5z922x2SOAYuCW\nkOvIiu8kdEyPnzafD1aui7AaEZHwhBYOZnYAcANwM6mPGoYCO4B47627VwLTYq/lnMM7H8iQgUGr\nmLvz0JgpEVckIhKOMI8cbgP+5u4f7uH1MmCFu+9Mmr8M6GBm7UOsJWuuG7a7Y/rRcdPYUZn89kRE\nck8o4WBmPYCLgDv3slh7YGOK+Rtij+3CqCXbzjnxGDod2BaA1es28OykORFXJCKSvrCOHO4Gfu3u\nqb7881pJSTFXn39ifHrk6EkRViMiEo60w8HMTgL6ACNTvZzwfC3QOsUybWKPOdube9X5J1JUFLzV\n/5v5Du8sWxVxRSIi6SkJYRunE5xxNNMsngUdY4/jzWwHwdlIc4B+ZlaS1O/QHVjl7mtTbXz48OHx\n5+Xl5ZSXl4dQcrg6d9ifc08qizcpPTR2Cr//0UURVyUihaCiooKKiorQt2uZuBeymd0O3A50c/fl\nsXmDgYnAKe4+KTavCbAKeNLdv59iO54r92qeMH0+Z/33vQDs36YFH4+/m+ax24qKiGSLmeHuqc4Y\nrZdMja1kSY+4+8sE4TDCzJrHZt8KVAJ3ZaiOrBkysDfdDwlOuFq/YQujXnkz4opERBou7CukzzKz\nt4FrCIbPGG9mbyUsMgyYC8w2s4XAQKDc3VeGWUcUioqKuOaCk+LTI8foimkRyV0ZaVYKSy41KwGs\n+XQDnc++icqduwCY/fefc+yRXSKuSkQKSWNvVipIBx3QhqGnfiE+/aCumBaRHKVwCFniPab/Nv51\nNm7etpelRUQaJ4VDyE7+Qg96dT8YgE1btvP3CTMirkhEpP4UDiEzsxpHDw+MnkQu9ZuIiIDCISO+\nefZ/0bxpKQBzlnzEjAXLoi1IRKSeFA4ZsH+blnxtyID4tMZbEpFco3DIkMSmpademsXazzZFWI2I\nSP0oHDJkQJ9uHNczuMZh2/ZK7nx0fMQViYjUncIhQ8yMX1z15fj0n0dV8N5Hn0RYkYhI3SkcMuj8\n8r4MKjscgMqdu7j1/mcjrkhEpG4UDhlkZvz2B7tvj/3USzOZqTOXRCQHKBwybNCxh9cYUuPG/xmt\n6x5EpNFTOGTBXd87n5LiYFdPemsJL0ydF3FFIiJ7p3DIgiO7duCaC3ef2vqze8ewMzZyq4hIY6Rw\nyJLbrv4yrVs2A2Dheyt57LnpEVckIrJnCocsOeiANvzsW2fEp28bOY7NW7dHWJGIyJ4pHLLoR18/\nnU4HtgVg1boN/OHJVyKuSEQkNYVDFrVo1oRfXnNufPqev05k9boNEVYkIpKawiHLLj93EH0O6wQE\n93u44+HnI65IRKQ2hUOWFRcXcc9/XxiffmjsFN5ZtirCikREalM4ROCsE47mlP49Adi1q4qb/zw2\n4opERGpSOETAzGocPYx9dTbTZr8bYUUiIjUpHCLSv3c3Ljlj9w2BbvyThtUQkcZD4RChO797Pk1K\nSwB4be57jPm/tyOuSEQkoHCIUPdD2nP9xeXx6ZvuG0ulhtUQkUZA4RCxW688m7atWwDw7odreHD0\n5IgrEhFROETugP1acuuVZ8Wn73j4eTZs2hphRSIiCodG4fqLT+HQjgcAsPazTdzz14kRVyQihU7h\n0Ag0a1rKnd89Lz79hydf4eM16yOsSEQKncKhkbj0zOM5rmcXALZur+S2kc9FXJGIFDKFQyNRVFRU\n437Tjz8/nXnvfhxhRSJSyBQOjchpx/fizEF9AKiqcm66d0zEFYlIoUo7HMysr5k9bGYLzWyumS0w\ns/8xs/ZJy7Uys/vMbHFsmYlm1jvd359v7v7+hZgZAOOnzef/Zi6OuCIRKURhHDk8BbQF+rl7GTAY\nGAJMM7NmCcuNAsqAvu7eB3gDqDCzTiHUkDfKenTmsnMGxqdv/J/RVFVVRViRiBSiMMKhCviZu28F\ncPcVwG+BHsDZAGY2GDgDuM3dt8XWGwEUA7eEUENeGXHteTRrWgrAW4uX89RLsyKuSEQKTRjhUObu\n7yXNWxl7bBt7HArsAKZWL+DulcC02GuSoHOH/fnRpafFp2/58zNs31EZYUUiUmjSDgd335li9pGA\nA9VjQZQBK1IsuwzokNw/IfCzy86k3X4tAfhg5True7oi2oJEpKCEfraSmRUD3wYecffqmxS0Bzam\nWLz6Bsrtwq4j1+3Xqjm3XX1OfPr2B59j7tKPIqxIRApJJk5l/QWwHfhhBrZdUK4dejI9u3YAYPPW\n7Zz3k/tZ+9mmiKsSkUJQEubGzOwKYBhQXt1BHbMW6JhilTaxx3V72ubw4cPjz8vLyykvL0+7zlzR\npLSE0fdcy8ArfsOmLdtZtmIdF9/0EBPv+wGlJcVRlycijUBFRQUVFRWhb9fCuvuYmX0T+Clwmruv\nSXptJHAF0DKx38HMniM4BTbl6axm5ro7GjxbMZvzb3ggPn39xeXc+9NLIqxIRBorM8PdLd3thNKs\nZGbfICkYzOwcM7s6tshooBQ4IWGdJrHp0WHUkM/OK+/LiGu/Ep++7+kKHnlm6l7WEBFJT9pHDmb2\ndeAR4OfA6oSXTiI4Q+mO2HIvAi2BM9x9q5ndAVxLcFHcSlLQkcNu7s7FNz3Ev/79FgClJcW8OvLH\nnND3iIgrE5HGJKwjhzDCYR3B9QzJxThwh7v/MrZcS+BugiuodwEfAj9090V72bbCIcHmrdsZdOU9\n8bOWOrRrw8wnbqZL7F4QIiKNJhwySeFQ27IVaxnwrV/Hz1rq1+tQpjx8I82bNYm4MhFpDBpVn4Nk\nT7dO7fnX3ddQUhz80725aDlX/epvKERFJEwKhxz0pX5H8qcbvxaf/vuEGfz2ry9FWJGI5BuFQ466\ndujJfOeCk+LTN903lvFT50VYkYjkE/U55LAdlTs57bo/MnV2MEpJm5bNmPHEzfTslup6QxEpBOpz\nkNgV1NfQpcP+AGzYvI2v/Ph+Ptu4JeLKRCTXKRxy3EEHtOHZ33+X5rH7PyxZvppLb32UXbt0gyAR\naTiFQx447qhDeez2y+LTL06fzy1/HhthRSKS6xQOeeKrQwZw8xVnxqfv+etL/H3CjAgrEpFcpg7p\nPFJVVcV5P7mf56cEZy01a1rKlIdvoH/vbtEWJiJZow5pqaWoqIj/HfFtjoqdrbRteyUX3DiSVWs/\nj7gyEck1Coc8s1+r5oz7w3dp27oFAB+tXs/Qnz6oe1CLSL0oHPJQj0M78NRdV1FUFBxZTp/7H753\n9z80xIaI1JnCIU+d8cU+3PPfQ+PTjz47jZ/9aQw7d+6KsCoRyRXqkM5j7s5ltz/O38a/Hp93Yt8j\n+MedV9E5duGciOQXDdktdRJ0Sj/AhOkL4vPa7deSv/3ySs464egIKxORTFA4SJ1VVVXxm8cn8IuR\n46iq2r0/f3bZGYy47jxKS4ojrE5EwqRwkHqb/NYSLrn1UVZ88ll83qCyw3nqrqt0NzmRPKFwkAb5\nZP1GvnnbX5j42sL4vAP2a8kTwy/nnJPKIqxMRMKgcJAGq6qq4p6/vsTPH3i2xgB9N3xzMHd97wI1\nM4nkMIWDpG3q7Hf52i0P8/Ga3c1MA4/pzlN3XU3Xg9tFWJmINJTCQUKx9rNNXHb7Y4yfNj8+b/82\nLXhi+OWce/KxEVYmIg2hcJDQVFVV8fv/fZmb//xMjWamH3/9dH59/QU0KS2JsDoRqQ+Fg4Ru+pz/\n8LVbHubD1evj847v041//vpqunVqH2FlIlJXCgfJiHWfbeLyOx6PD/sN0LZ1Cx67/TLOL+8bYWUi\nUhcKB8kYd+cPT77CTfeOYWdCM9P3LirnlivPotOBbSOsTkT2RuEgGff6vPf46s0Ps3zVp/F5xcVF\nXFDel+uGfYlT+vfELO3PoIiESOEgWfHp55u54o4nGDd5Tq3XjurWkWuHnsxl53wxfv8IEYmWwkGy\nxt0Z/e+3uO/pCia9taTW682blnLpmcdz3bAv0a9X1wgqFJFqCgeJxIL/rGDk6Ek88cLrbNy8rdbr\nx/fpxnXDvsRXB/enebMmEVQoUtgUDhKpTVu28fcJM7j/X5OYs+SjWq/v36YFV37lBK4dejJHdDko\nggpFCpPCQRoFd+f1ee9x/6hJPP3Km+yo3FlrmSEDe3PdsC9xzonHUKJxm0QyKifDwcwOAv4I9IvN\nmgf80N0/3sPyCocc8sn6jTw2bjojx0zm/Y/X1nq9c4f9ufCU4xjQuysDenejx6EHUVSkO9WKhCnn\nwsHMmgAzgcXAJbHZfwEGAce5++YU6ygcclBVVRUTX1vI/f+q4IWp89nTv2Gbls3o16trPCwG9OnG\noR0P0OmxImnIxXC4GngQOMzdl8XmdQA+Bm5y99+lWEfhkOOWrVjLQ2On8Mgz0/hk/cZ9Ln/g/q3p\n36srA/oEgdG/V1c6tt8vC5WK5IdcDIcJQE937540fy6w2d2/mGIdhUOe2L6jkn/PWMyMBcuYuXAZ\nMxd+UKewgKA5akDvrvTv1ZV+vbrS9eB2dGzXhv1aNddRhkiSXAyHFcBidz81af444FR3b5ViHYUD\nUFFRQXl5edRlhMrdWb7qU2YuWMasRR8wc+EyZi38gA0pTo+tYeMKaN0JgKZNSuhwQBs6tmtDx3b7\nBY/tE57H5ndo14YWeXhabT5+LhpK+2K3sMIhm2MxtwdS/am4AWhhZk3dfXsW68kZ+fjBNzO6HtyO\nrge3Y9jpwfkJVVVVLF2+Jh4WMxd8wNvvLGfr9srdKyaEw/YdO1m+6tMaw3vsSeuWzeKB0b5tK1o0\naxL/ad60SdJ0acrnics2KS2mpLiYkuIiiouLIjmCycfPRUNpX4Qvm+GgQwDZq6KiInp260jPbh35\n+ln/BcDOnbtY+P5KZi4ImqLG/+s9mnY5iFXrPmfTlrr/LbFx8zY2bt7G0uVrMlS7UVJcTHGRUVIS\nC42ionh4xIOkqIiSkuCxqMgwDLNgfbPd02ZW43Wz2HJWFJ9+f9ZrTF3zRwCM3eFUnVPVgZWYW/F5\nKZZvqDCCMd1NvPPaDN78/L6065DdshkOa4HWKea3Iehz0FGD1FJSUkxZj86U9ejMt88/keE7ljJ8\n+HAguBBv9acbWbX2c1at28DqTzewat0GVq37nFVrE56v20Dlzl0ZrbOqytlRFbvGI/FIJ5NWrGXZ\njMXZ+V2N3YrVLN05b9/LSZ1ls8/hReCoFB3S84CN7j4oxTo62hARqadc63MYAzxoZl3d/QOIn8p6\nFHBTqhXCeIMiIlJ/2TxyKAVmAYuArxP0QTzK7ovgtmSlEBER2aesjV3g7pXAYGAXsDD204rgNFYF\ng4hII9KoB94TqS8z+xVwC3CFuz8RdT0iYTGzg4HHgCHunvE/7LM+6pmZHWRmT5rZ4tjPKDM7pI7r\nlprZCDNbZGbzzGyamZ2Q6ZozpaH7wswONrM7zGxObD8sMrPRZnZ0NurOhHQ+Fwnb6Az8mKDJMmf/\n6kl3X5jZsWb2rJm9GftsLDazuzNZc6ak+X1xsJk9EtsHc8xsvpndbGbZ7GsNhZldCEwHDqOen+0G\nf2+6e9Z+gCbAHOCfBMFUBDwOLAFa1mH9kQQD97WLTX8b2Awcm833EfW+iO2Hd4BDYtNNgadj++Lo\nqN9btj8XCdv5K/AcUAV8K+r3FcW+IOjD+xj4YsK87wLvRf3esrkvYsu+DcwF9o/N6wtsAX4b9Xtr\nwL54HTg89v6r6rlug743s/0Gr479x+2WMK8DsBO4YR/r9iTor7g8af584Pmo//GyvC8eAK5MmndY\nbHt/ivqkEMQQAAAEZElEQVS9ZXNfJCzfD3gXGJLj4ZDO58IITvj4SdL8EuCMqN9blvdF79i6P0ia\n/wywIur31oB9URR7rFc4pPO9me1mpaHABx4blRXA3VcTdE4P3ce6FxB8+F9Nmv8qMMTMcu0O9+ns\ni+sJ2h4TrYw9tg2rwCxKZ19U+z1BX8OO0KvLrnT2xYkEXwbPJ850953uPjHkOrMhnX1Rfdep0qT5\npUTQnJ4ud69q4KoN/t7M9k4qA95PMX8ZcEwd1t0FLE+a/z7BX0a90y0uyxq8L9x9l8fiP8GRsceK\ntCvLvnQ+F5jZ+UBTd3865LqikM6+qL6QtG2sz2F+rK19hJk1C7PILEnn/8gS4O/ANWbWFcDMTgVO\nB+4Nt8xGrcHfm9numEln8L32wJYUX4obYo/tQqoxW8IeiPA7BIeKfwujuCxr8L6IXT/zG+DKDNaX\nTel8LrrEHv8BXOzus8ysDHgBOB44I/RqMyvd/yOXEdx5cqmZfUJw6vwP3H1k+KU2Wg3+3sz2kUPO\nnkGSAaHtCzM7DbiY4AshSwP7hCqdfXEdMN/dp4dVTMTS2RfVRwePuPssAHefC9wNDDazk9MtLssa\nvC9iR0oVQH+gq7sfApQDt5jZLaFUl+eyHQ7pDL63FmhptYeAbBN7XBdCfdkUykCEZnYsQSfVue6e\nq6OwNWhfmFlbgqFXUg2/kqtDr6Tzuaj+K3t20vzq6f5p1pZt6eyLK4ETgBvdfSWAu78N/A4YEft/\nUwga/L2Z7XCYC3RPMb87sK8hFecQ1NslaX53oJKgkyqXpLMvAIg1GYwFvurur4dYW7Y1dF8MJOh4\nHGVmb5vZ28DDsdd+GZv383BLzbh0PheLYo/J/6937WF+Y5fOvqjuk1iaNH8pwR8OuRaUDdXg781s\nf1jGAF2rO4igxuB7oxMXNLMOSWk3luAw85SkbZ4CvOS5NwRHOvuiOhieAb5R3aQSu+gnF9tTG7Qv\n3H2Cux/q7sdV/wBXxRb9RWzer7L0HsKSzudiPEEQJP9VXH1x5Mzwy82odPbF6thjV2qqns61loZE\ne2xuC/V7M8vn6pYSJNlTQDFBOD1GcEFXi4TlTiD4kN+ftP4D1LyY4wqCiznKsvk+ot4XBH8VfRLb\nH99I+Pkh8GrU7y3bn4ukbZUTnN9+WdTvK4p9QXBK7wrgiNj0IQR/LU+I+r1lc18A3YDPgYlAq9i8\nQwmuhVlKcHZb5O+xAfvkcfZwnUPY35tZPXLwug++txFYT/AhT/R9YBQwzYL7QFxFMM7I3EzXHrY0\n98Vw4ADgGoKrgqt//kAOdvqH8LmoHmahulnJ2d2s1C/T9YcphH1xI8GpmuPNbBEwieCv7PMyXHro\n0tkXHlwbcTzwKTDTzOYAE4AXgRM8x24uZmb3xD7f5wIe+2y/FTtbr1qo35saeE9ERGrJtQ4qERHJ\nAoWDiIjUonAQEZFaFA4iIlKLwkFERGpROIiISC0KBxERqUXhICIitSgcRESkFoWDiIjU8v8BOf+N\ndB9FDUEAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f20761b0310>"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T1 = T_analytical(x, .2, 100, alpha, L)\n",
      "T2 = T_analytical(x, .2, 200, alpha, L)\n",
      "numpy.sqrt(numpy.sum((T1-T2)**2)/numpy.sum(T2**2))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "6.9279171182600926e-13"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That looks like it should. We'll now use this result to study the convergence of the Crank-Nicolson scheme."
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Time convergence"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We said this method was second-order accurate in time, remember? That's in theory, but we should test that the numerical solution indeed behaves like the theory says.\n",
      "\n",
      "Leaving $\\Delta x$ constant, we'll run the code for different values of $\\Delta t$ and compare the result at the same physical time, say $t=n_t\\cdot\\Delta t=10$, with the analytical expression above.\n",
      "\n",
      "The initial condition of the rod problem has a very sharp gradient: it suddendly jumps from $0{\\rm C}$ to $100{\\rm C}$ at the boundary. To resolve that gradient to the point that it doesn't affect time convergence, we would need a very fine mesh, and computations would be very slow. To avoid this issue, we will start from $t=1$ rather than starting from $t=0$.\n",
      "\n",
      "First, let's define a function that will compute the $L_2$-norm of the error:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def L2_error(T, T_exact):\n",
      "    \"\"\"Computes L2 norm of error\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    T      : array of float\n",
      "        array with numerical solution\n",
      "    T_exact: array of float\n",
      "        array with exact solution\n",
      "    Returns:\n",
      "    -------\n",
      "    e: L2 norm of error\n",
      "    \"\"\"\n",
      "    e = numpy.sqrt(numpy.sum((T-T_exact)**2)/numpy.sum(T_exact)**2)\n",
      "    return e"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For fun, let's compare the Crank-Nicolson schem with the implicit (a.k.a., backward) Euler scheme. We'll borrow some functions from [notebook 2](http://nbviewer.ipython.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/04_spreadout/02_Heat_Equation_1D_Implicit.ipynb) to do this."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def generateMatrix_btcs(N, sigma):\n",
      "    \"\"\" Computes the matrix for the diffusion equation with backward Euler\n",
      "        Dirichlet condition at i=0, Neumann at i=-1\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    T: array of float\n",
      "        Temperature at current time step\n",
      "    sigma: float \n",
      "        alpha*dt/dx^2\n",
      "    \n",
      "    Returns:\n",
      "    -------\n",
      "    A: 2D numpy array of float\n",
      "        Matrix for diffusion equation\n",
      "    \"\"\"\n",
      "   \n",
      "    # Setup the diagonal\n",
      "    d = numpy.diag(numpy.ones(N-2)*(2+1./sigma))\n",
      "    \n",
      "    # Consider Neumann BC\n",
      "    d[-1,-1] = 1+1./sigma\n",
      "    \n",
      "    # Setup upper diagonal\n",
      "    ud = numpy.diag(numpy.ones(N-3)*-1, 1)\n",
      "    \n",
      "    # Setup lower diagonal\n",
      "    ld = numpy.diag(numpy.ones(N-3)*-1, -1)\n",
      "    \n",
      "    A = d + ud + ld\n",
      "    \n",
      "    return A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def generateRHS_btcs(T, sigma):\n",
      "    \"\"\" Computes right-hand side of linear system for diffusion equation\n",
      "        with backward Euler\n",
      "    \n",
      "    Parameters:\n",
      "    ----------\n",
      "    T: array of float\n",
      "        Temperature at current time step\n",
      "    sigma: float\n",
      "        alpha*dt/dx^2\n",
      "    \n",
      "    Returns:\n",
      "    -------\n",
      "    b: array of float\n",
      "        Right-hand side of diffusion equation with backward Euler\n",
      "    \"\"\"\n",
      "    b = numpy.zeros_like(T)\n",
      "    \n",
      "    b = T[1:-1]*1./sigma\n",
      "    # Consider Dirichlet BC\n",
      "    b[0] += T[0]\n",
      "    \n",
      "    return b"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def implicit_btcs(T, A, nt, sigma):\n",
      "    \"\"\" Advances diffusion equation in time with implicit central scheme\n",
      "   \n",
      "    Parameters:\n",
      "    ----------\n",
      "    T: array of float\n",
      "        initial temperature profile\n",
      "    A: 2D array of float\n",
      "        Matrix with discretized diffusion equation\n",
      "    nt: int\n",
      "        number of time steps\n",
      "    sigma: float\n",
      "        alpha*td/dx^2\n",
      "        \n",
      "    Returns:\n",
      "    -------\n",
      "    T: array of floats\n",
      "        temperature profile after nt time steps\n",
      "    \"\"\"\n",
      "    \n",
      "    for t in range(nt):\n",
      "        Tn = T.copy()\n",
      "        b = generateRHS_btcs(Tn, sigma)\n",
      "        # Use numpy.linalg.solve\n",
      "        T_interior = solve(A,b)\n",
      "        T[1:-1] = T_interior\n",
      "        # Enforce Neumann BC (Dirichlet is enforced automatically)\n",
      "        T[-1] = T[-2]\n",
      "\n",
      "    return T"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now, let's do the runs!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx = 1001\n",
      "dx = L/(nx-1)\n",
      "\n",
      "dt_values = numpy.asarray([1.0, 0.5, 0.25, 0.125])\n",
      "error = numpy.zeros(len(dt_values))\n",
      "error_ftcs = numpy.zeros(len(dt_values))\n",
      "\n",
      "t_final = 10.0 \n",
      "t_initial = 1.0\n",
      "\n",
      "x = numpy.linspace(0,L,nx)\n",
      "\n",
      "Ti =  T_analytical(x, t_initial, 100, alpha, L)\n",
      "T_exact = T_analytical(x, t_final, 100, alpha, L)\n",
      "\n",
      "for i,dt in enumerate(dt_values):\n",
      "    sigma = alpha*dt/(dx*dx)\n",
      "\n",
      "    nt = int((t_final-t_initial)/dt)\n",
      "    \n",
      "    A = generateMatrix(nx, sigma)\n",
      "    \n",
      "    A_btcs = generateMatrix_btcs(nx, sigma)\n",
      "\n",
      "    T = Ti.copy()\n",
      "    T = CrankNicolson(T, A, nt, sigma)\n",
      "    \n",
      "    error[i] = L2_error(T,T_exact)\n",
      "    \n",
      "    T = Ti.copy()\n",
      "    T = implicit_btcs(T, A_btcs, nt, sigma)\n",
      "    \n",
      "    error_ftcs[i] = L2_error(T,T_exact)\n",
      "    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And plot,"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(6,6))\n",
      "plt.grid(True)\n",
      "plt.xlabel(r'$\\Delta t$', fontsize=18)\n",
      "plt.ylabel(r'$L_2$-norm of the error', fontsize=18)\n",
      "plt.xlim(1e-2,1)\n",
      "plt.ylim(1e-4,1)\n",
      "plt.axis('equal')\n",
      "plt.loglog(dt_values, error, color='k', ls='--', lw=2, marker='o')\n",
      "plt.loglog(dt_values, error_ftcs, color='k', ls='--', lw=2, marker='s');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAGXCAYAAAANqX+7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XucFOWV//HPAbkoXlARwesoId4vERIjiTKyKlE2xvsl\nroI/TXbdDchPUYxRGTbGFTUqkMTNRRnXn4mLt9WIEiQwamJMNiqCGoMSRoOKeCMS5DI65/dH91Az\nQ/dMX6q6qru+79erX9pV1dVnzquYM/U8Tz2PuTsiIiJJ1iPuAERERLqjYiUiIomnYiUiIomnYiUi\nIomnYiUiIomnYiUiIom3RdwBVJqZ7QDcBPwdcGBv4BJ3fzXWwEREJK/UFStgd+Bjd58AYGbfAu4A\njow1KhERySt1xcrdXwC+1W7TcmDXmMIREZECJLrPyswGm9lcM2uN8Gu+CvwgwvOLiEiZEluszOwU\n4GkyfUp554Qys4FmdreZvZJ93Wtmu7bbP8HMXs2+RnX67Bhga3e/OaqfQ0REypfYYgVcDhxDpmBZ\nrgPMrDfwOJnmzP2zr7XAQjPrB+DuM9x9aPa1oN1nx5C5qzo30p9CRETKluRiNcLdl3VzzFjgIGCy\nu7e6eyswmczd2EX5PmRmpwPHuvu/uLub2fTQohYRkdBZ0mddN7NG4Dx336ywmtlcYB9336vT9sXA\nWnc/IsdnDgaeBd4luGPb1t37hR27iIiEo9pHAx4MvJJjezMwKsd23H0x0CvCmEREJGRJbgYsxABg\nTY7tHwFbmVmfCscjIiIRqPY7q0jbMM0s2W2kIiIJ5e45B8aVqtrvrN4DtsmxfVsyfVYbyv0Cdw/1\nNWXKlFCPzXdModu7ep/v/5UL5aJWc1HItkrlotjzJSkXUaj2YrUY2CvH9r2AJRWOpSD19fWhHpvv\nmEK3d/W+/f83Nzd3G0uxlIv8313uscpF98fk2l7Itkrlopg8FHp8pXIRibD/Kgr7BTQCn+bZ9w2g\nFdiz3badgRbg0hC+26dMmeILFy70tBs7dmzcISSGchFQLgLKhfvChQt9ypQpnikt4daCah+63gv4\nI/An4BwyfVi3AyOAz7n7x2V+tyc9P5XS1NQU/V9OVUK5CCgXAeUiYGZ4yH1WiS1WZnYDcCywB9Cf\nTJOfA4e7e0u74wYCtwDDs/uXABPd/c0QYlCxEhEpUhTFKrF9Vu5+ubt/zt13dPee2f8/rH2hyh63\nyt3Pcfd93H1fdz89jELVpqGhgaamprBOV7WUg4ByEVAuAspFJgcNDQ2RnLvah65HLqrEi4jUmvr6\neurr65k6dWro505sM2ASqBlQRKR4qWoGFBERaaNi1Q31WWUoBwHlIqBcBJQL9VnFSn1WIiKFUZ9V\nTNRnJSJSPPVZiYhIKqlYSUHUHh9QLgLKRUC5iJaKVTc0wEJEpDBRDrBQn1UX1GclIlI89VmJiEgq\nqVhJQdQUGlAuAspFQLmIloqViIgknvqsumBmPmXKlE0PuomISH5NTU00NTUxderU9KxnlQQaYCEi\nUjwNsJDYqD0+oFwElIuAchEtFSsREUk8NQN2Qc2AIiLFUzOgiIikkoqVFETt8QHlIqBcBJSLaKlY\ndUNzA4qIFEZzA8ZEfVYiIsVTn5WIiKSSipUURE2hAeUioFwElItoqViJiEjiqc+qC+qzEhEpnvqs\nREQklVSspCBqjw8oFwHlIqBcREvFqht6zkpEpDB6ziom6rMSESme+qxERCSVVKykIGoKDSgXAeUi\noFxES8VKREQST31WXVCflYhI8dRnJSIiqaRiJQVRe3xAuQgoFwHlIloqViIiknjqs+qC+qxERIqn\nPisREUklFSspiNrjA8pFQLkIKBfRUrHqhuYGFBEpjOYGjIn6rEREiqc+KxERSSUVKymImkIDykVA\nuQgoF9FSsRIRkcRTn1UX1GclIlI89VmJiEgqqVhJQdQeH1AuAspFQLmIloqViIgknvqsuqA+KxGR\n4kXRZ7VFmCerFmZ2K7AN8CFwKDDL3e+ONyoREcknrc2AG939AnefBFwJ/MzM0pqLgqg9PqBcBJSL\ngHIRrVT+gnb3y9u93Rd4yd1b44pHRES6lug+KzMbDMwCjnP3UAurmR0MXAXsDZzu7stzHKM+K5FO\nxo0bR3Nz82bb6+rqaGxsrHg8kjyp6rMys1OA7wMtQN6KYWYDgVuAYdlNS4CJ7v5mdv8EYHx23z+7\n+wIAd18MnGFmhwNPmtkh7v5BJD+MSA1pbm7miSeeiDsMSZkkNwNeDhwDPA3krNBm1ht4nEzR3T/7\nWgssNLN+AO4+w92HZl8LzKxH277s/t8D64CjIv1pqpza4wPKheSi6yJaSS5WI9x9WTfHjAUOAia7\ne2u232kymaa9i/J8Zg/gZ21vzGxHYBDQ3XeJiEhMEtsMWOCAh1OB1929ud3n3jGzl7P7bsrxmfcB\nM7NZwAfAfsDF7r6k/KhrV319fdwhJEbac7F27dq4Q0iktF8XUUtssSrQwcArObY3A6NyfcDd1wBn\nRRiTSM164YUXWLRoUdxhSAoluRmwEAOANTm2fwRsZWZ9KhxPzVJ7fCCtuXjuuecYNWoUn3zyCdtv\nvz1HHnkkhxxyCCNHjmTkyJHU1dXFHWKs0npdVEq131lFPq583Lhxm/4R9u/fn0MPPXTT7X7bxan3\n6XrfJinxVOr9c889R0tLCyeeeCKzZ8/md7/7HYsWLWLixImbjm9qakpMvJV+33bHmZR4Kvm+qalp\n02MLUf3RkujnrADMrBE4L9dzVmb2JvBndx/VafvDwNHuvk2Z363nrETaWbp0KXV1dfTu3TvuUCTB\nUvWcVYEWk5mBorO9yDxvJSIh+uxnPxt3CJJS1dJnle/25gFgTzPbs22Dme1MpoDdH8YXNzQ0bNb8\nk0bKQUC5CCgXAeUik4OGhoZIzl0txSrf7WQjmTuoaWbWMzsZ7fXAX4DbwvjihoaGTW20Imkxf/58\nbrstlH9CkiL19fWRFavE9lmZ2Q3AsWQe4u1PpsnPgcPdvaXdcW3TLQ3P7u8w3VKZMajPSlJn7ty5\nnHTSSWzYsIGFCxfqjzUpWqr6rDrNjN7VcauAc6KKo+3OSv9gJQ0eeeQRTj31VDZu3MhFF13EUUdp\nFjIpXNuI0Cgk9s4qCXRnFWg/JDntajUXDz74IGeeeSYtLS2MHz+e6dOnY9b1H8e1motSKBeBKO6s\niu6zMrNZZnaHmR0TZiAiEp9169Yxfvx4WlpauOSSSwoqVCKVVPSdlZm1AncD09z9xUiiSgjdWUma\nvPjii/zP//wP3/nOd1SopCxJ6bNa5e7nhhlEkqnPStLiwAMP5MADD4w7DKliieqzMrP5ZGaUeKuL\nY+5097HlBhc33VkF1B4fUC4CykVAuQgkos8KmADcZmbDuzjmuBLjEZEKWLx4cdwhiBSllDur5cDW\nwI5kVth9F2i/9pQBu7l7r7CCjIvurKQWjBs3jubm5k3v33zzTV577TU+//nP84c//CG+wKRmJaXP\nahe6WGo+a3Bp4YhI2Jqbm3niiSc22/7BBx/EEI1IaUppBvzA3Y929/p8L+DDkOOMjeYGzFAOArWS\ni912263sc9RKLsKgXEQ7N2Apd1aFjAT8SgnnTaSoEi8iUmvaRk5PnTo19HNrBosuqM9KakF9fX3O\nZsCRI0fqbkAikZQ+KyzzxOB5wBnAkOzmV4HZ7n5XSLGJiIgAJRQrM+sNPASM7rTrs8AYM/s6cGL7\nmdGl+ukZkkC15SLfMuNhLD9ebbmIknIRrVLurK4APgdcBjwMvJ3dvgtwIjAJ+Dbw72EEGDfNYCHV\nrrGxMe4QJCWSNoPFUuAsd38uz/5hwD3uPjSE+GKlPiupNu7OVVddRb9+/bjyyivjDkdSKoo+q1KK\n1Up3H1TuMdVAxUqqibvz7W9/m2nTptGzZ09eeukl9tlnn7jDkhRKynRLLWa2c76dZjYI+KT0kCSJ\nNGoskMRcuDuXX34506ZNY4sttuC///u/K1KokpiLuCgX0SqlWM0F7jOzwzrvyDYB3gc8Vm5gIlIY\nd+fSSy/lpptuYosttmD27NmceuqpcYclEqpSmgEHAb8HdgdW0nGAxc7AG8AX3X1liHHGQs2AUg3e\ne+89DjvsMFauXMl9993HiSeeGHdIknKJ6LPKBjIQuA44Ddg2u/kj4F7gSnd/N7QIY6RiJdVi2bJl\nLF26lOOPPz7uUEQS02eFu69y9wuBHchMWjsY2MHdv1ErhaqN5gbMUA4CSczFkCFDYilUScxFXJSL\nhM0NaGazAAd+7u7zgXdCjypBNDegiEhhEjU3oJm1AncD09z9xdAjShA1A0rStLa28sgjj6hfShIt\nKc2Aq9z93FovVCJJ09rayje/+U2+9rWv8b3vfS/ucEQqqpRi9aKZ7dLVAWZ2Z4nxSEKpPT4QRy4+\n/fRTLrjgAm6//Xa23HJLjjjiiIrHkIuui4ByEa1SitUE4DYzG97FMceVGI+IdPLpp59y/vnn09jY\nyFZbbcWjjz7KqFGj4g5LpKJK6bNaDmwN7AisA94FWtsfAuzm7r3CCjIu6rOSOIwbN47m5uZN75ct\nW8aKFSvo2bMnCxYs4KijjoovOJECJGU9q12B35IpSvkMLi0cEWlubs65WOKBBx6oQiWpVUoz4Pvu\nfrS71+d7AR+GHGds9JxVhnIQiCsX/fv3j+V7u6LrIqBcRPucVSnF6hMzu8PMjunimK+UGlDStK1n\nJSIiXauvr09UsdoV6EVmXsCc3P2FkiOSRFLBDkSZi08//ZRq6ifVdRFQLqJVSp/VKnc/N/RIRFKu\npaWFs88+m+XLl8cdikji6DkrKYja4wNR5KKlpYWzzjqL+++/n5UrV3L44YczcuTIDq+6urrQv7dc\nui4CykW0SrmzanvO6rvu/sc8x+g5K5ECbdy4kbPOOosHH3yQ/v378/jjjzN8eFePMYqkj56z6oKe\ns5Kobdy4kTPOOIOHHnqI/v37M3/+fIYNGxZ3WCJlScpzVrsAT6PnrETK9uGHH/Lyyy+z/fbbM3/+\nfA47bLMFuEWE0vqsPkjTc1aSofb4QJi52HnnnVm4cCELFiyoykKl6yKgXESrlDurQkYC1sxzViJR\n23XXXdl1113jDkMk0Upa1j4t1GclIlK8pKxnhWWcbGY/NrP7s9uGZrf1DDNAkVqxfv16fvrTn1bV\nQ78iSVF0sTKzLYFfA/cD3wCOze7aGvgZMN/MtgktwphpbsAM5SBQSi7Wr1/PySefzDe/+c3IpqOJ\ng66LgHKRvLkBrwb2Af4NOBz4GMDdnwd2A94HrgwrwLhpbkAp1/r16znppJOYO3cuAwYM4LTTTos7\nJJFIRDk3YCnPWb0K/JO7/z77/m13H9xu/wDgaXf/bKiRxkB9VlKudevWcdJJJzFv3jx22mknFixY\nwIEHHhh3WCKRSkqf1Q5thSoXd3+PTJOgSOpdeumlzJs3j4EDB7Jw4UIVKpESlVKsNppZ3od+zWwI\nHWe0kBqg9vhAMbm45pprqK+vZ+HChRxwwAHRBRUTXRcB5SJapTxn9Sgw28zOd/fX2u8wsy8DM4Ff\nhhGcSLUbNGgQCxcujDsMkapXSp/VYOAZYHdgOZlBFS+RWedqJ6AZ+KK7rwo10hioz0oKMW7cOJqb\nmzfbXldXR2NjY8XjEYlbIuYGdPe3zezzwHXAaWQWYjwU+BvwU+A72X4rkVRobm7miSeeiDsMkZpW\n0kPB7r7K3S8EdiAzae1gYEd3/2cVqtqk9viAchFQLgLKRbRK6bPaxN1bgXdCikVERCSnku6saoWZ\nTTIzjVwsgB6MDigXAeUioFxEK7XFyswOBOoBjaAQEUm4spoBq5WZ9QK+C3wbOCHmcKpCU1OT/nLM\n6pyLurq6nMfl215LdF0ElItoJbpYZYfJzwKOc/cw7wKnANOBNSGeU1JKw9NFopfY9azM7BTg+0AL\nMMTdcy49YmYDgVuAYdlNS4CJ7v5mdv8EYHx2378A64DT3P0SM6sD/pKvEOo5KxGR4iVlbsBKrWd1\nOXAM8DSQ84c2s97A42TuEPfPvtYCC82sH4C7z3D3odnXr4ETge3N7Dbg2ux5fmRmmgpbRCShkrye\n1Qh3X9bNMWOBg4DJ7t6aHUo/GdgbuCjXB9z9Cnc/390vAq7KbvtXd78vhJhrlp4hCSgXAeUioFxE\nK7HrWWULT3dOBV539+Z2n3sHeDm7Ly8zGwk0AG5mM83soNKjFRGRKCV+PSszawTOy9WvZGZvAa+4\n+6hO2x8GRrl7WUuVqM9KRKR4SemzStJ6VgPIPaLvI2ArM+tToThERCRCpQxd32hmg9397Vw7K7ye\nVeS3PePGjdv0vEz//v059NBDNz1L0dZGnYb37dvjkxBPnO/btiUlnjjfL1q0iIkTJyYmnjjf33rr\nran+/dD2CEdUzxeW0gx4O/BZ4Hx3f619M2C79ayeyQ5gKD/ArpsB3wT+nKcZ8Gh3L2ugh5oBA016\n4HET5SKgXASUi0AUzYCJX8+qm2L1GLCvu+/VafsSYI27jyjzu1WsRESKlIg+q2zz3+eBO4AdCdaz\n6k1mPasvRLDwYr6K8QCwp5nt2bbBzHYG9iUztL5sDQ0NHZp/REQkt6amJhoaGiI5d1kzWJhZDzJ3\nUwDvFjjcvNjvaCT/nVUv4I/An4BzyBS124ERwOfc/eMyv1t3Vllq4ggoFwHlIqBcBBJxZ9Ve9kHc\nd7KvTYXKzPYrNzAzu8HMnge+SuZZqOfN7LlsgWr7/hYyDyV/SubZqpfJjEQcVW6hEhGR5IhkbkAz\ne8vddwn9xBVmZj5lyhTq6+v1F5OISDeamppoampi6tSp8Q+wADCzk4CvAbuQ6bPqsJvMVElV/4yT\nmgFFRIqXiGZAM5tIZmDDOWTm5du702uvUs4ryaZBJgHlIqBcBJSLaJXyUPC3gCuAW919Y64DzCzn\nA8MiIiKlKOU5q3eBgV21j5lZvbs3lRlb7NRnJSJSuET1WZlZE3Cqu7/fxTGj3H1BmbHFTn1WIiLF\nS0SfFZlmwB+Z2SFdHHN3ifFIQqk9PqBcBJSLgHIRrW77rMxsOZvPILENcLqZfQy8R8eJa43MbOgi\nIiKh6LYZ0Mw20MXS8nl80d37lhNYEqgZUESkeFE0AxYyGvADdz+6mJPW0mjAhoYGDbAQESlA2wCL\nKBRyZ3WMu88v6GRmW7n7x2Z2sLsvDiXCGOnOKqB5zwLKRUC5CCgXgVgGWHQuVGb2oy4Ov8HM3iOz\nfIiIiEgoShm6vmmxxRz7BgKjgSvdvezJbOOmOysRkeIlZfHFvMUqu9+Ad9x9YLnBxU3FSkSkeHEN\nsMDMZpEZvm7AdmZ2Rxfn2xd4K5zw4qcBFhlqjw8oFwHlIqBcRDvAotC5Acd2ej8uz3Efk1kI8aJS\nA0qaqFa9FBGpNW1/2E+dOjX0c4feDFhL1AwoIlK8pEy3dGWYAYiIiHSn6GLl7rOiCESSTfOeBZSL\ngHIRUC6ipUUSRaRoc+bMYfTo0UycOJHRo0czZ86cuEOSGlfK4oupotGAGWn/+dtLey7mzJnDxRdf\nzLJlyzZta/v/MWPGxBVW7NJ+XUDM0y2lmQZYiGxu9OjRzJs3L+f2uXPnxhCRJE1SBlhICqk9PpD2\nXGzYsCHn9vXr11c4kmRJ+3URtW6LlZldZWYLzGzHSgQkIsnWp0+fnNv79q36VYEkwQq5szof+Cnw\nIYCZdX5AWFJA7fGBtOdi7NixZGZVCwwZMoTx48fHFFEypP26iFohAyz6ufsv2r2/Hrizqw+Y2SJ3\nP7SsyEQkkZYsWYK7M2DAAA444AD69u3L+PHjUz24QqJXyHpWfwHOcfffZd93O4NFrcxyoQEWAc17\nFkhzLtauXcvgwYNZs2YNzzzzDOvWrUttLjpL83XRWVwT2d4N/NbM3gHWAztlC1g+BgwIIzgRSZZ+\n/frxm9/8hkceeYTDDz9cgwqkYgq5s+oB/CvwD8D2wBHA77o57xfdvep7W3VnJSJSvFjurNy9FfhB\n9tXWxFff1WfM7O1QoksAPRQsIlKYRD0UbGbHdF7qvpRjqoHurAJqjw8oFwHlIqBcBGJbfLE9d5+f\nXQ34POAMYEh216vAbHe/qxYKlYiIJEcpd1a9gYeA0XkO+RVworu3lBlb7HRnJQKrV6/miiuuYNKk\nSXzmM5+JOxypAkmZbukK4HPAZcA+wLbZ177A5dl93w4rQBGJ14wZM/jxj3/MRRfVzALgUoVKKVb/\nBJzg7t9391fd/e/Z11J3vwkYA5wbbpgSNw1RDqQpF3/729+45ZZbALjqqqs225+mXHRHuYhWKcVq\nW3d/Lt9Od38W2Kb0kEQkKWbOnMnq1asZOXIkI0eOjDscSbFS+qz+Cgx393fy7B8E/NHddwshvlip\nz0rS7KOPPqKuro4PP/yQBQsWcPTRR8cdklSJpPRZzQXuM7PDOu8ws2HAfcBj5QYmIvH67W9/y9//\n/neOPPJIDcmW2JVSrK4G9gD+aGZvmdmz2dfbwP8Cu2WPkRqi9vhAWnJx/PHHs2zZMm677bbNZllv\nk5ZcFEK5iFYpz1mtNLPPA9cBp5EZ/QfwEXA7cKW7vxteiCISl9133z3uEESAMpe1z84buFP27bvZ\nqZlqhvqsRESKl4gZLNrLFqecAy1qheYGFBEpTKLmBkwT3VkFNO9ZoJZz4e55+6dyqeVcFEu5CCRl\nNKCI1Kibb76ZE044gRdeeCHuUEQ60J1VF3RnJWny8ccfs9dee7Fq1SrmzJnDCSecEHdIUqV0ZyUi\nkfnP//xPVq1axfDhwzn++OPjDkekAxUrKYieIQnUYi7WrVvHDTfcAMCUKVMK7reqxVyUSrmIViTF\nysyujeK8IhKNn/zkJ7zzzjsMGzaMMWPGxB2OyGbKfc5qALBV583A/7r7wHICSwL1WUlazJgxg6uv\nvpq77rqLE088Me5wpMpF0WdVykS2A4DpwClAnzyHubv3LDO22KlYSZqsXr2a7bbbrqih6yK5JGWA\nxU+ArwCzgeuBf8/x+jisAKNgZo1m9na71z1xx5R0ao8P1Gou+vfvX3ShqtVclEK5iFYpM1gcDXzB\n3V/Nd4CZJb0dwd19cNxBiIhIYUppBlzi7gdFFE9FmNksYBXQk8zd5U3u/laO49QMKCJSpKQ0A15v\nZt/o6gAz+2OJ8XQ+z2Azm2tmYU+Q+zBwi7tPIrP21m/MbMuQv0Mk0VasWIH+GJNqUXSxcve7gQ1m\n9jsz+4GZNZjZNe1eU4B9yg3MzE4Bngb2BvL+izKzgWZ2t5m9kn3da2a7tts/wcxezb5GZX+GB919\nZfb/HyfTHPrlcmOuZWqPD9RCLjZs2MARRxzB4YcfzsqVK0s+Ty3kIizKRbSK7rMysxOAH5MZCXh4\nnsPC+HPtcuAYMgs5fiZPLL2Bx4FXgP2zm+8AFprZ59x9rbvPAGZ0+txn3X1pu00bgL4hxCxSFRob\nG1mxYgXbbbcdAwdW/VMmkgKl9Fm9CPwJ+AHwNrAxx2F/KPc5KzPr4e6tZtYInOfum90FZpsjfwzs\n7e7N2W07A28CV7j7TXnO/Qd3/0L2//cGngX2cfdVnY5Tn5XUnI0bNzJ06FDeeOMN7rnnHs4888y4\nQ5Iak5T1rHYBDnX3T/IdYGY5i0QxClzI8VTg9bZClf3cO2b2cnZfvjiWmNnPgZVk7tpO71yoRGrV\nnXfeyRtvvMF+++3HaaedFnc4IgUpZYDFEqC7wQhPl3DeUhwMLM+xvRnIO2LR3S9w96+7+yXufqK7\nz48qwFqh9vhANeeipaWF6667DoCrr76anj3Le3a/mnMRNuUiWqUUq38DfmBm+3dxzOwS4ynWAGBN\nju0fAVuZWb4ZNkRSyd2ZNGkSxxxzDGeccUbc4YgUrJQ+q+XA1sCOwFrgfaB9k50Bu7t7KU2Mub6v\nkfx9VhuAue7+tU7b/x/wdWBLd99Qxnerz0pqwpw5c5gxYwYbNmygT58+TJgwQRPWSmSS1Gf1NJmi\nlE+lZod4D9gmx/ZtgbXlFKo248aNo66uDshMR3PooYduWrq67bZf7/U+ye/Xrl3LxRdfzLJly2iz\nbNkyFi9ezBFHHBF7fHpf/e+bmppobGwE2PT7Mmyl3Fm93d1URYUcU8T3NZL/zuoxYF9336vT9iXA\nGncfUeZ3684qq6mpadNFmnbVlovRo0czb968nNvnzp1b1rmrLRdRUi4CSZnB4p8KOOYrJZy3K/kq\nxgPAnma2Z9uG7ND1fYH7Q45BpCpt2JC7gWH9+vUVjkSkdCUVKzO7w8yOyXeAu79QRky55KvQjWRG\nJ04zs55m1oPMTPB/AW4L44sbGho23e6mmf5iDFRbLvr0yT3OqG/f8p+Dr7ZcREm5yNxdNjQ0RHLu\nUpoBW4G7gWnu/mIkUWW+5wbgWGAPoD+wmMwd1uHu3tLuuIHALcDw7P4lwER3fzOEGNQMKFXvyiuv\n5D/+4z86bBsyZAjTp0/XIAuJRFIWX1zp7oPCDCKpVKwCao8PVFMu3J1hw4bx/PPPs99++zFw4ED6\n9u3L+PHjQylU1ZSLqCkXgaSMBnzRzHbJtaRGGzO7093HlhFXYjQ0NFBfX6+LUKrSY489xvPPP8+g\nQYN49tln2XJLLS4g0Wlqaoqs26SUO6v9gf8AvuvuOZcCCXM0YJx0ZyXVzN0ZMWIEzzzzDDfddBOX\nXnpp3CFJSiSlGbD9Q8HrgHfZ/KHg3dy9V1hBxkXFSqrZU089xVFHHcWOO+5Ic3MzW2+9ddwhSUok\nZej6LsCLwJPA/5KZh++Ndq/XgU9Dik8SQiMiA9WSixEjRjB79mxuvPHGyApVteSiEpSLaJXSZ/WB\nux/d1QFm9naJ8SSO+qykWvXs2ZPTTz897jAkRZLWZ3VMd7OUm9khETxrVXFqBhQRKV4i+qzSRMVK\nRKR4SemzwjLGmtkcM3sl+/qlmZ0bZnCSHGqPDygXAeUioFxEq+hiZWa9gUeBWcDxwGezrzHAnWb2\nmJlV/UhAkWr00ksvceONN7JmTa5l3kSqVyl9VtcA/wrcCDwMtA2m2AU4EZgE/Mjd/z3EOGNhZj5l\nyhQNsJCRoAM+AAAVE0lEQVSq8fWvf51f/OIXXHLJJXz/+9+POxxJmbYBFlOnTo2/z8rMlgJnuftz\nefYPA+5x96EhxBcr9VlJNVm6dCn77bcfPXv25LXXXmOPPfaIOyRJqaT0WW2br1ABuPuz5F4QUaqY\n2uMDSc3F9ddfT2trK2PHjq1YoUpqLuKgXESrlGLVkl0zKiczGwR8UnpIIlKs119/nbvuuosePXow\nefLkuMMRCV0pxWoucJ+ZHdZ5R7YJ8D7gsXIDk2RRn10gibm49957+eSTTzj77LP5zGc+U7HvTWIu\n4qJcRKuUGSyuBn4P/NHMVtJxgMXOZKZcOi2c8OKnGSykGkyaNIkvfOELDB5c9fNHSxVL1AwWsGnB\nw+vIFKVts5s/Au4FrnT3d0OLMEYaYBHQWj0B5SKgXASUi0BS1rPC3VcBF5rZN4GdspvfdfdWADPb\nz93/FFKMIiKScpFMt2Rmb7n7LqGfuMJ0ZyUiUrzE3FmZ2UnA18j0U3WercLIrHUlIiISilKmW5oI\nPACcAxwE7N3ptVcp55Vk0zMkgaTk4qabbuKyyy5j5cqVscWQlFwkgXIRrVLurL4FXAnc7O4bcx2g\n9axEorVmzRquu+46PvzwQ/7xH/+RQYMGxR2SSLJGA5rZu8DArjpzzKze3ZvKjC126rOSpLrhhhuY\nPHkyX/rSl3jqqacwC7V7QKQsSZlu6SVghwjOKyIF+PjjjzdNUnvVVVepUEkqlFJUvgX8yMwO7uKY\nu0uMRxJK7fGBuHPxs5/9jFWrVjFs2DBGjx4dayxx5yJJlItoldJn9Utga+B0M/sYeA9obbffgAEh\nxCYiObz22muA7qokXUrps9oAPE2mKOXzRXfvW05gSaA+K0mqP//5zwwdOpQePdTiLskTRZ9VKcXq\nbXfvcgKyQo6pBipWIiLFS8oAi3MLOOYrJZxXEkzt8QHlIqBcBJSLaBVdrNx9fvv3ZvbFHMe8UE5Q\nIiIi7ZU9N2CtNPnlYmY+ZcoUPRQssWtpaaFXr84zm4kkS9tDwVOnTo2/z2qzE9R4sVKflcTN3Rk+\nfDgHHHAAN998MwMGaLCtJFtS+qwkhdQeH6h0Lh555BGee+45fv3rX7P11ltX9Lu7o+sioFxES8VK\nJMHcnWuvvRaAyy67jL59q/6JEJGShNEMOMLdnw4pnkRRM6DE7fHHH+e4445jp512Yvny5fTr1y/u\nkES6lchmwFotVCJxmjNnDqNHj+a0004D4IQTTlChklRTM6AURO3xgahzMWfOHC6++GLmzZvHRx99\nBMCTTz7JnDlzIv3eUui6CCgX0VKxEkmYGTNmsGzZsg7bli9fzsyZM2OKSCR+ZfdZbTqR2deBLwF/\nAm5393Vm9hngGGCVuz8QyhdVkPqsJA719fU88cQTm20fOXKk/nqXqhBFn1Ups65vxsymAP8H+ANw\nGDDezEa7+2tmtg74K7qLEylInz59cm7XSEBJs7AKyH7APu5+ursfAZwJTDez3YFPQ/oOiZH+og9E\nnYsJEyYwZMiQDtuGDBnC+PHjI/3eUui6CCgX0Qrlzgr4vbuvb3vj7ovM7EzgO8CjIX2HSCqMGDGC\n6dOnM3PmTNavX0/fvn0ZP348Y8aMiTs0kdiE0mdlZqcA/YGpwPHu/mK7ff8M/NDdwyqMFaO5AaXS\nNmzYQF1dHQceeCCzZ89m++23jzskkYIlem7ATScy2xs4CHjU3Vs67fuyu/8mlC+qIA2wkEq78847\nGTduHAcffDCLFi3SSsBSlRL5UDCAmR1Lpt/ql50LFUA1FirpSO3xgahy4e7ceuutAEycOLEqCpWu\ni4ByEa1QipW7P5491y/N7AYz2y+M84qkyZNPPsmiRYsYOHAgZ599dtzhiCRKwc2AZvY5d3++m2N6\nAbOBr1ZjH1VnagaUSjrppJN46KGHmDJlCg0NDXGHI1KyuJsBL+vugGwT4IXAxpIjEkkhd+ewww5j\nt91246KLLoo7HJHEKaZY7VzIQe7+PqDJbWuM2uMDUeTCzLjmmmtobm5m550L+qeWCLouAspFtIop\nVkeb2UIz+46ZjTCznl0c+065gYmkUc+eXf2zEkmvYvqs/ga8Ceyb3fR34LfAQmAB8GxbB4+Z3eXu\n54YfbmWpz0pEpHhxzw34K3c/w8wGAaPavUZn9//NzJ4kU7z2DDPIsJlZX+B7ZO4sewGHuPuR8UYl\nIiL5FNMMOA3A3Ve6+8/d/UJ33xvYC7gAeAQYBtxMZvb1JLseWODu/9fdvwVcHndASaf2+ECYuXj/\n/fep5rt3XRcB5SJaBRcrd382z/bX3X2Wu5/r7rsB+wOvhhVg2MxsS+B8YDsz+56Z3QasiTksSamT\nTz6ZQw45hJdeeinuUEQSLbTpljqc1Owedz8rhPMMBmYBx7l7WLNt7Ae8BFzj7tdm3z8B7OvuH3Q6\nVn1WEplnn32W4cOHs+2227JixQq22WabuEMSCUXcz1kV45pyT5CdHPdpYG8gb8Uws4FmdreZvZJ9\n3Wtmu7bbP8HMXs2+RgFtvxEeBHD3P5EZOHJsuTGLFGP69OkAXHjhhSpUIt2IpFi5+9IQTnM5mVWG\nnwZyVmgz6w08TmagyP7Z11pgoZn1y8Yyw92HZl8LgBXZj7e2O9VGoHcIMdcstccHwsjF22+/zT33\n3EOPHj0SuU5VoXRdBJSLaCV59d4R7r6sm2PGkpnpfbK7t7p7KzCZzN1YzmkA3P0tMkPujwIws+2B\nIWSaAkUq4kc/+hEtLS2cfPLJ1NXVxR2OSOJF0mcVJjNrBM7L1WdlZnPJrFC8V6fti4G12VWLc51z\nD2AGsBzYDfgvd/9ljuPUZyWRuOuuu7j22mu54447+NKXkj54VqQ4UfRZVXuxegt4xd1Hddr+MDDK\n3bcu87tVrCQyra2tmFlVLAUiUoxqGmBRKQPIPez8I2ArM+tT4XhqltrjA2HlokePHlVfqHRdBJSL\naFV7sdJtj4hIClT7mlPvEQxFb29bMn1WG8r9gnHjxm3qAO/fvz+HHnoo9fX1QPCXVBre19fXJyoe\nvU/O+zZJiSeu923bkhJPJd83NTXR2NgIENmAoWrvs3qMzMO8nQdYLAHWuPuIMr9bfVYSmtWrV9Ov\nXz969eoVdygikUpzn1W+ivEAsKeZbZo418x2JjMz/P2VCCwtOv8VnWal5uLKK6+krq6Oxx57LNyA\nYqTrIqBcRKtailW+Ct0ILAGmmVlPM+tBZpLavwC3hfHFDQ0NugilbB988AF33nknb731FnvumehF\nCURK1tTURENDQyTnTmwzoJndQGYKpD2A/sBiMndYh7t7S7vjBgK3AMOz+5cAE939zRBiUDOghGLa\ntGlcccUVHHfccfzqV7+KOxyRSKXyOas4qVhJGFpaWth7771ZsWIFjz76KMcff3zcIYlEKs19VrFR\nM2CGchAoNhcPPPAAK1asYN9992X06NHdf6CK6LoIKBfRNgNW+9D1yEWVeEmPgQMH8uUvf5lzzjmH\nHj3096HUrvrsYy5Tp04N/dxqBuyCmgElTK2trSpWkgpqBhSpYipUIqXTv55uqM8qQzkIKBcB5SKg\nXKjPKlbqsxIRKYz6rGKiPisp1bp169i4cSPbbbdd3KGIVJz6rESqxKxZs9h999354Q9/GHcoIjVB\nxUoKovb4QHe5aG1tZfr06axZs4aBAwdWJqiY6LoIKBfRUrHqhgZYSLEee+wxli5dyh577MHJJ58c\ndzgiFZPKuQGTQH1WUopjjz2W+fPnc+ONNzJp0qS4wxGpOM0NWGEqVlKsJUuWcPDBB9OvXz/++te/\nsv3228cdkkjFaYCFxEZNoYGucrHFFltwyimncMEFF6SiUOm6CCgX0dJzViIhmDNnDjNmzGDDhg30\n6dOHY489Nu6QRGqKmgG7YGY+ZcqUTQ+6ieQyZ84cLr74YpYtW7Zp25AhQ5g+fTpjxoyJMTKRympq\naqKpqYmpU6eqz6qS1GclhRg9ejTz5s3LuX3u3LkxRCQSL/VZSWzUHh/onIsNGzbkPG79+vUViCZe\nui4CykW0VKxEytSnT5+c2/v27VvhSERql5oBu6BmQCmE+qxEOoqiGVCjAUXK1FaQZs6cyfr16+nb\nty/jx49XoRIJkZoBu6HpljKUg0CuXIwZM4a5c+fS1NTE3LlzU1OodF0ElAutZxUrrWclIlIYrWcV\nE/VZiYgUT0PXRUQklVSspCBqjw8oFwHlIqBcREvFSkREEk99Vl1Qn5WISPHUZyUiIqmkYiUFUXt8\nQLkIKBcB5SJaKlbd0EPBIiKFifKhYPVZdUF9ViIixVOflYiIpJKKlRRETaEB5SKgXASUi2ipWImI\nSOKpz6oL6rMSESme+qxERCSVVKykIGqPDygXAeUioFxES8VKREQST31WXVCflYhI8dRnJSIiqaRi\nJQVRe3xAuQgoFwHlIloqVt3Q3IAiIoXR3IAxUZ+ViEjx1GclIiKppGIlBVFTaEC5CCgXAeUiWipW\nIiKSeOqz6oL6rEREiqc+KxERSSUVKymI2uMDykVAuQgoF9FSsRIRkcRLZZ+Vma0DVrfbtC1wnrvf\n3+k49VmJiBRJfVbh+Ym7D257AYuAOXEHJSIiuaWyWLn7xW3/b2ZfAZ5y9/UxhpR4ao8PKBcB5SKg\nXEQr0cXKzAab2Vwza43way4Cfhjh+UVEpEyJ7bMys1OA7wMtwBB375nnuIHALcCw7KYlwER3fzO7\nfwIwPrvvn919QbvPDgW+5+5n5Dm3+qxERIqUtj6ry4FjgKeBnD+0mfUGHge2APbPvtYCC82sH4C7\nz3D3odnXgk6nGA/MjCh+EREJSZKL1Qh3X9bNMWOBg4DJ7t7q7q3AZGBvMs17eZnZNsCh7v5UKNHW\nOLXHB5SLgHIRUC6ildhilS083TkVeN3dm9t97h3g5ey+rowFZpUcYMosWrQo7hASQ7kIKBcB5SJa\nW8QdQJkOBl7Jsb0ZGNXVB939B1EEVKtWr17d/UEpoVwElIuAchGtxN5ZFWgAsCbH9o+ArcysT4Xj\n6VYxTQWFHJvvmEK3d/U+6mYN5SL/d5d7rHLR/TG5theyrVK5KPbctZwLqP5iVXVD9ar1l1Jzc3O3\nsRRLucj/3eUeq1x0f0zYv6DDzoWKVUeJHbrexswayUyFtFlhNbM3gT+7+6hO2x8Gjnb3bcr87mQn\nR0QkocIeul7tfVaLgX1zbN+LzPNWZQk72SIiUppqaQbMd4fzALCnme3ZtsHMdiZTwO7P8xkREaky\n1VKs8t3hNJK5g5pmZj3NrAdwPfAX4LYKxSYiIhFLbLEysxvM7Hngq4Cb2fNm9pyZ9Wo7xt1bgGOB\nT8k8W/UysDUwyt0/rmCst5rZ7WZ2k5nNN7NzKvXdSWJmO5jZHWY2w8ymm9kvs1NapZaZ9TKzyWb2\ndzPbI+54Ks3MdjOzB83sNjN7xMyGdf+p2pT2a6FNqb8nEj/AohqY2Q3ufnn2/78APAH0K/DB5pph\nZocA33D3b2Xffws4092PjDey+GRz8L/A74A6d38j5pAqKjvY6R53/7mZHQw8CAxN278N0LXQptTf\nEypWITOz84AJ7j487ljiZmZjgJnuvnfcscQtu3JAqn5BmdkOwLvA9u7+UXZbMzDW3Z+IM7Y4pfFa\n6EqhvycS2wwYhSiXHDGzg81sNjABOD3s84etQsuvfBWoiplCKpSPqhFSPuqAjW2FKmtldnvV0LUR\niCgXBf2eSE2xyi458jSZSW7z3k6a2UAzu9vMXsm+7jWzXdvtn2Bmr2Zfm57vcvfF2aVG/g14MvtX\nZSJFnYvsvjHA1u5+c1Q/R1gqkY9qElY+aoFyEYgiF8X8nkhNsSKiJUfMrEfbvuz+3wPrgKMi/WnK\nE+nyK9kL8KvAuZH+FOGpxHI01SSUfACvA73NbNt2Hx1EZu7OahFWLmpBqLko+veEu6fiBfTI/rcR\naM1zzDeAtvbktm07A58Ak/J8pg74Rbv3O5KZm/CguH/mSucie8zpwK3t3k+P++eNMx/tjm0F9oz7\nZ610PoCHgXOy/38IsIxsX3k1vKK4NqrpWojwuij690Rq7qw8uiVH3gfMzGaZ2feBu4CL3b3sGTSi\nElUusqO9fg6cYWZvm9nbwIXlRxytCK8NzOxIM5tJptmkwcy+Vma4kQs5H/8KnG5mtwHXAWd59rdT\nNQgzF9V4LbQXVi5K/T1R7dMtha3oJUfcfQ1wVoQxxaWUXCwGeuXaVwNKWo7GM4t7PkVmVepaUlA+\n3H0FcFKFYopLobmo1WuhvW5zUervidTcWRWo6pYciZBy0ZHy0ZHyEVAuApHlQsWqo6ppnqgA5aIj\n5aMj5SOgXAQiy4WKVUfvAbmWFdkWWOvuGyocT5yUi46Uj46Uj4ByEYgsFypWHS0ms7xIZ6EsOVJl\nlIuOlI+OlI+AchGILBdpLVZaciSgXHSkfHSkfASUi0DFc5HWYqUlRwLKRUfKR0fKR0C5CFQ8F6kp\nVlZFS45ETbnoSPnoSPkIKBeBuHOhWddFRCTxUnNnJSIi1UvFSkREEk/FSkREEk/FSkREEk/FSkRE\nEk/FSkREEk/FSkREEk/FSkREEk/FSkREEk/FSkREEk/FSkREEk/FSqQKmNlQM3vXzOoKOHY/M/uH\n6KMSqRwVK5HqcBWwI/DdAo69ATg52nBEKkvFSiThzGwY8BrwMPB1Mzuoi2N7AF8CFlQoPJGK0BIh\nIglnZv8F/AtQB7wA/Mrd/zHPsZ8D/gjs7O7vVSxIkYjpzkokwcxsFPA7d//Y3V8G/gs4wcy+3Om4\nU8zsLuBO4F3gFjP7YeUjFomG7qxEEszM7gbOc/dPs+93A5YCz7n7l3Mc/z/AX919fGUjFYmW7qxE\nEsrMTgMeaitUAO6+AvgBMMLMTux0fE9gJPBERQMVqQDdWYkkULbw3Onu/5RjX3/gL8CbwMGe/Uds\nZl8AngEGufuqSsYrEjXdWYkk0/nAHbl2uPtqYBpwAHBuu12jgD+rUEktUrESSRgz6wsc4e5dDT+f\nDrwFTDWz3tlto8g2AZpZbzP7TrSRilTOFnEHICKbGQ+0mtnEbo57CTgWuIhM8RpE5lksgIuB/xdZ\nhCIVpj4rkQQxs37A68AOBX7EgfeAPcnMWnE28Bwwz92fjiRIkRioWImISOKpz0pERBJPxUpERBJP\nxUpERBJPxUpERBJPxUpERBJPxUpERBJPxUpERBJPxUpERBJPxUpERBJPxUpERBLv/wOzbloQE78r\nKAAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f20760b4b90>"
       ]
      }
     ],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "error"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 18,
       "text": [
        "array([  3.81125927e-05,   9.41813943e-06,   2.25089054e-06,\n",
        "         4.63970974e-07])"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See how the error drops four times when the time step is halved? This method is second order in time!\n",
      "\n",
      "Clearly, Crank-Nicolson (circles) converges faster than backward Euler (squares)! Not only that, but also the error curve is shifted down: Crank-Nicolson is more accurate.\n",
      "\n",
      "If you look closely, you'll realize that the error in Crank-Nicolson decays about twice as fast than backward Euler: it's a second versus first order method!"
     ]
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "Spatial convergence"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "To study spatial convergence, we will run the code for meshes with 21, 41, 81 and 161 points, and compare them at the same non-dimensional time, say $t=20$. \n",
      "\n",
      "Let's start by defining a function that will do everything for us"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx_values = numpy.asarray([11, 21, 41, 81, 161])\n",
      "\n",
      "dt = 0.1\n",
      "error = numpy.zeros(len(nx_values))\n",
      "\n",
      "t_final = 20.0 \n",
      "\n",
      "x = numpy.linspace(0,L,nx)\n",
      "\n",
      "for i,nx in enumerate(nx_values):\n",
      "    \n",
      "    dx = L/(nx-1)\n",
      "    x = numpy.linspace(0,L,nx)\n",
      "   \n",
      "    sigma = alpha*dt/(dx*dx)\n",
      "\n",
      "    nt = int(t_final/dt)\n",
      "    \n",
      "    A = generateMatrix(nx, sigma)\n",
      "\n",
      "    T = numpy.zeros(nx)\n",
      "    T[0] = 100\n",
      "      \n",
      "    T = CrankNicolson(T, A, nt, sigma)\n",
      "    \n",
      "    T_exact = T_analytical(x, t_final, 100, alpha, L)\n",
      "    \n",
      "    error[i] = L2_error(T,T_exact)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And plot!"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plt.figure(figsize=(6,6))\n",
      "plt.grid(True)\n",
      "plt.xlabel(r'$n_x$', fontsize=18)\n",
      "plt.ylabel(r'$L_2$-norm of the error', fontsize=18)\n",
      "plt.xlim(1e-4,1)\n",
      "plt.ylim(1e-4,1)\n",
      "plt.axis('equal')\n",
      "plt.loglog(nx_values, error, color='k', ls='--', lw=2, marker='o');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAGXCAYAAAANqX+7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmY1OWV9//3AVlVRgm2LC4sgugTV8wYGZcGNa32iImJ\nWTQoLshoAmYyRp088aFJ5sloxgkiueKARiGj6KjRazStxLUYHZP8YtwdcQGBCBhEnMcF2ezz+6Oq\nmaaprq7lrvreVfV5XVddsb71rW+dPin69L1879vcHRERkZj1SDoAERGR7qhYiYhI9FSsREQkeipW\nIiISPRUrERGJnoqViIhEb5ekA6g0MxsIXAd8BDgwEviuu7+RaGAiItKluitWwL7ARnefAWBm3wZu\nAY5LNCoREelS3RUrd38B+HaHQ28BwxIKR0RE8hD1mJWZDTGzxWbWVsaPOR34WRmvLyIiJYq2WJnZ\nmcDTpMeUulwTyswazOx2M1uaedxtZsM6vD7DzN7IPCZ2em8zsJu7/7RcP4eIiJQu2mIFXAGcRLpg\nWbYTzKw38Ajp7syDM4+PgSfMbFcAd7/B3UdnHo93eG8z6VbV5LL+FCIiUrKYi9V4d1/WzTnnAYcA\nV7p7m7u3AVeSbo1d0tWbzOws4GR3/xt3dzObEyxqEREJzmJfdd3MFgDnuvtOhdXMFgMHuvuITsdf\nBD5292OyvOdQ4I/Au/xPi22Au+8aOnYREQmj2mcDHgoszXJ8BTAxy3Hc/UWgVxljEhGRwGLuBszH\nIODDLMc/APqbWZ8KxyMiImVQ7S2rsvZhmlncfaQiIpFy96wT44pV7S2r9cDuWY4PID1mtbnUD3D3\nsj9mzpxZ9vd1d26u17O9ls+x7p7HlMtK5bOQ4/WSz9DfzXxzpXzm91oxuSuHai9WLwIjshwfAbxU\n4ViK1tjYWPb3dXdurtezvZbPsc7PV6xYkTOGEIrNZaHvLTafhRyvl3yG/m5mO55PfpXPcP/Wy6Fa\nZgNOdveeWV6bCswDRrj7ysyxvYG3gavc/Z9L/GyPPT/VZMqUKSxYsCDpMGqG8hmW8hmOmeF12g3Y\n1Q+9gHQL6loz62lmPYBrgOXAjSE+uKWlhVQqFeJSdW/KlClJh1BTlM+wlM/SpVIpWlpaynLtaFtW\nZvYT4GRgP2AP0l1+Dhzt7ls7nNcAzAaOyrz+EvAdd18dIAa1rEREClSOllW0xSoGKlZhpVKpivRt\n1wvlMyzlM5x67gYUEZE6ppZVDmbmM2fOpLGxUX9xiYh0I5VKkUqlmDVrlroBK0ndgCIihVM3oFQ1\nzaoMS/kMS/mMm4qViIhET92AOagbUESkcOoGTIBuChYRyU9d3hQcA7WswtJ9LGEpn2Epn+GoZSUi\nInVJLasc1LISESmcWlYiIlKXVKy6oQkW4SiPYSmfYSmfpSvnBItq39a+7MqVeBGRWtO+NN2sWbOC\nX1tjVjlozEpEpHAasxIRkbqkYiUVozGBsJTPsJTPuKlYiYhI9DRmlYPGrERECqcxqwRo6rqISH60\nNmBC1LIKS2uvhaV8hqV8hqOWlYiI1CW1rHJQy0pEpHBqWYmISF1SsZKK0USVsJTPsJTPuKlYScW1\ntrbS1NREY2MjTU1NtLa2Jh2SiEROY1Y5aMwqvNbWVi677DKWLVu2/dioUaOYM2cOzc3NCUYmIqFo\nzCoBus8qrBtuuGGHQgWwbNky5s6dm1BEIhJKOe+zUrHqRktLi+69CCSVSrF58+asr23atKnC0VQ/\n/REVlvJZusbGRhUrqQ19+vTJerxv374VjkREqok2X5SKaWxs5OOPP2bZsmU7jVlNnz49wciqk1r8\nYSmfcdMEixw0waI8Wltbueaaa3jqqafo0aMHixYt4mtf+1rSYYlIIJpgIVWtfUygubmZJ598kokT\nJ9LW1sa7776bbGBVSmMsYSmfcVOxksRMmzaN3r1788477yQdiohETt2AOagbsLy2bNnCBx98wKBB\ng5IORUQCKkc3oIpVDipWIiKF05iVVDWNCYSlfIalfMZNxUpERKKnbsAczMxnzpxJY2Oj7sEQEelG\nKpUilUoxa9YsjVlVksasKueJJ57g+uuv5/LLL+e4445LOhwRKYHGrKSq5RoTeOyxx7j//vv5l3/5\nl8oFVOU0xhKW8hk3FSuJwtSpUzEz7rnnHtavX590OCISGXUD5qBuwMpqbm7mwQcf5J/+6Z+4/PLL\nkw5HRIqkbkCpadOmTQNg/vz56I8EEelIxUoqprsxgdNOO41hw4axbt26nTZolJ1pjCUs5TNu2iJE\norHLLrvw61//mjFjxtC/f/+kwxGRiGjMKgeNWYmIFE5jViIiUpdUrKRiNCYQlvIZlvIZNxUrERGJ\nnsasctCYVXI2bdrEXXfdxbJly5g1a1bS4YhIAbSfVYWpWCVnzZo17LfffgD86U9/YsiQIQlHJCL5\n0gSLBLS0tKgvO5BC8jh06FAmTZrEp59+yi233FK+oKqYvpdhKZ+lS6VStLS0lOXaKlbdaGlp0fYg\nCWlf0eKmm27i008/TTgaEelOY2Nj2YqVugFzUDdgstra2jjggAN46623ePDBBzn11FOTDklE8qBu\nQKkrPXr04OKLLwbg7rvvTjgaEUmSipVUTDFjAhdccAH3338/8+fPDx9QldMYS1jKZ9y0NqBEraGh\ngdNPPz3pMEQkYRqzykFjViIihdOYlYiI1CUVK6kYjQmEpXyGpXzGTcVKqsaGDRuYPXs2q1atSjoU\nEakwjVnloDGruEyePJnbbruNq6++mh/+8IdJhyMiXdDagBWmYhWXJUuW0NjYyJAhQ1i5ciW9evVK\nOiQRyUITLKSqlTomcPzxxzN27FjWrl3Lr3/96zBBVTGNsYSlfMZNxUqqhpltXy9w3rx5CUcjIpWk\nbsAc1A0Ynw0bNjB06FC2bNnCihUrtm8jIiLxKEc3YF2uYGFm1wO7A+8DhwO3uvvtyUYl+Rg4cCBz\n587ls5/9LPvuu2/S4YhIhdRrN+AWd7/Q3S8Hvg/cbGb1mouKCTUmMHXqVI455hjMgv7hVnU0xhKW\n8hm3gltWZnYr4MAid380fEjl5+5XdHg6FnjF3duSikdERHIreMzKzNqA24Fr3f3lskT1P581BLgV\n+IK7B235mNmhwA+AkcBZ7v5WlnM0ZiUiUqBYpq6vc/fJFShUZwJPky4mXVYMM2sws9vNbGnmcbeZ\nDevw+gwzeyPzmNh+3N1fdPevAt8C/sPMBpbxxxERkRIUU6xeNrOhuU4ws4VFxtPRFcBJpAtW1gpt\nZr2BR0h3Zx6ceXwMPGFmuwK4+w3uPjrzeNzMerS/lnn998AnwPEBYpYcyjEmsGLFChYvXhz8utVA\nYyxhKZ9xK6ZYzQBuNLOjcpzzhSLj6Wi8uy/r5pzzgEOAK929LTPudCXp1tglXbxnP+Dm9idm9hlg\nMNDdZ0lkli5dysiRI5k8eTKbN29OOhwRKaNixqzeAnYDPkO6RfIu0HFyggH7uHuQtXDMbAFwbrYx\nKzNbDBzo7iM6HX8R+Njdj8nynt2BmzKxbwAOAu5291uznKsxq4i5O0cccQQvvPACixYt4hvf+EbS\nIYkI8dxnNZQcXXMZQ4oLp2CHAkuzHF8BTMxyHHf/EPh6GWOSCmlf0eLSSy9l3rx5KlYiNayYbsAN\n7j7B3Ru7epC+2bYSBgEfZjn+AdDfzPpUKA7JQznGBM455xx23XVXlixZwtKl2f5uqV0aYwlL+Yxb\nMcVqch7nnFLEdYuhPro6N2DAAM4++2wA5s+fn3A0IlIuBXcD5nMjsLu/UFw4BVtPetmkzgaQHrMq\nedR9ypQpDB8+HIA99tiDww8/nMbGRuB//hLT8/yetx8Lff1LL72U3XffncMOO6ws14/1efuxWOKp\n9uftx2KJp5qep1IpFixYALD992VoRS1ka+l1bs4FvgqMyhx+A7jL3f81XHjdTrB4CBibZYLFS8CH\n7j6+xM/WBAsRkQJFcVNw5t6mB0mvLHEqMCbzaAYWmtlDZhZ6V7yuKsa9wP5mtn+H+PYmvYTSrwLH\nICVq/0tMwlA+w1I+41bMmNVVwBHA94ADSXe5DSBdIK7IvPb3oQLM6KpCLwBeAq41s56ZxWivAZYD\nN4b44JaWFn2JRUTykEqlaGlpKcu1i7nP6nXg6+7+bBevjwPudPfRJQVm9hPgZNI38e4BvEi6hXW0\nu2/tcF4DMBs4KvP6S8B33H11KZ+fuba6AUVEClSObsBiitU77j641HOqgYpVdXr22WcZO3Ys/fv3\nTzoUkboUxZgVsDUzLpSVmQ0GthUfktSqSnSnNjU1MW7cOI488kiamppobW0t+2cmRd3TYSmfcStm\nBYvFwD1mdlnnrsBMF+Ac4KEQwYkUorW1leeeew6A1157jddee41ly9JLPjY3NycZmoiUqJhuwMHA\n74F9gXeAtZmXhgJ7A6uAz7v7OwHjTISZ+cyZM2lsbNzhXgyJU1NTEw8//HDW4/W6MrtIJaVSKVKp\nFLNmzUp+zAq2T2r4MfAV0jMBIb3E0d3A99393WARJkhjVtWlsbGRJUuW7HT8hBNOUBePSAXFMmaF\nu69z94uAgaQXrR0CDHT3qbVSqCS8cheMPn2yLwXZt2/fsn5uUlSAw1I+41bMTcG3mtktZnZSZg+p\nP2cebd2/W6R8ZsyYwahRo3Y4NmLECKZPn55QRCISSjFjVm3A7cC15d7aPmnqBqw+ra2tzJ07l02b\nNtG3b1+mT5+uyRUiFVY191nVCk2wEBHJX1QTLMzsUdILy67Jcc5Cdz+v1OCSppZVWB1XtJbSKZ9h\nKZ/hxDLBYgZwo5kdleOcLxQZj4iIyE6KaVm9BewGfAb4BHgX6Di5woB93D30yusVp5ZVbWhra2Pd\nunUMHlwXvdciiStHy6qYFSyGAf9J1yuhQ3oqu0jiXnnlFSZNmsTAgQP5wx/+kHQ4IlKkYroB33P3\nCe7e2NUDeD9wnInRFiHhJJHHkSNH8v777/PMM8/w7LNZNwqoWvpehqV8lq6cW4QUU6y2td9nleOc\nU4oNKDYtLS0adK1i/fr149xzzwVg/vz5CUcjUtsaGxuj2s9K91lJVXn11Vc5+OCD2W233VizZg27\n77570iGJ1LRYZgOuc/fJtV6opHYcdNBBHHfccXz00UfccccdSYcjIkUopli9bGZDc51gZguLjEdq\nWJJjAtOmTWPcuHHsvXeXW7FVHY2xhKV8xk33WUldOPvss3nmmWc444wzkg5FRIqg+6xy0JiViEjh\nYrnPaijwNHVyn1X7bEDNCBQRya19bcByKKZltdbdcxajfM6pBmpZhaW118JSPsNSPsOJZTbg5DzO\nqZn7rKQ2bdmyBf0hIlI9itrWvl6oZVWbfvzjHzN79mwefPBBPve5zyUdjkjNiaVlhaV9yczmmdmv\nMsdGZ471DBmgSGjvvfce69evZ968eUmHIiJ5KmZb+37AY8CvgKnAyZmXdgNuBh41My0RIDuJ5T6W\niy++GIA77riDDz74IOFoihdLPmuF8hm3YlpWVwMHAt8CjgY2Arj7c8A+wHvA90MFKBLagQceSGNj\nIxs3buT2229POhwRyUMxswHfAL7p7r/PPN9h5p+ZDQKedvcxQSNNgMasatedd97JN77xDQ477DCe\ne+45zIJ2r4vUtVjGrAa2F6ps3H096S7BmqAtQmrTl770JRoaGth777358MMPkw5HpCaUc4uQou6z\nAo5097Xtzzu1rEYBS9x9n6CRJkAtq7Biu4/lgw8+YMCAAUmHUbTY8lntlM9wYmlZPQjcZWYHdH7B\nzI4F7gEeKDUwkXKr5kIlUm+KaVkNAX4H7Au8RXpSxSukt7vfC1gBfN7d1wWNNAFqWYmIFC6KllWm\n++9zwC2kF7PtBRwO9AZuAv6yFgqViIjEo6ibgt19nbtfBAwkvWjtEOAz7j4tM8FCZCeaqBKW8hmW\n8hm3oopVO3dvc/c/Zx5t3b9DJD7Lly/n/PPP51vf+lbSoYhIF7Q2YA4as6oPb775JqNHj6Zfv36s\nWbOGPfbYI+mQRKpaFGNWIrXmgAMO4MQTT+STTz7htttuSzocEclCxUoqJuYxgWnTpgEwb968qtk6\nJOZ8ViPlM24qVt3QChb14YwzzqChoYGXX36Z3/72t0mHI1KVolrBop5ozKq+XHXVVcyePZs5c+bw\nN3/zN0mHI1K1yjFmVVSxsvSqn18kvSPwIHf/spmNBj4L3O/un4YMMikqVvVl/fr0XReDBg1KOBKR\n6hbFBAvtZyXFir07ddCgQVVVqGLPZ7VRPuOm/axERCR62s8qB3UDiogULopuQOpsPysREUleMcVq\nS2bl9awy+1lp6SXZSTWNCTz00EOceeaZbNiwIelQulRN+awGymfctJ+VSBbXX3899913H7/85S+T\nDkVE0H5WOWnMqn7de++9fPnLX+aggw7ilVdeIX23hojkI4oxK+1nJfXg9NNPZ/Dgwbz66qs89dRT\nSYcjUve0n5VUTDWNCfTq1YsLL7wQSK8XGKNqymc1UD7jVpb9rMzsoNJDE0nWRRddhJnx2GOPsXnz\n5qTDEalrZVkb0MzWuPvQ4BeuMI1ZycMPP8yxxx5L//79kw5FpGrEtDbgF4EzgKGkx6x2eBkY7+59\nSg8vWSpWIiKFi2KChZl9B7gXOAc4BBjZ6TGimOvGSluEhKM8hqV8hqV8li6qLULM7E1gPnC9u2/p\n4pwdlmCqVmpZhZVKpWhsbEw6jJqhfIalfIYTRTegmb0LNOT6LW5mje6eKjG2xKlYiYgULopuQNI3\nAA8sw3VFovXJJ5+wcOFC7rrrrqRDEalLxRSVbwM/N7PDcpxze5HxSA2r5jGBRx55hClTpnD11VcT\nS2u7mvMZI+Uzbrt0d4KZvQV0/te5O3CWmW0E1rPjwrUGVM8OdiJ5OO200xg2bBivv/46qVSKCRMm\nJB2SSF3pdszKzDYDT5MuQvn6vLv3LSWwGGjMSjqaOXMmP/zhD/na177GnXfemXQ4ItFKZIJFMTP7\nNBtQatGqVasYMWIEPXv25O2336ahoSHpkESilNQEi8n5XszM2m/zbyouHKll1T4msN9++3HkkUey\ndetWxo8fT1NTE62trYnFU+35jI3yGbdux6zc/dGOz83s5+5+aRen/8TMvg6cB7wYID6RaLS2trJ2\n7VoAli1btv0B0NzcnGRoIjWvmPusuuziM7MG0q2q77t71S9mq25A6aipqYmHH3446/HFixcnEJFI\nnGK5z6pLmX2sbiO9z5VITelq5fVNmzZVOBKR+tNtNyCAmd1Kevq6AX9hZrfkuN5YYE2Y8KSWVPty\nNn36ZF+buW/fZCa+Vns+Y6N8xi2vYkV6DKqjKV2ctxF4Fbik2IBEYjVjxowdxqkARo0axfTp0xOM\nSqQ+BB2zqjUas5LOWltbmTt3Lps2baJHjx4cdNBBzJ07lx49tMKYSLtYFrI9391vDRlEUszscuAn\n7p71N42KlXTF3RkzZgxvvvkmjz76KCeeeGLSIYlEI4oJFjVUqD4LNLLzUlJSJrV0H4uZ8c1vfhOA\n+fPnJxJDLeUzBspn3Oqy78LMegE/Av6ewpaREtnuggsuoEePHtx3332sW7cu6XBEalrUxcrMhpjZ\nYjNr6/7sgswE5gAfBr6u5FBrM6323XdfTjvtNLZu3cqCBQsq/vm1ls+kKZ9xi7ZYmdmZpBfQHUmO\nrjozazCz281saeZxt5kN6/D6DDN7I/M40czGA/1rYXNISd60adOAdFdgW1vov6lEpF20xQq4AjiJ\nHCu+m1lv4BHSU/APzjw+Bp4ws10B3P0Gdx+deTwGTAL2NLMbgX/IXOfnZvaVcv9A9a4WxwROOeUU\nLrvsMm677TbMKtujXIv5TJLyGbd89rP6ATAROMvd3yt/SNuNd/e2bn4BnAccApzh7m0AZnYlsJr0\nvV7XdX6Du1/V/t9mNhw4O8dahyI57bLLLlx//fVJhyFS8/LZImQZ8APg3zLF4zx3X1iR6NKfvwA4\nN9v0cjNbDBzo7iM6HX8R+Njdj8lx3ROA80mvKv9zYL67v9TpHE1dFxEpUFL7Wb3j7oM7PO/2pmAz\ne97dDw8SYO5itQZY6u4TOx2/H5jo7ruV+NkqViIiBUrqPquNZtZlC6ULexcTTBEGkX1G3wdAfzPL\nvpibJEJjAmEpn2Epn3HLZ23A24H/NLM/A5uAvcxseY7zjXQRqYSyN3umTJnC8OHDAdhjjz04/PDD\nt09xbf9y63l+z59//vmo4inH87Vr1zJu3DjGjBmjfFbZc+Wz+OepVGr77Rvtvy9Dy6cbsAdwKXAi\nsCdwDPDbbq77eXcPshR1N92Aq4HXuugGnODuu5f42eoGlLwtXLiQ888/n69+9avceeedSYcjkphY\n1gbMZ8wq2GK33RSrh4CxWSZYvAR86O7jS/xsFSvJ26pVqxgxYgQ9e/Zk9erV7LXXXkmHJJKIKNYG\nJD17LsQ5heiqYtwL7G9m+7cfMLO9Se+p9avAMUiJ2rsNatV+++3HqaeeWrEVLWo9n5WmfMatmIVs\nH7W088ystcPKEQ+Y2eT2cwLH2VWFXgC8BFxrZj0zXZbXAMuBG0N8cEtLi77EkreOK1qoVS71JpVK\n0dLSUpZrF9MN2Bv4d6Cpi1N+A0xy960lBWb2E+BkYD9gD+BF0i2sozte28wagNnAUZnXXwK+4+6r\nS/n8zLXVDSgF2bZtG8OHD2f16tU8/vjjTJgwIemQRCquHN2A+e4U3NFVwBHA94D7gbWZ40NJL2V0\nOenVzH9YSmDufkWe560Dzinls0RC2WWXXbj88stZvXo1I0eOTDockZpRTMvqdeDr7v5sF6+PA+50\n99EB4kuUWlZhpVKp7dNepXTKZ1jKZzixTLAY0FWhAnD3PwIlTRmPicasRETyE9uY1Z+Ao9z9z128\nPhh4xt33CRBfotSyEhEpXCwtq8XAPWZ2ZOcXMl2A9wAPlRqYiIhIu2KK1dWkZ+g9Y2ZrzOyPmcda\n4A/APplzRHZQr92pK1euLMt16zWf5aJ8xq2Y+6zeAT4H3AL0Jz0z8AigH/AL4HOZc0Tq2pYtWzj6\n6KMZM2YM7777btLhiFS1onYKdvd17n4RMBAYknkMdPep7l5T/yo1wSKceptp1bt3b/baay+2bNnC\nwoXht4Crt3yWm/JZuqgmWNQTTbCQUj3wwANMmjSJ0aNH89prr9HNztciNSGWCRYiRanHFuqpp57K\nsGHDeOONN1iyZEnQa9djPstJ+YybipVIGe2yyy5ceOGFAMybNy/haESql7oBc1A3oISwatUqTjrp\nJC655BL+9m//NulwRMouiv2s6omKlYTi7hqvkrqhMasEaDZgOPWcx3IUqnrOZzkon6Ur52zAshQr\nM/uHclw3CS0tLZrSKiKSh8bGxjinrpvZINI3Bu9wGPiDuzeUElgM1A0oIlK4KPazyhSoOcCZQJ8u\nTtNveJEutLW1sXXrVvr06eqfj4h0Vkw34HzgFOAu0lvI/zDLY2OoAKV2aEwA7rnnHg444ADmzp1b\n8rWUz7CUz7gVs1PwBOAv3f2Nrk4ws0nFhyRSu3r37s1bb73F/Pnz+bu/+zvNEBTJUzH7Wb3k7oeU\nKZ6oaMxKQtu2bRvDhw9n9erVPPHEE5q8IzUplqnr15jZ1FwnmNkzRcYTHU1dl5C0ooXUsugWsjWz\nc4FLgD8C64G2ji8Dl7t71W9tr5ZVWKlUSi0J0itaDB8+nF69evH222+z1157FXUd5TMs5TOcWGYD\nngbMIz0T8OguTtNveJEu7Lfffpx66qmsW7eOtWvXFl2sROpJMWNWLwOvAj8D1gJbspz2/+k+K5Gu\nffzxx+y6665JhyFSFlG0rIChwOHuvq2rE8zsuuJDEql9KlQihSlmgsVLpLewz+XpIq4rNU4TVcJS\nPsNSPuNWTLH6FvAzMzs4xzl3FRmPiIjITooZs3oL2A34DPAx8B47zwbc192L6WKMisaspFK0hYjU\nkpjGrJ4mXZS6MqS4cOLTvuq6prRKOXz44YfMmjWL3/3udzz55JMqWFLVUqlU2bpTi2lZrXX3nMUo\nn3OqgVpWYek+lp2VsqKF8hmW8hlOLCtYfDOPc04p4roidafjihbz589POBqReBXTsrqV9E2/i9z9\n0bJEFQm1rKQSOq5osXr1agYNGpR0SCIliaVldR7QC3gnZCAi9ap9RYstW7awcOHCpMMRiVIxxWqd\nu09295eDRyM1TfexdG3atGmYGcuWLcv7PcpnWMpn3IqZDfiymQ119zVdnWBmC939vBLiEqkrp512\nGsuXL2f48OFJhyISpWLGrA4G/hH4kbtn3QpEswFFROpXOcasSr0p+BPgXXa+KXgfd+8VKsikqFiJ\niBQulgkWQ4GXgf8A/gCsAFZ1eKwEPg0Un9QQjQmEpXyGpXzGrZgxqw3uPiHXCWa2tsh4oqMVLERE\n8hPbChYndXd/lZkd5u4vlBRZBNQNKEmYPXs21113HQ0NDTQ0NDBjxgyam5uTDkskb1GsDZjPjcC1\nUKhEktDa2sqPfvQj3n//fdasSU+4bZ/OroIl9ayYMSss7TwzazWzpZnHA2Y2OXSAUjs0JtC9G264\ngffff3+HY8uWLWPu3Lk7nat8hqV8xq3glpWZ9Qb+HWjq9NIYoNnMzgYmufvWAPGJ1JXNmzdnPb5p\n06YKRyISl2JaVlcBRwDfAw4EBmQeY4ErMq/9fagApXZokkr3+vTpk/V43759dzqmfIalfMat2FXX\nT3P3f3b3N9z9o8zjdXe/DmgG1B0oUoQZM2YwatSoHY7tu+++TJ8+PaGIROJQTLEa4O7PdvWiu/8R\n2L34kKRWaUyge83NzcyZM4empiYOP/xwJk6cyI033ph1coXyGZbyGbdi7rPaamZ7u/ufs71oZoOB\nbaWFJVK/mpubNfNPpJNi7rO6ifT41GWdW1hmNg6YA7zq7lODRZkQ3WclIlK4WNYGHAz8HtiX9J5W\n7atVDAX2Jr3k0ufdver3u1KxEhEpXBRrA2aK0OeAW4D+pGf/HQH0A34BfK4WCpWEpzGBsJTPsJTP\nuBV1U7Bd+bFgAAAS2klEQVS7r3P3i4CBwJDMY6C7T3X3d83soJBBitSzbdu28cADD/Czn/0s6VBE\nElNwN2BeFzVb4+5Dg1+4wtQNKDF4/fXXOfDAA+nfvz9r165lwIABSYckklMU3YCZQL5oZrea2W/M\n7PFOjydI73UlIgGMGTOGE044gY0bN7Jo0aKkwxFJRMHFysy+A9wLnAMcAozs9BhRzHVj1dLSor7s\nQJTH4l188cUA3HTTTduPKZ9hKZ+lS6VStLS0lOXaxcwGfBO4Gfipu2/p4hxtay87SaVSWtKmSJs2\nbWLYsGFs2LCBZ555hnHjximfgSmf4cTSDfgXwLVdFaqMbxQZj9Qw/SIoXt++fTn33HMBWLhwIaB8\nhqZ8xq2YFSxeIT0L8L0c59RMN6BILC699FKOPPJIvvKVryQdikjFFVNUvg383MwOzXHO7UXGIzVM\nYwKlGT16NJMnT6Zfv36A8hma8hm3YlpWDwC7AWeZ2UZgPdDW4XUDBgWITUREBChugsVm4GnSRakr\nn3f3nTfgqTKaYCEiUrhyTLAopmW1wd0n5DrBzNbmel1ERKQQxYxZ5bOx4ilFXFdqnMYEwtm4cSNX\nXXUVGzduTDqUmqHvZ9yKWcj20Y7PzezzWc55oZSgRCS3SZMmce2113LPPfckHYpIRZS8NmCt3ACc\njcasJFY333wzU6dO5a/+6q946qmnkg5HZAdR7Ge10wVUrEQq7qOPPmLIkCF89NFHvPLKKxx88MFJ\nhySyXSwrWIgURWMC4ey2227bV1zouF6gFE/fz7ipWIlUqb/+678G4Je//CWbNm1KOBqR8grRDTje\n3Z8OFE9U1A0osfvud7/LSSedRFNTEz179kw6HBEg0jGramRmC4CmDoeWuPvXs5ynYiUiUiCNWYXj\n7j6kw2OnQiXhaUwgLOUzLOUzbsWsYFETzOxaoCfpgn2du69JOCQREelCsG5AMzsb+CvgVeAX7v6J\nmR0AnASsc/d7i7jmEOBW4AvuHqwVaGZfAn7r7u+Y2cnAPOB/ufsnnc5TN6CISIGi7QY0s5nAPwIN\npLe7f97Mhrv7m6RXaS/4NnszO5P0grkjgS4rhpk1mNntZrY087jbzIZ1eH2Gmb2ReUwEcPf73P2d\nzH8/QrqFeWyhMYrEZP369bz22mtJhyFSFqFaKwcBB7r7We5+DPA1YI6Z7Qt8WuQ1ryDdKutyhXcz\n6w20F5uDM4+PgSfMbFcAd7/B3UdnHo9n3jem06U2A1W/SnzsNCYQVsd8Ll68mGHDhvHtb387uYCq\nnL6fcQtVrH7v7ttv9HD350kXrIuBEUVec7y7L+vmnPOAQ4Ar3b3N3duAK0m3xi7J8b7b2v/DzEaS\n3n/r90XGKZK4o48+mh49evDoo4+yfPnypMMRCS5UsVppZheY2Z/M7LMA7r7J3a8GDmXHzRnzkik8\n3fkysNLdV3R435+B/8q81pWXzGyRmf0UuB44y93XFRqjFKZ9xQUJo2M+99xzT8466ywAfvGLXyQU\nUXXT9zNuISdYjCTdynnQ3bd2eu1Ydy9qtc3MPVHnZptgYWZrgKXuPrHT8fuBie6+WzGf2eE6mmAh\nVePJJ5/k+OOPZ/DgwaxatYpevXolHZLUqZgnWJxMetzqgc6FCqDYQpWHQcCHWY5/APQ3sz5l+lwp\ngsYEwuqcz2OPPZaxY8fyzjvv0NramkxQVUzfz7gFuc/K3R8xs9OBB8zsFeBWd381xLW7++hyf8CU\nKVMYPnw4AHvssQeHH3749u6C9i+3nuf3/Pnnn48qnmp/3jmfS5Ys4ZRTTmH8+PEccMABicdXbc/1\n/Sz+eSqVYsGCBQDbf1+Glnc3oJkd4e7PdXNOL+Au4HR3D1IIu+kGXA281kU34AR3373Ez1Y3oIhI\ngZLuBvxedydkugAvArYUHVFhXiT7bMMRwEsVikFERMqskGK1dz4nuft7pO+NCqmr5s29wP5mtn/7\nATPbGxgL/CpwDFKi9m4DCUP5DEv5jFshxWqCmT1hZv/bzMabWa79CP5camCddNWcXEC6BXWtmfU0\nsx7ANcBy4MYQH9zS0qIvsYhIHlKpFC0tLWW5diFjVv8PWE261QLwEfCfwBPA48Af2wd4zOxf3X1y\nSYGZ/QQ4GdgP2IN0l58DR3eccWhmDcBs4KjM6y8B33H31aV8fubaGrOSqrZu3ToaGhqSDkPqTKL7\nWZnZXe7+VTMbDEzs8BieOeX/Af9Bunid6e7Hhww0CSpWUq3a2to4/fTTefjhh1m5ciVDhw5NOiSp\nI0lPsLgWwN3fcfdF7n6Ru48kPZnhQuDXwDjgp6RXXxfZgbpTw8qVzx49etCvXz+2bdvGrbfeWrmg\nqpi+n3HLu1i5+x+7OL7S3W9198nuvg/pxWTfCBVg0jRmJdVq6tSpANx88820tRW84plIwaIYsyro\nomZ31sLuu+oGlGrW1tbGyJEjWblyJb/5zW/4whe+kHRIUieS7gYsxP8p03VFJE89evTgoosuAmD+\n/PkJRyNSmrIUK3d/vRzXleqm7tSw8snn+eefz5AhQxg7dizqJchN38+4BVkSSUTiNGzYMP70pz/R\ns2eu2yJF4leWMataYWY+c+ZMGhsbty/eKCIi2aVSKVKpFLNmzUruPqt6pAkWIiKFq6YJFiI70ZhA\nWMpnWMpn3FSsROrMp59+mnQIIgVTN2AO6gaUWvLQQw/xgx/8gHPOOYfvfve7SYcjNUzdgCJStM2b\nN/Pss89y0003aRq7VB0Vq25ouaVwlMewCs1nc3MzgwcPZunSpTz11FPlCaqK6ftZunIut6Ri1Y2W\nlhZNW5ea0KtXLy644AIAbrrppoSjkVrU2NhYXWsD1gqNWUmtWb58OaNGjaJv376sWbOGPffcM+mQ\npAZpzEpESjJy5EhOPvlkhg0bxvLly5MORyRvKlZSMRoTCKvYfC5atIjXX3+dcePGhQ2oyun7GTcV\nK5E6M2jQIB566CGamppobGykqamJ1tbWpMMSyUkL2XajfYKFJlmUTjkMq9h8tra2ctlll7Fs2bLt\nx9r/u7m5OURoVUnfz9K1rw1YDppgkYMmWEgtampq4uGHH856fPHixQlEJLVGEyykqmlMIKxi87l5\n8+asxzdt2lRCNNVP38+4qViJ1Jk+ffpkPd63b98KRyKSP3UD5qBuQKlF2casRo0axZw5c+p6zErC\nKUc3oCZYiNSZ9oI0d+5cNm3aRN++fZk+fboKlURNLasc1LIKK5VKacZVQMpnWMpnOJpgISIidUkt\nqxzMzGfOnKn7rERE8tB+n9WsWbOCt6xUrHJQN6CISOHUDShVTfexhKV8hqV8xk3FSkREoqduwBzU\nDSgiUjh1A4qISF1SsZKK0ZhAWMpnWMpn3FSsREQkehqzykFjViIihdOYlYiI1CUVK6kYjQmEpXyG\npXzGTcWqGy0tLfoSi4jkIZVK0dLSUpZra8wqB41ZiYgUTmNWIiJSl1SspGLUnRqW8hmW8hk3FSsR\nEYmexqxy0JiViEjhNGYlIiJ1ScVKKkZjAmEpn2Epn3FTsRIRkehpzCoHjVmJiBROY1YiIlKXVKyk\nYjQmEJbyGZbyGTcVKxERiZ7GrHLQmJWISOE0ZiUiInVJxaob2iIkHOUxLOUzLOWzdOXcImSXsly1\nhpQr8SIitaaxsZHGxkZmzZoV/Noas8pBY1YiIoXTmJWIiNQlFSupGI0JhKV8hqV8xk3FSkREoqcx\nqxw0ZiUiUjiNWYmISF1SsZKK0ZhAWMpnWMpn3FSsREQkehqzykFjViIihdOYlYiI1CUVK6kYjQmE\npXyGpXzGTcVKRESiV5djVmbWF/i/pIt1L+Awdz8uy3kasxIRKVA5xqzqddX1a4BH3L0VwMyOSTge\nERHJoe66Ac2sH3A+8Bdm9n/N7Ebgw4TDqgsaEwhL+QxL+Yxb1MXKzIaY2WIzawt42eHA7sBId//f\nwA3A42Y2MOBniIhIQNGOWZnZmcA/A1uBUe7es4vzGoDZwLjMoZeA77j76szrM4DpmdemAR8BvwMO\ncfdXMuc8B1zj7v/W6doasxIRKVC93Wd1BXAS8DSQ9Yc2s97AI6TH3g7OPD4GnjCzXQHc/QZ3H515\nPA68nXl7x9baFqB3WX4KEREpWczFary7L+vmnPOAQ4Ar3b3N3duAK4GRwCXZ3uDua4D/BI4HMLM9\ngVHAklCBS3YaEwhL+QxL+YxbtMUqU3i682Vgpbuv6PC+PwP/lXmtK+cAp5rZbGA+cL67ryohXMnD\n888/n3QINUX5DEv5jFu1T10/FFia5fgKYGJXb8oUpi+WKSbpwn//938nHUJNUT7DUj7jFm3LKk+D\nyD7t/AOgv5n1qXA8RSm2+6GQ93V3bq7Xs72Wz7EkulVK+cxK5LOQ4/WSz9DfzWzH8/0Ol1s15jOW\n72a1F6uamKpXL8VqxYoVOWMIoZ6KVa3kM5ZipXzGXayinbrezswWAOe6+06F1cxWA6+5+8ROx+8H\nJrj77iV+dtzJERGJlJZb2tGLwNgsx0eQvt+qJKGTLSIixamWbsCuWjj3Avub2f7tB8xsb9IF7FeV\nCExERMqvWopVVy2cBaRbUNeaWU8z60F6kdrlwI0Vik1ERMos2mJlZj/JLIN0OuBm9pyZPWtmvdrP\ncfetwMnAp6TvrfovYDdgortvrECM+5jZfWZ2o5n92szGdf8u6YqZ9TKzK83sIzPbL+l4qpmZDTSz\nW8zsBjObY2YPmNnopOOqVmZ2vZn9wsyuM7NHzeycpGOqBWZ2eb5rv0Y7ZuXuV+R53jrSN/km4efA\nne6+yMwOBe4zs9F53tAsO5sGpIB/TDiOWrAvsNHdZwCY2beBW4Cd9m2TvGxx9wsBzOwvgSVmdof+\nrRfPzD4LNJLnrO7oZwPGKrNK+7vAnu7+QebYCuA8d9fSTSXI/KU1XKuKhGNmzcBcdx+ZdCzVzszO\nBWa4+1FJx1KtMj1kdwH/B3gh22zvzqLtBiynQFuPDCf919YHHY69kzleV8q0lUvdKlM+Twd+FvB6\nVSFkLs3sUDO7C5gBnFV6dNUnYD5nAnMoYC/BuitWma1Hnia92G2XzUozazCz281saeZxt5kNq1ig\nVUL5DKsc+cy0qnZz95+WJ+o4hc6lu7/o7l8FvgX8R73tgRcqn2Y2Hujv7qlCPr/uihWBth4BVgK9\nzWxAh7cOJr0uYT0JlU9JC5rPTKE6HZhcxphjFSSXZtajY17d/ffAJ2R2bqgjob6bk4A9Lb1L+z9k\n3vdzM/tKzk9397p6AD0y/7sAaOvinKmk97sa3uHY3sA24PIOx+4Hzsn892HAMjLjgPXyCJnPDq+1\nAfsn/bNVez5Jd1Vd3+H5nKR/vmrMJemu/Ts6vP4Z0uuPHpL0z1iN+ex0/vCurtX5UXctKw+79cil\nwFmZvxB+DHzdM/8P1IuQ+TSz48xsLukuhhYzOyNwuNELlc/M7NRFwFfNbK2ZrQUuCh9xvAJ+N98D\nzMxuNbN/Bv4VuMzdS14lp5oE/t2JmZ0AtJC+NWmumR2S68LRTl1PWF5bj7j722irkXzkm88ngSeB\n6ZUJq2p1m093fxHoleUc2VE+ufwQ+HoFY6pmeW/b5OlZ00uAKflcuO5aVnmqia1HIqJ8hqV8hqNc\nhlW2fKpYZVdXXXkVoHyGpXyGo1yGVbZ8qlhltx7Itr3IAOBjd99c4XiqnfIZlvIZjnIZVtnyqWKV\n3YuktxnpLMjWI3VI+QxL+QxHuQyrbPms92KlrUfCUj7DUj7DUS7Dqng+671YaeuRsJTPsJTPcJTL\nsCqez7orVlYFW49UE+UzLOUzHOUyrKTzqVXXRUQkenXXshIRkeqjYiUiItFTsRIRkeipWImISPRU\nrEREJHoqViIiEj0VKxERiZ6KlYiIRE/FSkREoqdiJSIi0VOxEhGR6KlYiYhI9HZJOgARyc7MZgD/\nC9gCzAQuAtqAo4D73P3fEgxPpKJUrEQiZGZjgfeAuaR3X90CfN/dN5vZJOBWQMVK6oa6AUXidDTw\nIHA48N/AP7j75sxrg0i3sETqhlpWIhFy94UAZjYBeMzd3+/w8gQglURcIklRy0okbhOBx9ufmFlf\n4K9RF6DUGRUrkUiZ2XBgf+CJDoe/CGwF/t3MjjOz8QmEJlJx6gYUiddEYK27L+1wbDxwj7tvNbMv\nAd8zsxOBBuDzwNOkZxAud/cFlQ5YpFzUshKJ1wjg9k7HFgHDzOwnwAJgINDP3e8A/gLoCSwF3q1g\nnCJlZ+6edAwiUiQz6wdscnc3s5eAE9x9Q9JxiYSmlpVIFXP3TzKFai+gh7tvsLT+SccmEpLGrESq\nmJl9EdiNdPffC5nDpwF/BDYmFZdIaGpZiVS3zwAHAZuB981sCrDN3d9JNCqRwDRmJSIi0VPLSkRE\noqdiJSIi0VOxEhGR6KlYiYhI9FSsREQkeipWIiISPRUrERGJnoqViIhET8VKRESi9/8DUH5MWreG\ne8wAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f20760a7d90>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "That looks good! See how for each quadrant we go right, the error drops two quadrants going down (and even a bit better!)."
     ]
    },
    {
     "cell_type": "heading",
     "level": 5,
     "metadata": {},
     "source": [
      "Dig deeper"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Let's re-do the spatial convergence, but comparing at a much later time, say $t=1000$. "
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nx_values = numpy.asarray([11, 21, 41, 81, 161])\n",
      "\n",
      "dt = 0.1\n",
      "error = numpy.zeros(len(nx_values))\n",
      "\n",
      "t_final = 1000.0 \n",
      "\n",
      "x = numpy.linspace(0,L,nx)\n",
      "\n",
      "for i,nx in enumerate(nx_values):\n",
      "    \n",
      "    dx = L/(nx-1)\n",
      "    x = numpy.linspace(0,L,nx)\n",
      "   \n",
      "    sigma = alpha*dt/(dx*dx)\n",
      "\n",
      "    nt = int(t_final/dt)\n",
      "    \n",
      "    A = generateMatrix(nx, sigma)\n",
      "\n",
      "    T = numpy.zeros(nx)\n",
      "    T[0] = 100\n",
      "      \n",
      "    T = CrankNicolson(T, A, nt, sigma)\n",
      "    \n",
      "    T_exact = T_analytical(x, t_final, 100, alpha, L)\n",
      "    \n",
      "    error[i] = L2_error(T,T_exact)\n",
      "\n",
      "plt.figure(figsize=(6,6))\n",
      "plt.grid(True)\n",
      "plt.xlabel(r'$n_x$', fontsize=18)\n",
      "plt.ylabel(r'$L_2$-norm of the error', fontsize=18)\n",
      "#plt.xlim(1e-4,1)\n",
      "#plt.ylim(1e-4,1)\n",
      "plt.axis('equal')\n",
      "plt.loglog(nx_values, error, color='k', ls='--', lw=2, marker='o');"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAGNCAYAAACv5F5GAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt0VPW99/H3lzuoVAVRKSiaIlauCt6oxxMRSiVaPFpv\nbVUQFVsFtbVesOcxPtZVUWlVKBwoFTitVq34tIdSaSs6HK/YggJGkYo3WgFFqCKSkJDv88dMaIBk\nktmZmb1n9ue11qzl7Nmz50PWmG9+123ujoiISJS1CjuAiIhIU1SsREQk8lSsREQk8lSsREQk8lSs\nREQk8lSsREQk8tqEHSDKzEzz+kVEAnB3y+b11LJqgrtn/KitrWXBggXs2LEj0Ptz8bjtttuK4jOz\ncc0g18jkPc09t6nzWvp6oTzC+ndE8ftZKN/Nps7JBRWrHLjuuus466yzmDx5cthRdiktLS2Kz8zG\nNYNcI5P3NPfcps5r6vV33323WZ8TdWF8N3P1uS29ZqF8NzP93GywXFXBYmBmHuTns3jxYoYPH07b\ntm1Zvnw5/fr1y0E6ibsxY8Ywd+7csGOI7MXMcHUD5ld5eTmJRCKj95x++ulceeWVVFdXM3bsWGpq\nanITTmJtzJgxYUcQ2U0ikaC8vDwn11bLKo2gLSuATz/9lH79+rFu3TruuusubrrppiynExGJJrWs\nCkjnzp2ZNWsWAI888ohaV5J1mbb4RQqZilUOfe1rX+Oxxx7jxRdfpE0brRIQEQlK3YBptKQbUEQk\nrtQNKCIisaRiJVKgNGYlcaJilWefffYZS5YsCTuGiEhB0ZhVGtkes9q0aRNDhgxh06ZNrFq1iiOO\nOCJr1xYRiQqNWRW4rl27csIJJ7Bt2zauuOKKnO2hJSJSbFSs8mzatGl07dqVxYsXM3v27LDjSAHT\nmJXEiYpVnnXr1o2pU6cC8P3vf59169aFnEhEJPpUrEJwwQUXcPbZZ7N161bmz58fdhwpUGHtVi4S\nBk2wSCOXi4I3bNjAsmXLKCsry8n1RUTCkosJFipWaWgHC4myRCKh1pVEkmYDiohILKlllYZaViIi\nmVPLqsitWLGCjRs3hh1DRCRyVKwi4uGHH2bIkCF897vf1WJhaRats5I4UbGKiFNOOYWOHTvyxBNP\n8Jvf/CbsOCIikaIxqzTyPWY1c+ZMrrrqKrp27crrr7/OQQcdlLfPFhHJFo1ZFbkrr7ySYcOGsWnT\nJiZMmBB2HBGRyFCxihAzY/bs2eyzzz48+eST2opJ0tKYlcSJilXEHHHEETzyyCNUVFTQs2fPsOOI\niESCxqzS0DorEZHMacxKRERiScVKpEBpzEriRMWqANTW1vLMM8+EHUNEJDQqVhHn7px55pkMGzaM\nP/7xj2HHkQjRjusSJypWTSgvLw+1u8XMOPXUUwG44oor+PTTT0PLIiKSTiKRoLy8PCfX1mzANKIy\nG7CmpoaTTz6Zv/71r4wfP57/+q//CjuSRIDuZyVRpdmAMdWmTRvmzJlD27ZtmTlzJk8//XTYkURE\n8krFqkD069eP//zP/wRgxowZIaeRKFCrSuJE3YBpRKUbsE51dTVz585l7NixtGnTJuw4IiINykU3\noIpVGlErViL1acxKokpjViIiEktqWaWhlpWISObUspK9bNq0iddffz3sGCIiOaViVcBWrFjBMccc\nw7nnnktlZWXYcSTPtDegxImKVQHr06cPXbp0YfXq1dx+++1hxxERyRmNWaVRCGNWL774Il/5yldo\n1aoVL730EkOGDAk7kojEnMasZC8nn3wy119/PTt37mTs2LFUVVWFHUlEJOtUrIrAHXfcwZe+9CUq\nKipYvHhx2HEkTzRmJXGibRCKQKdOnfjlL39JTU0Np5xySthxRESyTmNWaRTCmJWISNRozEpERGJJ\nxUqkQGnMSuJExaqI/eUvf6GmpibsGCIiLaZiVaTuvPNOTjzxRKZMmRJ2FMkR7bgucaJiVaSGDBmC\nu3PbbbexevXqsOOIiLSIilWRGjlyJJdddhlVVVWMHTuWnTt3hh1JskxjVhInKlZFbMqUKXTv3p2X\nXnqJ+++/P+w4IiKBaZ1VGsWwzur3v/89Z511Fr169eLNN9+kXbt2YUcSkSKn29rnWTEUK4CZM2fy\njW98gy5duoQdRURiQMUqz4qlWElxSiQSmhEokaQdLKRFFi5cyMiRIyktLWXkyJEsXLgw7EgiIs2i\njWxjYuHChVx77bWsXbt217G6/y4rKwsrlrSAWlUSJ2pZxcQDDzywW6GCZLGaOnVqSIlERJpPxSom\nGrspY2VlZZ6TSLZonZXEiYpVTLRv377B4x06dMhzEhGRzKlYxcTEiRMpKSnZ6/hRRx0VQhrJBo1Z\nSZxo6noaxTZ1feHChUydOpXKykq2bNnCypUr6dy5MxUVFfTo0SPseCJSJDR1XVqkrKyMRYsWkUgk\nePXVVxk9ejSffvopM2fODDuaBKAxK4kTTV2PKTNjxowZnHnmmYwbNy7sOCIiaakbMI1i6wYUEckH\ndQOKiEgsqViJFCiNWUmcqFjJbtatW8fGjRvDjiEishsVK9ll0aJF9O3bl6uuugqN1UWf1llJnMSy\nWJnZfWb2CzO718yeMrNvhZ0pCvr27QvAb3/7Wx599NGQ04iI/EssixWww93HufsNwCRgtpnF9Wex\nS8+ePZkyZQoA11xzDR9++GHIiSQdjVlJnGT8C9rM5pjZg2Y2PBeB8sHdb6z39Gigwt1rw8oTJZdf\nfjnDhw/n448/5uqrrw47jogIEGCdlZnVAg8Bk939tZyk+tdnHQrMAb7q7llt+ZjZAOCHwJHAee7+\nTgPnxHKd1XvvvUe/fv2orKxk1apVHH300WFHEpECEonb2pvZBnc/JJshGvmcc4ApQDVQ4u6tGzmv\nG/BTYHDq0CrgOnf/R+r1icCE1Gvj3f3pPd5/IvA4MNDdN+/xWiyLFcD8+fPp3bs3AwYMCDuKiBSY\nqCwKfs3Muqc7wczmBcxT343AcOAFoMF/tJm1A/5MctuoY1KPbcAzZrYPgLs/4O69U4+nzaxV3Wup\n15cC24FTs5C5aJx77rkqVBGnMSuJkyDFaiIww8yGpDnnqwHz1DfU3dc2cc6lQH/gJnevTY073USy\na+87jbznMGB23RMz6wIcAjT1WSIiEpIg3YDvAPsCXUi2SD4C6k9OMKCHu7fNSkCzucAlDY1Zmdki\noI+7H7HH8ZXANnc/uYH37Af8PJV9M/Bl4DfuPqeBc2PbDSgiElQuugGD7LrenTRdcymHBouTsQHA\n6gaOvwsMa+gN7r4VuDCHmYrW8uXLOe6448KOISIxFKQbcLO7n+bupY09gC1ZztmYrsDWBo5/CnQy\ns4bv5S4Zu+qqqxg8eDALFiwIO4qkaMxK4iRIsbq4Ged8LcB1g1AfXZ7UTV8fP348W7bk628REZGk\njLsB3f2pZpyzIlicjG0C9mvgeGeSY1ZVLf2AMWPG0KtXLwD2339/Bg0atGtPtrq/bOPwfMKECcye\nPZuKigquv/565s6dG6l8cXxedywqefQ8vs8TiQRz584F2PX7MtsC3XzRzAy4BDgfKEkd/hvwmLv/\nMnvxmpxg8SRwdAMTLFYBW919aAs/WxMs6nnzzTcZNGgQlZWVLFy4kFGjRoUdSUQiKBLrrFJrm/5A\ncmeJM4CjUo8yYJ6ZPWlmWZkJWE9jFeMJ4HAzO7xevoNJbqE0P8sZYq9Pnz7ccccdAPz4xz/Wzuwh\nq/vLViQOgoxZ3QwcC/wA6EOyy60zyQJxY+q1W7IVMKWxCj2X5I4Vk82sdWoz2ruAt4EZWc4gwPXX\nX8/dd9/Nk08+SbKBLSKSe0HWWa0BLnT35Y28Phh4xN17tyiY2d3ACJKLePcHVpJsYZ3o7tX1zqvb\nbmlI6vXdtltqYQZ1A4qIZKhg9gbM1/6BuaZiJSKSuUiMWQHVqXGhBpnZIUBN8Egi0hwas5I4CVKs\nFgGPm9leWxmkugAfB55sabCoKC8v1y+FJmzfvp3VqxvaSERE4iSRSFBeXp6TawfpBjwEWAr0BDYA\n61MvdQcOBt4HTnL3DVnMGQp1Azbtvffe46tf/So7duxg1apV7LvvvmFHEpGQRaIbMFWEjgceBDqR\nnP13LNAR+AVwfDEUKmme7t27s88++/Duu+9y8803hx1HRIpUoEXBu96cnCp+UOrpR8V2a3i1rJpn\nxYoVDBkyhJqaGp555pnddliQ3Km/e4VIlESiZWVmc8zsQTMbnrqH1MbUo6gKlTTfwIEDufXWWwEY\nN24c27ZtCzmRiBSbIBMsLgXakhyvEgFg0qRJDBgwgLfffptHH3007DixoFaVxElO1lkVC3UDZuaV\nV17hrbfe4rzzzgs7ioiEKBLdgMBrZtY93QlmNi9gHilgxx57rApVHmlJhcRJkGI1EZhhZkPSnPPV\ngHlERET2EuS29guBfYGXzWw78BFQf3KFkbyDb1EoLy+ntLRU4wMSOfpOStQkEomctfiDjFntAJ6n\n8Z3QIbkouENLgkWBxqxabs2aNfTs2ZOOHTuGHUVE8iQqY1Yfu/tp7l7a2APQfc+FefPmMWDAgJxt\nvxJ3GrOSOAlSrGrq1lmlOedrQQNJ8ejTpw/V1dXce++9LF26NOw4IlLAghSrL9LEOit3XxE4kRSN\nk046ie9973vU1tZy2WWXUVVVFXakoqIxK4kTrbNKQ2NWLbd9+3YGDRrEmjVrmDRpEnfeeWfYkUQk\nx6IyZqV1VtJsHTt25MEHH8TMmDZtGlu2aDgzWzRmJXGidVaSc1/5yleYMWMGy5cv54ADDgg7jogU\noCDdgO+QXGfVBWhsnVUPd2+brZBhUTegiEjmctENGGRRcHfgBdKvszo0WBwREZG9BSlWm939tHQn\nmNn6dK8XEu1gIVGl+1lJ1ERtB4vh7v5UE+cMLIbp6+oGzB1357XXXqN///5hRylYKlYSVbnoBmzR\nnYKLnYpVbtTU1HDOOeewaNEili1bpoIlUmSiMnUdS/oPM5tpZvNTx3qnjrXOZkApPm3atKFHjx5U\nV1czZswYqqurw44kIhEX5Lb2HYHFwHzgCmBE6qV9gdnAU2a2X9YSSlGaPHkyhx9+OMuXL+eee+4J\nO05B0joriZMgLav/BPoAVwMnAp8DuPsrQA/gY2BStgJKcdpvv/2YPXs2ALfffjsVFRUhJxKRKAsy\nweJvwLfdfWnq+Xp3P7Te612BF9z9qKwmDYHGrHJv/PjxzJo1i4svvpj//u//DjuOiGRBJCZYmNnH\n7t6l3vPdilXq2AfunnZLpkKgYpV7n376KdOnT+f666+nffv2YccRkSyIygSLHWbW6KJfMyth9x0t\nRBrVuXNnbr75ZhWqADRmJXESpFj9AXjMzL605wtmdgrwOLCgpcFERETqBOkGPBR4CegJvENyUkUF\nyftcHQS8S/K29h9mNWkI1A0oIpK5SHQDuvt64HjgQZKb2bYFBgHtgJ8DJxRDoapTXl6u7pY8++ST\nT1i3bl3YMUQkQ4lEgvLy8pxcu0U7WJhZK5KtKYCP3L2oxqrUssq/5cuXM3r0aHr16sWSJUto1SrQ\nuvVY0HZLElWRaFnV5+617r4x9SiqQiXh6NWrFzU1NTz33HNMmzYt7DgiEhHaGzANtazC8bvf/Y6z\nzz6bjh07smrVKkpKSsKOJCIZiFzLSiQXRo8ezTe/+U22b9/OuHHjqK1Vo10k7lSsJJIeeOABunXr\nxpIlSzTBpRH6uUicBLn5okjOdenShTlz5tCmTRuGDRsWdhwRCZnGrNLQmJWISOYiM2al+1mJiEg+\n6X5WIgVKY1YSJ7qflRSUN954A3XNisSP7meVhsasomXy5MnccsstHHPMMXTt2pX27dszceJEysrK\nwo4mIvXkYswqyGzAA+sKVUPcfZOZ7duCTCIN2rRpE+6+212F165dC6CCJVLkdD+rJmgj2+hYuXLl\nXsfWrl3L1KlTQ0gTPn0vJWpyuZGt7mfVhPLycm0WGhFVVVUNHq+srMxzEhFpSGlpaaSK1Q+Bw4A1\nZrYWONDMlpvZRuB/gS8At2UxowhAo3cT7tChQ56TRIP+iJI40f2spGBMnDhxr01tS0pKmDBhQkiJ\nRCRfdD+rNDQbMHoWLlzI1KlTqayspEOHDkyYMCG2kyt0PyuJqlzMBszJdktm9mV3fyPrF84zFSuJ\nMhUriapCKlYfuHv3rF84z1SsREQyF5V1VpjZ2cBooDvJMavdXiY5liWSN7W1taxZs4ajjz467Cgi\nkgNB9ga8DngC+BbQHzhyj8cRQa4rEtQnn3xCaWkpQ4cOZcOGDWHHyRuts5I4CVJUrgFuBvZ19+7u\n3mvPB7ApqylF0ujcuTP77LMPW7Zs4aqrrtLegSJFKMjegB8B3dIN5phZqbsnWpgtdBqzKhzr1q2j\nb9++bN26lYcffpiLLroo7EgisRWV+1lVAAfm4LoigfXs2ZOf/OQnAFxzzTVs3Lgx5EQikk1BuwGn\nm9nANOc8FDCPSGDjxo1jxIgRbN68ORb7BWrMSuKkydmAZvYOsGdf2H7AeWb2OcnxqfqLgQ3omrWE\nIs1kZsyePZsnnnhCu1qIFJkmx6zMrAp4gWQRaq6T3L3gN2zTmJWISObCWme12d1Py+SiZrY+YB4R\nEZG9NGfM6uLmXszMOqX+c2SwONGj+1lJVOl7KVGTy/tZBZm6Pt3dv9vIa9OAC4FL3X1hFvKFSt2A\nxWHjxo20bt2arl2LayhVewNKVEVib0AzW+/uDd4p2My6kWxVTXL3L2chX6hUrArf4sWLOf/88xkx\nYgSPPPJI2HFEYiEq66walbqP1a/Q3oASEUceeSRVVVU8+uijzJ8/P+w4IhJQs1pWZjaH5PR1Ay4A\nGvsTtQ1wNNDO3QdlK2RY1LIqDtOnT+fqq6+mW7duVFRUFE13oLoBJapC6wY0s+beVPFz4A1goru/\n2JJgUaBiVRxqa2s5/fTTSSQSXHTRRTz88MNhR8oKFSuJqsiPWRUbFavi8fbbb9O/f38qKytZuXIl\nffv2DTuSSNGKSrEa6+5zshkiqlSsisuvf/1rjjzySE488cSwo4gUtUgUqzhRsZIoUzegRFXkZwOK\niIjkglpWaahlJSKSObWsRLKsoqIi7Agi0gwqVhJb3//+9+nXrx8LFiwIO0og2htQ4qTJYmVmPzSz\np81Mu1JIUenRowcA48ePZ8uWLSGnEZF0mnM/q7XAD4FH3b3WzC5193l5SRcyjVkVt507d3Lqqafy\nwgsvcOmllzJ37tywI4kUhVCmrpvZBnc/pN7zJhcFm9mr2m5JCsGbb77JoEGDqKys5Pe//z1lZWVh\nRxIpeGFNsPjczE7O8LoHBwkjkm99+vThRz/6EQC33347hfTHicasJE6ac6fgh4DnzWwjUAkcZGZv\npznfgOLYKVRi4brrrmP79u1cffXVmGX1j0ERyZLmdAO2Ar4LnA4cAJwMNLVJ7Unu3iErCUOkbkAR\nkcxFYrulZo5ZFcVmtypWIiKZi8qi4IuzdI6ItIDGrCROMi5W7v6UJV1qZgvNbHXqscDMLq47J/tR\nw1FeXq5fCjG0Y8cO1qxZE3YMkYKSSCQoLy/PybWDdAO2A34HjGzklD8CX3f36hZmC526AePpgw8+\nYNSoUWzevJnXXnuNzp07hx1JpKBEpRvwZuBY4AdAH6Bz6nE0cGPqtVuyFVAk37p160a7du1Yt24d\nP/jBD8KOIyIEa1mtAS509+WNvD4YeMTde2chX6jUsoqviooKjjvuOHbs2MGf//xnhg8fHnakveh+\nVhJVUWlZdW6sUAG4+zJgv+CRRMLXt29fbrvtNgDGjRvH1q1bQ04kEm9BilW1mTW6Q4WZHQLUBI8k\nEg033ngjgwcP5v3334/kvoFqVUmcBOkG/DnJ8alr92xhpboA7wfecPcrspYyJOoGlFWrVvHXv/6V\nMWPGaHcLkWaKyqLgQ4ClQE9gA7A+9VJ3knsCvk9yB4sNWcwZChUriTKNWUlURWLMKlWEjgceBDqR\nnP13LNAR+AVwfDEUKhERiY6MW1a7vTm5b+BBqacfuXttVlJFhFpWIiKZi0TLqj53r3X3jalHURUq\nkca88847fPbZZ2HHEImVFhUrkbh5/PHH6d+/P7fcEv66d20DJnGiYiWSgd69e1NVVcW0adNYsmRJ\n2HFEYkPFSiQDAwcO5NZbbwXgsssuY9u2baFl0UxAiRMVK5EMTZo0iQEDBvD2228zadKksOOIxIKK\nlUiG2rVrx9y5c2ndujWzZs1i/fr1Tb8pBzRmJXGSk2JlZj/KxXVFouLYY49l5syZLFu2jEMPLfib\nYotEXkvXWXUluTB4t8PAX9y9W0uCRYHWWYmIZC4X66zaBAjRleT+f+cA7Rs5Tb/hRUQkazIuVsAs\n4N+Bx4B/ADsaOOf7LQklIk3T3oASJ0GK1WnACe7+t8ZOMLOvB48kUrhWrVpF//79w44hUnSCTLD4\ne7pCBeDuxwXMI1KQ3J1vf/vbDBw4kBdeeCEvn6lWlcRJkGJ1l5mlvVeVmf01YB6RgmRmHHbYYbg7\nY8eOZfv27WFHEikqgWYDmtklwHeAZcAmoP4mtgbc4O4Ff2t7zQaUTFRWVnLcccfxxhtvcMMNN3DP\nPffk9PM0ZiVRFZWbL44C5tP4TEAAd/fWLQkWBSpWkqmlS5cydOhQAJ5//nlOOumknH2WipVEVVSK\n1WvAG8A0kncJbmg24MtaZyVxddNNN3H33Xdz5plnsmDBgrDjiORdJNZZkbx9/SB3r2nsBDO7N3gk\nkcJ2++2384UvfIHrrrsu7CgiRSNIy2oJcKa7b01zzinu/lxLw4VNLSuJMnUDSlRF5U7BVwPTzOyY\nNOc8FjCPiIjIXoK0rN4B9gW6ANuAj9l7NmBPdw/SxRgpalmJiGQuSmNWL5AsSo0pmm2oy8vLKS0t\nVXeLtMi2bdvYsGEDJSUlYUcRyZlEIpGzW9cEaVmtd/e0xag55xQCtawkG9asWcOoUaPo0KEDy5Yt\no337dKs+mk9jVhJVURmz+nYzzvlagOuKFKUePXrQqlUrKioquOOOO8KOI1KQgrSs5pC8BcjD7v5U\nTlJFhFpWki3PPfccp556Kq1atWLp0qUMHjw47EgiOROVltWlQFtgQzaDiBSzU045hWuvvZadO3cy\nduxYduxoaC29iDQmSLH60N0vdvfXsp5GpIjdeeedlJSUsGrVKn7729+2+Hq5GsgWiaIgswFfM7Pu\n7v5BYyeY2Tx3v7QFuUSKTqdOnZg3bx4ffvgh7dq1Y+TIkVRVVdG+fXsmTpxIWVlZ2BFFIivImNUx\nwI+BO9y9wVuBaDagSOMWLlzItddey9q1a3cdKykp4f7771fBkqIQlY1s6y8K3g58xN6Lgnu4e9ts\nhQyLipXkwsiRI/nTn/7U4PFFixaFkEgku7QoWKQIVFVVNXi8srIyo+tonZXESZBitdndT0t3gpmt\nD5hHpOg1tii4Q4cOeU4iUjiCzAa8uBnnaFGwSCMmTpy417ZLBxxwABMmTMjoOmpVSZxk3LJqzkJg\nd18RLI5I8aubRDF16lQ2bNjAihUr+OSTT+jRo0fIyUSiK+MJFgBmZsAlwPlA3Z+IfwMec/dfZi9e\nuDTBQvLhmmuu4Wc/+xmDBg3i5Zdfpm3b5s1N0piVRFUkdrAws3bAH4A5wBnAUalHGTDPzJ40s4Kf\nCSiSL3fddRe9evXi1Vdf5a677go7jkgkBZm6/n+A7wL3AP8D1E2m6A58HbgBmO7u/zeLOUOhlpXk\nyzPPPMOwYcM48MADeeedd+jcuXPYkUQCi8o6qzXAhe6+vJHXBwOPuHvvLOQLlYqV5NOMGTMYNWoU\nhx9+eNhRRFokKsVqg7sf0tJzCoGKlUSZxqwkqiIxZgVUm9nBjb1oZocANcEjiYiI7C5IsVoEPG5m\nx+35QqoL8HHgyZYGE5H01KqSOAnSDXgIsBToSfKeVvUnWBwMvA+c5O4Ff78rdQNKmNydlStXMnDg\nwLCjiGQkEt2AqSJ0PPAg0Ak4NvXoCPwCOL4YCpVImKqrqznjjDM4/vjjee21hm8dp/tZSZwE6QbE\n3T9098uBA0luWnsocKC7X+HuH5nZl7MZUiRu2rZtS69evaiurmbMmDHU1GgYWOIt0A4WTV7U7AN3\n7571C+eZugElTFu3bqVfv368//773HnnnUyaNCnsSCLNEomp66kgZwOjSY5T7blbhQFD3b3hraUL\niIqVhO2pp55ixIgRtG3bluXLl9OvX7+wI4k0KRJjVmZ2HfAE8C2gP3DkHo8jglxXRPY2fPhwrrzy\nSqqrq5kyZcpur2nMSuIkyP2srgEmAT9x9x0NnaD7WYlkzz333MNRRx3FxIkTw44iEpogU9c/Arql\n6x8zs1J3T7QwW+jUDSgikrlIdAMCFSRnAWb7uiIiIg0KUlSuAaab2YA05zwUMI+INJPGrCROgoxZ\nLQD2Bc4zs8+BTUBtvdcN6JqFbCLSiM2bN7N+vYaGJT6CjFlVAS+QLEqNOcndO7QkWBRozEqiaPny\n5ZSVlXHYYYfx/PPP06ZNkL85RXInF2NWQb7lm939tHQnaDagSO6UlJTQpk0bXn75ZaZMmcJNN90U\ndiSRnAvSshru7k81cc5Ad1/RomQRoJaVRNWiRYs444wzaNeuHa+88grHHHNM2JFEdonMDha7XcDs\nJHd/KUt5IkXFSqJs1KhRPPnkk5xwwgnqDpRIiWqxWu/uh2YpT6SoWEmUffLJJ/Tr14+///3v/OlP\nf2LEiBFhRxIBVKzyTsVKou7pp5/GzDjttLTDyCJ5FZUJFiISAYlEgmHDhoUdQyQvtNOEiIhEXja6\nAYe6+wtZyhMp6gYUEclcVPYG3E0hFyozu8HMaps+U6RwvPLKK+zcuTPsGCJZFdtuQDPrB5QCajpJ\nQWpob8DJkyczZMgQ7rvvvvwHEsmhWBYrM2sL3AHcQvpto0QKSr9+/aitreWHP/whb775ZthxRLKm\nxWNWuy5k9k3gK8AbwC/cfbuZfQkYDnzo7k8EuOahwBzgq+6etcJqZj8CngLeBd5u7Noas5JCNGbM\nGObNm8ezxo8JAAANb0lEQVTJJ5/Ms88+S+vWrcOOJDETyTErADO7Dfgx0I3k7e5fNbNe7v4WyV3a\nHw9wzXNIbph7JGm66sysm5k9ZGarU4/fmNkX670+0cz+lnqcbmZDgU7FcHNIkYb89Kc/pXv37rz4\n4ovcf//9YccRyYpstVa+DPRx9/Pc/WTgAuB+M+sJBB3pvZFkq6zRHd7NrB3wZ5LrxY5JPbYBz5jZ\nPgDu/oC79049FgNfBw4wsxnAj1LXmW5m3wiYUyQUjd3P6oADDmDWrFkA3HfffVRVVeUxlUhuZGtR\n8FJ3r6x74u6vmtkFwK3AHwJec6i715qlbUleCvQHRrt7LYCZ3QT8A/gOcO+eb3D3m+v+28x6Ad90\n9+8GzCgSSWVlZUyfPp1zzz2X9u3bhx1HpMWyMmaV6rLbH7gdOMPdX6v32njgZ+4eqDCa2VzgkobG\nlcxsEckW3RF7HF8JbEu18hq77r8DY4GLgenALHdftcc5GrMSEclQJPcG3HUhsyNJtnL+4O7Ve7x2\nirs/F/C6c2m8WH0ArHb3YXsc/x9gmLvvG+Qz611HxUpEJENRnmAxguS41YI9CxVA0ELVDF2BrQ0c\n/xToZGbq/5Ci1diYlUgxysqYlbv/2czOAhaYWQUwx93fyMa1m/roXH/AmDFj6NWrFwD7778/gwYN\norS0FPjXLws91/Mwnr/66qsZnb948WLeeustxo8fH4n8el48zxOJBHPnzgXY9fsy25rdDWhmx7r7\nK02c0xZ4DDgr6BhVA9ecS+PdgP8A3mykG/A0d9+vhZ+tbkApCp9//jmnn346K1asYMWKFfTu3Tvs\nSFLEwu4G/EFTJ6S6AC8HdgROlJmVwBENHD8CWNXAcZFY6tSpEyUlJWzfvp3LLruM2lptiSmFJZNi\ndXBzTnL3j0mujcqmxpo3TwCHm9nhdQfM7GDgaGB+ljOIREpdN0xzPfDAAxx88ME899xzTJs2LTeh\nRHIkk2J1mpk9Y2a3mtlQM0u3h8vGlgbbQ2PNybkkW1CTzay1mbUC7gLeBmZkOYNIQTvwwAOZOXMm\nADfffDNvvfVWyIlEmi+TYrWVZOvqDuA5YIuZPWlmN5rZENt99W6L+xjM7G4zewU4C3Aze8XMlqfG\nxYBd3Y4jSO6S8XrqsS/JaeuftzSDSJTVDXRnYvTo0XzrW99i+/btzJihv+ekcGQyweIxdz/fzA4B\nhtV79Eqd8gnwv8AzwDnufmr24+aXJlhIMfr44495/PHHueKKK2jVKpY3XpAcC3VRsJkNdvdlDRw/\nnH8VrtOA7oC7e8Fv9axiJVGWSCQCta5Eci3U2YANFarU8ffcfY67X+zuPUhuJvu3bAUMW3l5ecYD\n2SIicZRIJCgvL8/JtbO23dJuFzV7xN0vzPqF80wtKxGRzEV6b8DdLmp2lLuvyfqF80zFSuLigw8+\noKamhsMOOyzsKFIEwl4U3GzFUKhEoi5b3dNLliyhb9++XHLJJVosLJGlqUAiMde3b1/atWvHkiVL\nNJ1dIisn3YDFQt2AEhdPPPEE5557Lvvssw8rV67kyCOPDDuSFLCC6QYUkcJyzjnncOGFF7Jt2zbG\njRun7kCJHBUrkQKV7SUVU6dO5aCDDuLZZ5/lL3/5S1avLdJSWbmNRzErLy+ntLRUiy+l6HXt2pVf\n/epXdOnShcGDB4cdRwpQIpHI2bpUjVmloTErEZHMacxKRERiScVKpEBpGzCJExUrEUnrxRdfRN3h\nEjYVK5EClY9JPzfeeCNDhw5l1qxZOf8skXRUrESkUccffzwAN9xwA++9917IaSTOVKxEClQ+xqzO\nO+88vvGNb/DZZ59x+eWXqztQQqNiJSJp/exnP6Nr16489dRT/PznPw87jsSU1lmloXVWIkmPPfYY\nF1xwAX379mXFihW0bl3wNwKXHMrFOivtYNEE7WAhkuwO3LJlCxdddJEKlTRKO1iERC0ribJEIqE/\noiSStIOFiIjEklpWaahlJSKSObWsRCQyqqqqePnll8OOITGhYiVSoMLcG3Dz5s0cd9xxlJaWcuqp\np1JaWsrIkSNZuHBhaJmkuGk2oIhk7IADDqBz585s376dZ599dtfxtWvXAlBWVhZWNClSalmJFKgw\nZwKaGe3bt9/r+Nq1a5k6dWoIiaTYqViJSFZVVlaGHUGKkIqVSIEK+35WDbWsADp06JDnJBIHKlYi\nEsjEiRMpKSnZ7VhJSQkTJkwIKZEUM02waIK2W5KoCvs7WTeJYurUqVRWVtKhQwcmTJigyRUxpu2W\nQqJFwSIimdOiYBHZJewxK5F8UrESEZHIUzdgGuoGFBHJnLoBRUQkllSsRAqUxqwkTlSsREQk8jRm\nlYbGrEREMqcxKxERiSUVK5ECpTEriRMVKxERiTyNWaWhMSsRkcxpzCoE5eXl6m4REWmGRCJBeXl5\nTq6tllUaallJlCUSidB3XhdpiFpWIiISS2pZpaGWlYhI5tSyEhGRWFKxEilQmvgjcaJiJSIikacx\nqzQ0ZiUikjmNWYmISCypWIkUKI1ZSZyoWImISORpzCoNjVmJiGROY1YiIhJLKlYiBUpjVhInKlYi\nIhJ5GrNKQ2NWIiKZ05iViIjEkopVE3TzRYkqfS8lanTzxZCoG1CiTDdflKjKRTegilUaKlYiIpnT\nmJWIiMSSipVIgdKYlcSJipWIiESexqzS0JiViEjmNGYlIiKxpGIlUqA0ZiVxomIlIiKRpzGrNDRm\nJSKSOY1ZiYhILKlYiRQojVlJnKhYiRSoV199NewIInmjYiVSoP75z3+GHUEkb1SsYiKMLqNcfGY2\nrhnkGpm8p7nnNnVeXLr5wvp3RvH7WSjfzUw/NxtUrGJCxapl14hisXr33Xeb9TlRp2LVsvfHpVhp\n6noaZqYfjohIALqflYiIxI66AUVEJPLahB1ARLLPzA4E7gU+Axw4Evieu/8t1GAigJndB+wHbAEG\nAXPc/aF071HLqoXMrK2Z3WRmn5nZYWHnEUnpCXzu7hPd/Vrgj8CDIWcSqbPD3ce5+w3AJGC2maWt\nRypWLTceSACdQs4hsou7r3D3a+odegf4Ylh5ROpz9xvrPT0aqHD32nTvUbFqIXef5u5Lw84hxcPM\nDjWzRWaW9n/eDJ0FTMvi9SSGsvndNLMBZvYYMBE4r6nzY1uscvQLQaRFzOwc4AWSY0yNTtU1s25m\n9pCZrU49fmNmDbaczKwM2Nfdf5Kb1BIH2f5uuvtKdz8fuBr439Q4a6NiWaxy8QtBJEtuBIaT/H42\nuE7FzNoBfyY5QeqY1GMb8IyZ7bPHuWUkW1UX5zCzxENWvptm1qr+9zTVM7UdODXdh8eyWJHlXwgi\nWTTU3dc2cc6lQH/gJnevTfX130Tyj6/v1J1kZucBI9z9Knd3M7s/Z6klDrL13TwMmF33BjPrAhwC\npL12XItV1n4h7CGrK7YlfpoaZE45F3jP3d+t976NwOup1zCzAcDDwPlmtt7M1gOXZz+xxEW2vpvA\nx4CZ2RwzmwL8ErjW3Velu3As11m15IduZnU/9HsBzOzfgPNJdieWm9lv3f132U8tsssAYHUDx98F\nhkFyPABom8dMItC87+ZW4MJMLxzLYtVMTf7QAdz9WeBZYEJ+YonQFdjawPFPgU5m1t7dq/KcSQRy\n+N2MazdgczT5Q89zHpE62tBToipn300Vq8bpF4JE1SaSW9XsqTOwTa0qCVHOvpsqVo3TLwSJqpXA\nEQ0cPwJIO0gtkmM5+26qWDVOvxAkbI217p8ADjezw+sOmNnBJLetmZ+PYBJ7ef9uqljpF4JEV2NL\nIeaS/INpspm1Tm0AehfwNjAjT9kk3vL+3VSx0i8EiRAzu9vMXiG564Sb2StmttzMdk1Dd/dqYASw\nk+T6ldeBfYFh7v55GLml+IX93YzlnYLN7G6SP9DDgP1Jdvk5cGLqh113Xjfgp8CQ1OurgOvc/R95\nDy0iEmOxLFYiIlJY1A0oIiKRp2IlIiKRp2IlIiKRp2IlIiKRp2IlIiKRp2IlIiKRp2IlIiKRp2Il\nIiKRp2IlIiKRp2IlIiKRp2IlIiKRp2IlIiKR1ybsACLSMDObCPQFdgC3AZcDtSTvAvD/3P3REOOJ\n5JWKlUgEmdnRwMfAVJK3sNkBTHL3KjP7OjAHULGS2FA3oEg0nQj8ARgE/BP4kbtXpV7rSrKFJRIb\nalmJRJC7zwMws9OAxe6+pd7LpwGJMHKJhEUtK5FoGwY8XffEzDoAZ6IuQIkZFSuRiDKzXsDhwDP1\nDp8NVAO/M7N/M7OhIUQTyTt1A4pE1zBgvbuvrndsKPC4u1eb2X8APzCz04FuwEnACyRnEL7t7nPz\nHVgkV9SyEomuI4CH9jj2MPBFM7sbmAscCHR0918DXwBaA6uBj/KYUyTnzN3DziAiAZlZR6DS3d3M\nVgH/7u6bw84lkm1qWYkUMHffnipUBwGt3H2zJXUKO5tINmnMSqSAmdnZwL4ku/9WpA6PApYBn4eV\nSyTb1LISKWxdgC8DVcAWMxsD1Lj7hlBTiWSZxqxERCTy1LISEZHIU7ESEZHIU7ESEZHIU7ESEZHI\nU7ESEZHIU7ESEZHIU7ESEZHIU7ESEZHIU7ESEZHI+/9LvfphxbpK8gAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f2075c031d0>"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "error"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "array([  3.59569224e-03,   1.34923182e-03,   4.90905474e-04,\n",
        "         1.76021088e-04,   6.26681531e-05])"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Wait, convergence is not that great now! It's not as good as second order, but not as bad as first order. *What is going on?*\n",
      "\n",
      "Remember our implementation of the boundary conditions? We used\n",
      "\n",
      "\\begin{equation}\\frac{T^{n}_{N-1} - T^{n}_{N-2}}{\\Delta x} = q.\\end{equation}\n",
      "\n",
      "Well, that is a **first-order** approximation! \n",
      "\n",
      "But, why doesn't this affect our solution at an earlier time? Initially, temperature on the right side of the rod is zero and the gradient is very small in that region; at that point in time, errors there were negligible. Once temperature starts picking up, we start having problems.\n",
      "\n",
      "**Boundary conditions can affect the convergence and accuracy of your solution!**"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "---\n",
      "###### The cell below loads the style of the notebook"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.core.display import HTML\n",
      "css_file = '../../styles/numericalmoocstyle.css'\n",
      "HTML(open(css_file, \"r\").read())"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Arvo:400,700,400italic' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=PT+Mono' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Shadows+Into+Light' rel='stylesheet' type='text/css'>\n",
        "<link href='http://fonts.googleapis.com/css?family=Nixie+One' rel='stylesheet' type='text/css'>\n",
        "<style>\n",
        "\n",
        "@font-face {\n",
        "    font-family: \"Computer Modern\";\n",
        "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
        "}\n",
        "\n",
        "#notebook_panel { /* main background */\n",
        "    background: rgb(245,245,245);\n",
        "}\n",
        "\n",
        "div.cell { /* set cell width */\n",
        "    width: 750px;\n",
        "}\n",
        "\n",
        "div #notebook { /* centre the content */\n",
        "    background: #fff; /* white background for content */\n",
        "    width: 1000px;\n",
        "    margin: auto;\n",
        "    padding-left: 0em;\n",
        "}\n",
        "\n",
        "#notebook li { /* More space between bullet points */\n",
        "margin-top:0.8em;\n",
        "}\n",
        "\n",
        "/* draw border around running cells */\n",
        "div.cell.border-box-sizing.code_cell.running { \n",
        "    border: 1px solid #111;\n",
        "}\n",
        "\n",
        "/* Put a solid color box around each cell and its output, visually linking them*/\n",
        "div.cell.code_cell {\n",
        "    background-color: rgb(256,256,256); \n",
        "    border-radius: 0px; \n",
        "    padding: 0.5em;\n",
        "    margin-left:1em;\n",
        "    margin-top: 1em;\n",
        "}\n",
        "\n",
        "div.text_cell_render{\n",
        "    font-family: 'Alegreya Sans' sans-serif;\n",
        "    line-height: 140%;\n",
        "    font-size: 125%;\n",
        "    font-weight: 400;\n",
        "    width:600px;\n",
        "    margin-left:auto;\n",
        "    margin-right:auto;\n",
        "}\n",
        "\n",
        "\n",
        "/* Formatting for header cells */\n",
        ".text_cell_render h1 {\n",
        "    font-family: 'Nixie One', serif;\n",
        "    font-style:regular;\n",
        "    font-weight: 400;    \n",
        "    font-size: 45pt;\n",
        "    line-height: 100%;\n",
        "    color: rgb(0,51,102);\n",
        "    margin-bottom: 0.5em;\n",
        "    margin-top: 0.5em;\n",
        "    display: block;\n",
        "}\t\n",
        ".text_cell_render h2 {\n",
        "    font-family: 'Nixie One', serif;\n",
        "    font-weight: 400;\n",
        "    font-size: 30pt;\n",
        "    line-height: 100%;\n",
        "    color: rgb(0,51,102);\n",
        "    margin-bottom: 0.1em;\n",
        "    margin-top: 0.3em;\n",
        "    display: block;\n",
        "}\t\n",
        "\n",
        ".text_cell_render h3 {\n",
        "    font-family: 'Nixie One', serif;\n",
        "    margin-top:16px;\n",
        "\tfont-size: 22pt;\n",
        "    font-weight: 600;\n",
        "    margin-bottom: 3px;\n",
        "    font-style: regular;\n",
        "    color: rgb(102,102,0);\n",
        "}\n",
        "\n",
        ".text_cell_render h4 {    /*Use this for captions*/\n",
        "    font-family: 'Nixie One', serif;\n",
        "    font-size: 14pt;\n",
        "    text-align: center;\n",
        "    margin-top: 0em;\n",
        "    margin-bottom: 2em;\n",
        "    font-style: regular;\n",
        "}\n",
        "\n",
        ".text_cell_render h5 {  /*Use this for small titles*/\n",
        "    font-family: 'Nixie One', sans-serif;\n",
        "    font-weight: 400;\n",
        "    font-size: 16pt;\n",
        "    color: rgb(163,0,0);\n",
        "    font-style: italic;\n",
        "    margin-bottom: .1em;\n",
        "    margin-top: 0.8em;\n",
        "    display: block;\n",
        "}\n",
        "\n",
        ".text_cell_render h6 { /*use this for copyright note*/\n",
        "    font-family: 'PT Mono', sans-serif;\n",
        "    font-weight: 300;\n",
        "    font-size: 9pt;\n",
        "    line-height: 100%;\n",
        "    color: grey;\n",
        "    margin-bottom: 1px;\n",
        "    margin-top: 1px;\n",
        "}\n",
        "\n",
        ".CodeMirror{\n",
        "        font-family: \"PT Mono\";\n",
        "        font-size: 90%;\n",
        "}\n",
        "\n",
        "</style>\n",
        "<script>\n",
        "    MathJax.Hub.Config({\n",
        "                        TeX: {\n",
        "                           extensions: [\"AMSmath.js\"],\n",
        "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
        "                           },\n",
        "                tex2jax: {\n",
        "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
        "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
        "                },\n",
        "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
        "                \"HTML-CSS\": {\n",
        "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
        "                }\n",
        "        });\n",
        "</script>\n"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 23,
       "text": [
        "<IPython.core.display.HTML at 0x7f2075a99710>"
       ]
      }
     ],
     "prompt_number": 23
    }
   ],
   "metadata": {}
  }
 ]
}