07_Parabolic_PDEs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<script>\n",
       "code_show=true; \n",
       "function code_toggle() {\n",
       " if (code_show){\n",
       " $('div.input').hide();\n",
       " } else {\n",
       " $('div.input').show();\n",
       " }\n",
       " code_show = !code_show\n",
       "} \n",
       "$( document ).ready(code_toggle);\n",
       "</script>\n",
       "The raw code for this Jupyter notebook is by default hidden for easier reading.\n",
       "To toggle on/off the raw code, click <a href=\"javascript:code_toggle()\">here</a>.\n",
       "<style>\n",
       ".output_png {\n",
       "    display: table-cell;\n",
       "    text-align: center;\n",
       "    vertical-align: middle;\n",
       "}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "HTML('''<script>\n",
    "code_show=true; \n",
    "function code_toggle() {\n",
    " if (code_show){\n",
    " $('div.input').hide();\n",
    " } else {\n",
    " $('div.input').show();\n",
    " }\n",
    " code_show = !code_show\n",
    "} \n",
    "$( document ).ready(code_toggle);\n",
    "</script>\n",
    "The raw code for this Jupyter notebook is by default hidden for easier reading.\n",
    "To toggle on/off the raw code, click <a href=\"javascript:code_toggle()\">here</a>.\n",
    "<style>\n",
    ".output_png {\n",
    "    display: table-cell;\n",
    "    text-align: center;\n",
    "    vertical-align: middle;\n",
    "}\n",
    "</style>\n",
    "''')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![title](figs/title.png)\n",
    "\n",
    "<h1><center>Module 07: Numerical Solutions to Parabolic PDEs</center></h1>\n",
    "\n",
    "This module focuses on generating numerical solutions of **parabolic** partial differential equations (PDEs).  Recall that these types of second-order PDEs\n",
    "\n",
    "$$Au_{xx}+2Bu_{xy}+Cu_{yy}+Du_{x}+Eu_{y}+Fu+G=0 \\tag{1}$$\n",
    "\n",
    "are defined by the condition $B^2-AC = 0$ (where $y$ is considered the $t$ variable and thus $A=E=1$ and all others are equal 0).\n",
    "\n",
    "Three interesting examples (with increasing complexity) are the following:\n",
    "\n",
    "   * **Heat Flow / Diffusion (Homogeneous medium)**: Used to model time-varying heat flow (or diffusive) phenomena with the rate of evolution is determined by the thermal conductivity $\\kappa$ (or diffusivity). The solution $\\phi=\\phi(x,y,z,t)$ is known on the boundary of the solution domain $D$ (e.g., $\\left. \\phi\\right|_{\\partial D}=F(x,y,z)$) as is the initial distribution $\\phi_0=\\phi_0(x,y,z,t=0)$:\n",
    "   \n",
    "$$\\frac{\\partial \\phi}{\\partial t} = \\kappa \\nabla^2 \\phi. \\tag{3a}$$ \n",
    "   \n",
    "   \n",
    "   * **Heat Flow / Diffusion (Heterogeneous medium)**: This is like the example above, except now that the thermal conductivity (or diffusivity) is spatially varying, $\\kappa = \\kappa(x,y,z)$:\n",
    "\n",
    "$$\\frac{\\partial \\phi}{\\partial t} =  \\nabla \\cdot \\left(\\kappa \\nabla \\phi \\right). \\tag{3b}$$ \n",
    "   \n",
    "   \n",
    "   * **Convection - Diffusion**:  These types of problems involve the types of heat flow (or diffusion) described above along with material transport at velocity ${\\bf v} = [v_x,v_y,v_z]$ like we investigated in Module 4:\n",
    "   \n",
    "$$\\frac{\\partial \\phi}{\\partial t} + {\\bf v} \\cdot \\nabla \\phi =  \\nabla \\cdot \\left(\\kappa \\nabla \\phi \\right). \\tag{3c}$$ \n",
    "\n",
    "In the section below, we will look at generating solutions $\\phi$ of the 1D and 2D heat flow / diffusion equation and then the 2D convection-diffusion equation throughout a computational domain $D$.  We will also examine different boundary conditions (BCs) on boundary of D (defined by shortform $\\delta \\Omega$), namely the Dirchelet BC\n",
    "\n",
    "$$ \\phi(\\mathbf {x}) = f(\\mathbf {x}), \\quad \\forall \\mathbf {x}\\in \\delta \\Omega$$\n",
    "\n",
    "and the Neumann BC\n",
    "\n",
    "$${\\frac {\\partial \\phi}{\\partial \\mathbf {n} }}(\\mathbf {x} )=f(\\mathbf {x} )\\quad \\forall \\mathbf {x} \\in \\partial \\Omega $$\n",
    "\n",
    "where ${\\bf n}$ is a normal vector defined at the grid boundary that points inward."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1D Heat Flow\n",
    "\n",
    "Let's start with one of the more straightforward examples: 1D heat flow in a thin uniform metal bar with constant thermal conductivity $\\kappa$.  Let's say that the bar is of length $L$ in the $x$ direction, and is much thinner than $L$ in both of the $y$ and $z$ dimensions.  In this case, we can rewrite equation 3a with the following explict **partial differential equation (PDE)**:\n",
    "\n",
    "$$\\frac{\\partial \\phi}{\\partial t} = \\kappa \\frac{\\partial^2 \\phi}{\\partial x^2}. \\tag{4}$$ \n",
    "\n",
    "For simplicity, let's say that the ends of the bar are fixed at some $\\phi(x=0,t)=A^\\circ$C and $\\phi(x=L,t)=B^\\circ$C, which represent homogeneous Dirichelet **boundary conditions**. Finally, let there be an **initial condition** where the temperature given by $\\phi(x,t=0) = \\phi_0(x)$. We can also say that the bar is perfectly thermally insulated along its length, and thus can only lose or gain heat at its ends. Thus, we have sufficient information to solve this heat flow problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical solution domain\n",
    "\n",
    "To compute a numerical solution, let's first **discretize** the solution space.  In this case, we will segment the $x\\in[0,L]$ into $nx+1$ discrete element spaced $dx = L/nx$ apart.  We will also break our time axis $t\\in[0,t_{max}]$ into $nt+1$ steps where $dt=t_{max}/nt$.  This will give us the following example solution grid:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAacAAAGNCAYAAAC12L7IAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnXu4HVV5/z/vyQWMcOQWraIhGPECFJQEAUslXrnYVlSoRQlqfzaNgVqVqtAqEKvFtioXFRURKRcvFFFRULQgIioiIMQLggkQQVDkKohKIO/vjzWH7Ozsfc5Mzt5nvTPz/TzPPOfsNd939nftPTNr1qxZ7zZ3RwghhIjESG4DQgghRDdqnIQQQoRDjZMQQohwqHESQggRDjVOQgghwqHGSQghRDhq1TiZ2VwzczM7bQrey83skmG/j6iOmb3UzL5nZvcU39OXcnuqK2Z2s5nd3KN81MxOLNY/XHzOzy7WzTCzZWb2CzP7U7Fu/ynwqmOyCzM7pvhcFlaIucTMws8hmj7IjZnZNODvgYOBPwc2Be4Bfg1cAZzn7ucN8j03lLED0t3n5nUiqmBmc4EvA/cCnwZ+B/w8g49nAX8HPBt4DvCUYtUMd3+4T8zNwDYdRQ+T/N8GXAOcB3zZ3R+q4GMucFNX8R+B+4EbgSuBz7v7d8pus+C/gH8EvgqcATxCOo4BDgeOAi4FzgZWk+E7iECPk/wa0jlvOfApdz9r6l01g4E1TkXD9FVgH9KJ43zgVmALYB7wGuCZpAOwDjwLeDC3CbEeLwY2Bg53989k9LE36QT9CPALUoOwccnYE0jHyAgwCjwDeAXpou4XZrbI3X9Q0c99wPHF/9NJx93OwJuAQ83sG8Ah7v6brrgX9dneXwE3uPtf91n3APCSKg3pAIh8TC4r/s4gfZ/7Ay8ws/nu/rZ8turLIHtOB5EapmuBvdz9vs6VZjYL2G2A7zdU3L2VV4I14EnF39uyuoCvAd8Hlrv7H3r0isbjeHe/ubPAzB4H/DvwT8CFZrZ7xX3wXnc/prvQzJ4KfAp4KfB1M9vD3f84tt7dV/bZ3pNIPaN+6+6a4oYp9DHZ/dmb2YuAbwJvMbMTu79vUQJ3H8gCnAQ48JaKcRsBR5C6wQ+SbnN8B/jbHtq5xXuc1lV+SapKz+2/voh5ffF6YfG613JaR5wDl/TY3uOAY4HrSVfL9wAXAi/uoR17r2NIt3/OJ10xPwh8G3hexc/qb4CLgNuBP5FO0N8Glvb6PIrP9r2k2z5/AlYCRwMze2x7f+BM4Abg96Qr46uANwMjffzMAt5JunV0fxFzHXAi8IQe2iNJt6/Gtv994KCSdR/ve1vYodsOOB34FfBQ8RmdDmzXY5vHjMWTevY/KHzdvAH7/83FtqaX0MwdR/M/heZLJd93bqHv67n47K+jx/FZeLq54/XYvtO9XAKc1mfdzd37+zj1v7mrbGaxj11NOpYeLHRfpuuYIuYx6fQ/9/ysWH9g13d1GvB04PPAHaRbgYPYh18H/Aj4Q7HdU4E/6xFzyTie9wYuAO5k7Tnjv4HN+n2fwCbAccAtxXtfA+xfaKYD/8rauwsrgcPKfLaD7DndVfx9etkAM5tJ2on2It2z/ijpQDoA+LyZPdvd/3WAHiF9mMuAtxSvj+9Yd80EfjcDvgtsD/ywiN0K+FvgG2b2Jnf/RI/QBcA7SCfjU4A5wKuAi4o6Xj+RaTNbDHyCdN//K6Sd5/HATsAbSBcH3ZwN7AqcQxoXeDlpZ15gZn/jxd5T8H7SQfID0kHxOOCFpFtQuwKLuvxsDnyLdOvoetKB8BDpFu7fA+cCvym0mwEXk8Zmri60I6QD4TNmtoO7v2uCj+Bm0ve2kLS//E9RNrYOM9sV+D/SWOd5pJPDM4HXAi83sxe5+5U9tn048BLS5/qtou65WAYcAvyVmY26++8mu0F3f9DMPkDa917Luvt8N6eRTl5HA6uK15A+43uLv93Hzr2TsHca6a7LT0gn4D+QemZ7ku7E/N94wTmPyRJY8bd7XGoe6Ti7ATgLeAzponwy+/BbSb3jzwNfJ31+bwAWmtlu7v7bCc2aHUXa/+4mDdHcQTq//AuwX9Hr7t4fZ5B6iFuQLihmkr7PL5jZS4GlpDtmXyM1dgcCHzaz37r758c1VPUKcZwriOeQTk5rSAOorwS2mSDmSNIXdwEdV5ykk+7NxbrndZTPZZI9p/Gu4ia6SiM1Dl78ta4rnfuKD39uj6u0Xu//j0X5SSU/36uK7T++x7qten0epJ1/847yjUkHowOLumLm9djuCGuv5HfrWveZovxjdPWsSAfW4zpen1Zo39Gl25h0IK0Bnl3ycziGrt5SUW6s7R28tmvdq4vyn3d67djW74HnTHL/H9tfJ9VzKnS3FLoXlHjfsWOi77489v0WuodZ91jreRz02v9LxIzt78eUiSNdBKwh9byn9dBvOZEnMh6THZ68R/mLi7qtoTgPdnxXDvxHj5jJ7MMPde/DpN6Mkx7MWO/80FX2gkL7Pbp6Saw9hx7XZ3/+CrBRR/lfFuV3ky4YNutY99TC648m/GzLfgklv6i/Jd1y8o7lLuCLwF/30P+i+PKe2WPd/yviT+0oG/tyT5vow+7xwXbviDdToXEiXSH8nnT7aose+n8vYo7qcSBc1kM/g9SbubLkZ3tV8f6bl9BeQo8GqMvTt0q+7y496vV40oMAtwGPnSB+S9IJ8Yd91u9cbP+/SvoZOxgXdpX/xdjB1SfuO8X65/fY1nFl3nsCX2MH6iAap8sL3Xq3tntox46Jvvtyodu445h8fJen9WK79/8e9egVM7ZvHVMmjvQwiJN6Pjae/16ech+THZ682JeOAd5HulPxcFH+oR7f1a/pOJkPaB/+VA/940i92j+wbuNxCes3Tl8strNDn/f+EXBHn/2514XtjcW6F/ZY963i81nvgqRzGeij5O5+tpl9kdQK70nqTe1JGs/Y38xOJzUSbmabAk8DfuW9BzovLv4+Z5AeJ8EzSbccv+vud/dYfzHwLnr7Xa8b7u6rzew3wOYl3/8s4IPAT83s86T749/18bvr3+5R9h3SjrGOTzPbEng7sB/p6uaxXXFbd/y/K6lXdam7/34C37sC0wA3s2N6rJ9R/H3WBNuZiF2Kvxf3WX8xa/fJ7oH+Kyb53oOm3+2gQWxz0NvdYNz9d2b2FeCvgWvM7Auk/fMH7l7mqbzcx2QnR49thtQgfIfUYJzZQ3utu/+pR/lk9uH1jnV3v8/MriHdBn8W4w9b7EFqmA80swN7rJ8JzDazLd39ro7ye733QzW3AduSLqq7+RXpnPBnxf89GWjjBOkLBr5RLGOPmL+KNM5wCKmF/hJr7+vf3mdTY+WbDdrjBjIZv/3uyT9M+pImxN0/ZGZ3ku7hvpl039/N7NvA2733fejux4Zx90fM7C5S7wd49L79D0k70xWke/93F/42A/6Z9HDFGGN17LtjdbBl8XfXYunHJiW2NR6T+X5+3aMsJ2NPJE44TrAB23yE9MBAFF5NeqjmNax9HPuPZnYO8C++/qPvnWQ9Jjtxd5tY9Sj99rfJ1Kff5zT2XhONo25Jag+OnkC3CWufL4B067QXD0NqIPutY+2FaU+GniHC3R9x97NJ9z8hDbLD2kr9WZ/QJ3bpxmMNgJn1amwH1bgN0u8G4e6nu/vupB3pZaRHhJ9PevT48T1CntBdUFwsbEkxAFvwRlLDtMzdd3P3pe7+Lk+Px/YatBw7sLfusa6bsc/jOHe3cZYXlNhWmffZkO8nRE8CwMyeBjyZdAD3uurcUMY+36u8zyThAbCm+Nvvone9E6S7/8Hdj3H3p5MeSjgYuKz4e84E75f9mNxA+u1vk6nPesd617Ym+gzuA+6Z4Bg1d181wXYGxlSmL7q/+GsA7n4/6bHCrc1sux76sYPp6hLbHrsSfEqPdQv6xDxCtSuk60mPmz67eFKtmyp+J4W73+vuF7j7P5AeNtiCNAjZzV49yv6SdPL4UUfZ04q/Xyi5jStIJ6Lnm1n37b9+2l7+BslYfRb2WT9WPvTvZ5IcVfz9SnGMTJpijuHhxcthZizoexwWje64F4rufounjAp7k8aj9yxuN/cjzDE5ICazD693nBZz555NeoT7ugne+3JgczPbYUKXU8TAGiczO8jMXmJm623TzP4M+IfiZee90lNJjdV/F1f0Y/qtgHd3aCZibMzgHzoLi4lwB/WJuYt0D/UxJbaPpwmHZ5G6te/pep95pFttq0lPKg4cM9unT89wrMfU6x79uzsPWjPbmDQfBFLqnzFuLv4u7HrP55CeqFyHYpzrc6QruQ90f+dmtklxYODud5A+twVm9u5edTCzeWa2bQ//Vfgu6WS1p5kd0LX9A0g9zBtIV+XhGMtlR3pk/17S3L9BbHdb0lyeZ5JOfr0eqx4UPyf1yF/e2ZMvjrETe3ibbWa9JuY/lvTE58OkJ7t6kvuYHAKT2YcXFcdrJ8eQequf7TPG1cnYna1PmtmTulea2WPNbPeJqzA4BjnmtBtpbOLXZnYZa/N9bUu6BfUY0nPwnV31DwD7kubfXGtmF5AGOA8knXT/y93LnEw+TRrMP9LMdibNDXh6se0vksa8urmINAbydTO7lPTI6bXu/pVx3ucIUg/gsGI+wrdYO6diU9LkspvGiZ8MnyPdi7+M1JhY4WVX0u2fXvNBriM9QNE5z2ke6WTVecCeTvr8jjezF5CuWrcjpak5lzQu0M1hwI7AEtJcigtJJ5JtSVe+f0N6KmhMux3pBLKoqMNvSOMgzyrqcBDr54grTfGQzetIcy4+b2ZfJp0sx1LJ3E9K37NmnM2UpriA+kBH0VbF30915Ft7f5+Hfd5iZveSvsOx9EXPJ52UbwAOdvcbKlrarOOBk+mkQf2dSQPdI6RH9l9X4iS1wRQPFJxAurD8UfFw1HTSHLLbWD+rx9bA5WZ2Hak3cAvp8/gr0u2oE0v0HnMekwNlkvvw14DvmtnZpLGpPYvlZkpc6Lj7RWZ2BOni9RfFufgmUsO/Dalndhlp7tnUMN6jfFUWUlf+UFJjcD3pCuoh0gd1Aeke8nqZBkiPuP4raRLeH0hfwGX0yBxAn0fJi3U7FO8zlqngEtIH+np6P0r+WNIcnVtZ+9jnaR3rez5KS7o18Z+kE/ifSFe53wRe2kO7kAqP1k7w+S4pPtsbSb2ku0lXwu8ANu3SXlK8b3eGiBtJA569HmPdnjTp7w7S47lXkcaixvvMHwv8G2uze9xPujA4nq75WKSnfQ4jzaMYm3/yS9JFwlvomtMyzudwDD0eJe9Y/wxSw3s7qUG+nZT54hlVtzWBj7HPZbxlYVfMzV3rVxff448LzwfQI3vHBvj4Y/E9Xg58GNhznPie+2C//X+i/ZbU4B5BumX/UPEd/xfponOdONKxdBTpKbRfFfvE7cX+exBdj5f380SmY7LDk1f8rtY7lgaxD5POddeQzqO/JV20P7FHzCX9PJMatLNJFxIPFdu5BvgQsKDCfjDee5xGiSkVVohFg7D0swJ7ebUniIQQIgy1+j0nIYQQ7UCNkxBCiHCocRJCCBEOjTkJIYQIh3pOQgghwjHw3Ho52GqrrXzu3Lm5bQghRK246qqr7nT32bl99KIRjdPcuXO58speeU+FEEL0w8ymLFdeVXRbTwghRDjUOAkhhAiHGichhBDhUOMkhBAiHGqchBBChEONkxBCiHCocRJCCBGOVjZOK+9eydLzlzJ67Cgjy0YYPXaUpecvZeXdK3NbE0KICWnDOawRufUWLFjgZSfhfu0XX+OA/z2A1Y+sZvWa1Y+WzxiZwYxpMzjnwHPYd7t9h2VVCCEmxSDPYWZ2lbsvGJbXyTClPSczO9XM7jCzn/RZb2Z2opmtMLPlZrbLIN9/5d0rOeB/D+DB1Q+u86UCrF6zmgdXP8gB/3tA36uPs846i7lz5zIyMsLcuXM566yzxn2/tukjelKd8+sjeoqmLxsz2XNYrajyk9CTXYDnA7sAP+mzfj/ga6Sfet4d+EGZ7c6fP9/L8KavvslnvGeGcwx9lxnvmeGHnn/oerFnnnmmz5o1a52fw541a5afeeaZPd+rbfqInlTn/PqInqLpq8RM5hzWC+BKn8I2oMoy9W8Ic8dpnD4BHNTx+nrgiRNts2zjtOl/bDrulzq2jB47ul7sNttss86OM7Zss802Pd+rbfqInlTn/PqInqLpq8RM5hzWCwI3TlM+5mRmc4GvuvuOPdZ9FXi/u19WvL4IeKe7rzegZGaLgcUAc+bMmb9q1cT5C0eWjeBMXN8RG+GRox5Zt2xkhF6flZmxZs2a9bfRMn1ET6pzfn1ET9H0VWImcw7rhcacymM9ynp+E+5+srsvcPcFs2eXy/i+ycxNNlg3Z86cnlqV539v1TlueURP0cqrxEzmHFY7prqrRsbbehpzinHfvCn6iJ6i6SN6iqavEqMxp3yN08tY94GIK8pss2zjtOKuFT7rfbPG/WJnvW+Wr7hrRc/4M88807fZZhs3M99mm23G3dnaqI/oSXXOr4/oKZq+bMxkz2HdRG6cpnTMycw+CywEtgJ+AxwNzABw94+bmQEfAfYBHgTe4D3Gm7rRPCchRFtoyzyn1k3ChTRX4LjLj+OM5WfwwEMPsMnMTVi00yLeuvtbmbfFvCE6FUKIyTOoc5gapyFTtXESQggRu3GK9rTelNCGvFRCiObShnNY63pOGnMSQtSZtow5tarnpNx6MfKDNUkf0VM0fURP0fRlY5Rbr2aL5jnl10f0pDrn10f0FE1fJUbznGq2KLdefn1ET6pzfn1ET9H0VWKUW69mlB1zUm694ekjelKd8+sjeoqmrxKj3HoNRbn1hlce0ZPqnL88oqdo5VVilFuvZovGnPLrI3pSnfPrI3qKpq8SozGnmi3KrRdDH9GT6pxfH9FTNH3ZGOXWqxma5ySEaAua59RQ9t1uX5YvWc7i+YsZ3WiUERthdKNRFs9fzPIly9UwCSFC05ZzWOt6TkIIIRLqOQWjDXmphBDNpQ3nsNb1nDTmJISoMxpzaiDKrRcjP1iT9BE9RdNH9BRNXzZGufVqtmieU359RE+qc359RE/R9FViNM+pZoty6+XXR/SkOufXR/QUTV8lRrn1aoZy6+XXR/SkOufXR/QUTV8lRrn1Gopy6w2vPKIn1Tl/eURP0cqrxCi3Xs0WjTnl10f0pDrn10f0FE1fJUZjTjVblFsvhj6iJ9U5vz6ip2j6sjHKrVczNM9JCNEWNM+pobQlL5UQopm05RzWup6TEEKIhHpOwWhDXiohRHNpwzmsdT0njTkJIeqMxpwaiHLrxcgP1iR9RE/R9BE9RdOXjVFuvZotmueUXx/Rk+qcXx/RUzR9lRjNc6rZotx6+fURPanO+fURPUXTV4lRbr2aodx6+fURPanO+fURPUXTV4lRbr2Gotx6wyuP6El1zl8e0VO08ioxyq1Xs0VjTvn1ET2pzvn1ET1F01eJ0ZhTzRbl1ouhj+hJdc6vj+gpmr5sjHLr1QzNcxJCtAXNc2oobclLJYRoJm05h7Wu5ySEECKhnlMw2pCXSgjRXNpwDmtdz0ljTkKIOqMxpwai3Hox8oM1SR/RUzR9RE/R9GVjlFuvZovmOeXXR/SkOufXR/QUTV8lRvOcarYot15+fURPqnN+fURP0fRVYpRbr2Yot15+fURPqnN+fURP0fRVYpRbr6Eot97wyiN6Up3zl0f0FK28Soxy69Vs0ZhTfn1ET6pzfn1ET9H0VWI05lSzRbn1YugjelKd8+sjeoqmLxuj3Ho1Q/OchBBtQfOcGkpb8lIJIZpJW85hres5CSGESKjnFIw25KUSQjSXNpzDprznZGb7ACcA04BT3P39XesfB5wJzAGmAx9w90+Pt02NOQkh2oLGnIaAmU0DPgrsC2wPHGRm23fJDgV+5u47AwuBD5rZzEG8v3LrxcgP1iR9RE/R9BE9RdOXjVFuvWE9tw57ABd2vD4SOLJLcyRwEmDAtsAKYGS87WqeU359RE+qc359RE/R9FViNM9peI3TAaRbeWOvFwEf6dJsCnwLuB14AHjZRNtVbr38+oieVOf8+oieoumrxCi33pAwswOBvd39jcXrRcBz3f2fOjQHAH8BvA2YB3wT2Nndf9e1rcXAYoA5c+bMX7Vq1YTvr9x6w9NH9KQ659dH9BRNXyVGufWGx63AUzpePxm4rUvzBuDcomFfAdwEPLN7Q+5+srsvcPcFs2fPLvXmyq03vPKInlTn/OURPUUrrxKj3HrDu603HbiRNJY0E7gW2KFL8zHgmOL/JwC/ArYab7sac8qvj+hJdc6vj+gpmr5KjMachttA7QfcAKwE/q0oWwIsKf5/EvAN4MfAT4CDJ9qmcuvF0Ef0pDrn10f0FE1fNka59WqG5jkJIdqC5jk1lLbkpRJCNJO2nMNa13MSQgiRUM8pGG3ISyWEaC5tOIe1ruekMSchRJ3RmFMDUW69GPnBmqSP6CmaPqKnaPqyMcqtV7NF85zy6yN6Up3z6yN6iqavEqN5TjVblFsvvz6iJ9U5vz6ip2j6KjHKrVczyo45Kbfe8PQRPanO+fURPUXTV4lRbr2Gotx6wyuP6El1zl8e0VO08ioxyq1Xs0VjTvn1ET2pzvn1ET1F01eJ0ZhTzRbl1ouhj+hJdc6vj+gpmr5sjHLr1QzNcxJCtAXNc2oobclLJYRoJm05h7WucRpjja9Zpwu5xns/RSOEEBFp+jlMt/UKdFtPCFEHdFuvgSh9UYwULE3SR/QUTR/RUzR92RilL6rZokfJ8+sjelKd8+sjeoqmrxKjR8lrtih9UX59RE+qc359RE/R9FVilL6oZih9UX59RE+qc359RE/R9FVilL6ooSh90fDKI3pSnfOXR/QUrbxKjNIX1WzRmFN+fURPqnN+fURP0fRVYjTmVLNF6Yti6CN6Up3z6yN6iqYvG6P0RTVD85yEEG1B85waSltSfwghmklbzmGt6zkJIYRIqOcUjJV3r2Tp+UsZPXaUkWUjjB47ytLzlzZjVrUQovG04RzWup6TxpyEEHVGY04NRLn1YuQHa5I+oqdo+oieounLxii3Xs0WzXPKr4/oSXXOr4/oKZq+SozmOdVsUW69/PqInlTn/PqInqLpq8Qot17NUG69/PqInlTn/PqInqLpq8Qot15DUW694ZVH9KQ65y+P6ClaeZUY5dar2aIxp/z6iJ5U5/z6iJ6i6avEaMypZoty68XQR/SkOufXR/QUTV82Rrn1aobmOQkh2oLmOTWUtuSlEkI0k7acw1rXcxJCCJFQzykYbchLJYRoLm04h7Wu56QxJyFEndGYUwNRbr0Y+cGapI/oKZo+oqdo+rIxyq1Xs0XznPLrI3pSnfPrI3qKpq8So3lONVuUWy+/PqIn1Tm/PqKnaPoqMcqtVzOUWy+/PqIn1Tm/PqKnaPoqMcqt11CUW2945RE9qc75yyN6ilZeJUa59Wq2aMwpvz6iJ9U5vz6ip2j6KjEac6rZotx6MfQRPanO+fURPUXTl41Rbr2aoXlOQoi2oHlODaUteamEEM2kLeew1vWchBBCJNRzCkYb8lIJIZpLG85hres5acxJCFFnNOY0JMxsHzO73sxWmNkRfTQLzewaM/upmX17UO+t3Hox8oM1SR/RUzR9RE/R9GVjlFtvWM+twzRgJfBUYCZwLbB9l2Yz4GfAnOL14yfaruY55ddH9KQ659dH9BRNXyVG85z6Ny4zgd2BVwKvBfYG5laI3wO4sOP1kcCRXZqlwHur+FJuvfz6iJ5U5/z6iJ6i6avEKLdeB2Y2DXgF8EZgL1IDZR0SB34FfBb4pLuvGGdbBwD7uPsbi9eLgN3c/bAOzfHADGAHYFPgBHc/vce2FgOLAebMmTN/1apV49YDlFtvmPqInlTn/PqInqLpq8Qot16Bmb0K+DlwJvAn4F3AS4CdgaeTelGvAc4B9geuM7NPmtkT+m2yR1n3Jz0dmA+8jNQze7eZPX29IPeT3X2Buy+YPXv2eNV4FOXWG155RE+qc/7yiJ6ilVeJUW69tbfYbgPeAmxWphsG7AacC7yrz/oyt/WOAI7peP0p4MDx3ldjTvn1ET2pzvn1ET1F01eJ0ZjT2oZh4w3aKGzUp3w6cCOwLWsfiNihS/Ms4KJCOwv4CbDjeO+n3Hox9BE9qc759RE9RdOXjVFuvSFiZvsBx5Oe3DvV3d9nZksA3P3jhebtwBuANcAp7n78eNvUPCchRFtoyzyn0o1TMe6zmbtfUbx+DHAUsCPpVt1HhuZyAqqmL1p590qOu/w4zlh+Bg889ACbzNyERTst4q27v5V5W8wbolMhhJg8gzqHNaVx+gZwjbu/o3j9IeAw4MfATsBb3P2jwzI6HsqtJ4QQ1YncOFXJELET8F0AMxsBDgHe6e7zgfdSPNZdB9qQl0oI0VzacA6r0nP6I/Bid7/MzOYDV5Am4N5iZnsBX3X3TYfotS8acxJCtIW2jDlV6Tn9Bnha8f9LgZXufkvxehPg4UEaGwbKrRcjP1iT9BE9RdNH9BRNXzZGufV6Pwb+YeB24APAncD7OtYdAVyV65FDzXPKr4/oSXXOr4/oKZq+SozmOfVunB4LfJL0AMSpwKyOdd8Djs1VCeXWy6+P6El1zq+P6CmavkqMcuvVjLJjTsqtNzx9RE+qc359RE/R9FVilFuvoSi33vDKI3pSnfOXR/QUrbxKjHLrrb1d92XgOWW7YcDGwNuAJVPZ/dOYU359RE+qc359RE/R9FViNOa0trH5MCkb+Q+ANwO7ANO7NE8iZST/FHAPKRfec6eyEsqtF0Mf0ZPqnF8f0VM0fdkY5dbrwMzmkTKTvxZ4XNGq/65otDYn/faSkeY9fQw4w91731gdEprnJIRoC22Z51RlEu5M0k9e7EbqLW0M3EX6vadL3X3iX/sbEsqtJ4RoE8qtVxOUW08IIaoTuXFq1dN6Y7QhL5UQorm04RzWup6TxpyEEHWmLWNOreo5KbdejPxgTdJH9BRNH9FTNH3ZGOXWq9mieU759RE9qc759RE9RdNXidE8p5otyq2XXx/Rk+qcXx/RUzR9lRjl1qsZyq2XXx8aCeG6AAAaRUlEQVTRk+qcXx/RUzR9lRjl1uuDmT3HzM41szvN7GEz26Uo/w8z22c4FgeHcusNrzyiJ9U5f3lET9HKq8Qot16PBdgT+CPwM+BEYA2wS7HuvcCXcnX/NOaUXx/Rk+qcXx/RUzR9lRiNOfVunC4jJYI1YHpX4/RK4Je5KqHcejH0ET2pzvn1ET1F05eNUW69HpjZg8Ar3f3rZjYNWA0scPerzez5wIXu/phSGxswmuckhGgLmue0Pn8EZvVZ90TgvsnbGT77brcvy5csZ/H8xYxuNMqIjTC60SiL5y9m+ZLlapiEEKFpyzmsSs/pPGAz4AVF0Wpgvrv/yMy+Adzp7q8Zjs3xUW49IYSoTlN6Tu8m/Z7TtcX/DrzOzL4F7A4sG7y94dCGvFRCiObShnNYpXlOxaPj/w08H5hGeijiO8Db3P1HQ3FYAo05CSHagsaceuDuV7v7i4BNgScDo+7+gpwNUxWUWy9GfrAm6SN6iqaP6CmavmyMcuvVbNE8p/z6iJ5U5/z6iJ6i6avEaJ5TPzE8izTedDJwetfyP7kqodx6+fURPanO+fURPUXTV4lRbr0emNkhwKnFh3YH8FCXxN39qaU2NmCUWy+/PqIn1Tm/PqKnaPoqMcqt15t3kzJEzHb3rd19264lS8NUBeXWG155RE+qc/7yiJ6ilVeJUW693rf07gdelLur12vRmFN+fURPqnN+fURP0fRVYjTm1Ltx+jpwWG7DvRbl1ouhj+hJdc6vj+gpmr5sjHLr9cDMngacCxwLfAO4p1vj7r1vqA4ZzXMSQrQFzXNan1uBHwFnkh6IWN21dD8gEZK25KUSQjSTtpzDqvSczgBeDXwF+Dk9GiN3z5LCSLn1hBCiOpF7TtMraF8OvN3dTxiWmali5d0r+eD3P8iZy8/kgYceYJOZm3DwTgdz+B6HM2+LebntCSHEuLThHFal53Q7cIi7f3O4lqqjMSchRFvQmNP6fBrI8pMYg0K59WLkB2uSPqKnaPqInqLpy8Yot17vR8n/EbgJ+CZwOPD33UuuRw41zym/PqIn1Tm/PqKnaPoqMZrn1LtxWjPB8kiuSii3Xn59RE+qc359RE/R9FVilFuvB2a2zUQad19VamMDRrn18usjelKd8+sjeoqmrxKj3Ho9cPdVEy3DNDoIlFtveOURPanO+csjeopWXiVGufVqtmjMKb8+oifVOb8+oqdo+ioxGnMaWwk3AjsX/99UvO63rMxVCeXWi6GP6El1zq+P6CmavmyMcusVmNmngfe4+01mdlrRovfF3d8w3vphoXlOQoi20JZ5TqUfiIhM1fRFK+9eyXGXH8cZy894dHb1op0W8dbd39qY2dVCiOYyqHNY5MZp3AcizOxGM9t5qsxMJWt8zTpdyDV5EqoLIcQG0fRz2ES39dYAu7v7FVNnqTq6rSeEaAttua1XJX3RQDCzfczsejNbYWZHjKPb1cweMbMDBvXeSl8UIwVLk/QRPUXTR/QUTV82RumL1j6ttwbYdVBPXwDTgJXAU4GZwLXA9n10FwMXAAdMtF09Sp5fH9GT6pxfH9FTNH2VGD1Kvm7jdAFweonlfyZ8M9gDuLDj9ZHAkT10bwEOBU4bZOOk9EXD00f0pDrn10f0FE1fJUbpiwqKMadfA3/qK1qLu/tTxxMUt+j2cfc3Fq8XAbu5+2Edmq2BzwAvBD4FfNXdz+mxrcXAYoA5c+bMX7Vq4gQVSl80PH1ET6pzfn1ET9H0VWKUvmhd9nf3bUss4zZMBdajrPuTPh54p7uP+8m6+8nuvsDdF8yePbvEWyt90TDLI3pSnfOXR/QUrbxKjNIXrXtb77mD6qZR4rYeKRPFzcXyAHAHqYGc9G09jTnFuG/eFH1ET9H0ET1F01eJ0ZjT8Bqn6aRUR9uy9oGIHcbRn8YAx5yUvihGCpYm6SN6iqaP6CmavmyM0hcVDGOek5ntR7p1Nw041d3fZ2ZLANz9413a0+gz5tSJ5jkJIdpCW+Y5KX2R0hcJIWpGG9IXtbJxEkIIEbtxmvIMERFYefdKlp6/lNFjRxlZNsLosaMsPX9pM2ZVCyEaTxvOYa3rOWnMSQhRZ9oy5tSqnpNy68XID9YkfURP0fQRPUXTl41Rbr2aLZrnlF8f0ZPqnF8f0VM0fZUYzXOq2aLcevn1ET2pzvn1ET1F01eJUW69mlF2zEm59Yanj+hJdc6vj+gpmr5KjHLrNRTl1hteeURPqnP+8oieopVXiVFuvZotGnPKr4/oSXXOr4/oKZq+SozGnGq2KLdeDH1ET6pzfn1ET9H0ZWOUW69maJ6TEKItaJ5TQ9l3u31ZvmQ5i+cvZnSjUUZshNGNRlk8fzHLlyxXwySECE1bzmGt6zkJIYRIqOcUjDbkpRJCNJc2nMNa13PSmJMQos5ozKmBKLdejPxgTdJH9BRNH9FTNH3ZGOXWq9mieU759RE9qc759RE9RdNXidE8p5otyq2XXx/Rk+qcXx/RUzR9lRjl1qsZyq2XXx/Rk+qcXx/RUzR9lRjl1msoyq03vPKInlTn/OURPUUrrxKj3Ho1WzTmlF8f0ZPqnF8f0VM0fZUYjTnVbFFuvRj6iJ5U5/z6iJ6i6cvGKLdezdA8JyFEW9A8p4bSlrxUQohm0pZzWOt6TkIIIRLqOQWjDXmphBDNpQ3nsNb1nDTmJISoMxpzaiDKrRcjP1iT9BE9RdNH9BRNXzZGufVqtmieU359RE+qc359RE/R9FViNM+pZoty6+XXR/SkOufXR/QUTV8lRrn1aoZy6+XXR/SkOufXR/QUTV8lRrn1Gopy6w2vPKIn1Tl/eURP0cqrxCi3Xs0WjTnl10f0pDrn10f0FE1fJUZjTjVblFsvhj6iJ9U5vz6ip2j6sjHKrVczNM9JCNEWNM+pobQlL5UQopm05RzWup6TEEKIhHpOwWhDXiohRHNpwzmsdT0njTkJIeqMxpwaiHLrxcgP1iR9RE/R9BE9RdOXjVFuvZotmueUXx/Rk+qcXx/RUzR9lRjNc6rZotx6+fURPanO+fURPUXTV4lRbr2aodx6+fURPanO+fURPUXTV4lRbr2Gotx6wyuP6El1zl8e0VO08ioxyq1Xs0VjTvn1ET2pzvn1ET1F01eJ0ZhTzRbl1ouhj+hJdc6vj+gpmr5sjHLr1QzNcxJCtAXNc2oobclLJYRoJm05h7Wu5ySEECKhnlMw2pCXSgjRXNpwDpvynpOZ7QOcAEwDTnH393etfy3wzuLlA8Cb3P3a8bapMSchRFvQmNMQMLNpwEeBfYHtgYPMbPsu2U3AXu6+E/DvwMmDen/l1ouRH6xJ+oieoukjeoqmLxuj3HrDem4d9gAu7Hh9JHDkOPrNgV9NtF3Nc8qvj+hJdc6vj+gpmr5KjOY5Da9xOoB0K2/s9SLgI+Po/6VT329Rbr38+oieVOf8+oieoumrxCi33pAwswOBvd39jcXrRcBz3f2femhfAJwE7Onud/VYvxhYDDBnzpz5q1atmvD9lVtvePqInlTn/PqInqLpq8Qot97wuBV4SsfrJwO3dYvMbCfgFODlvRomAHc/2d0XuPuC2bNnl3pz5dYbXnlET6pz/vKInqKVV4lRbr0hLcB04EZgW2AmcC2wQ5dmDrACeF7Z7WrMKb8+oifVOb8+oqdo+ioxGnMabgO1H3ADsBL4t6JsCbCk+P8U4B7gmmKZ8MNTbr0Y+oieVOf8+oieounLxii3Xs3QPCchRFvQPKeG0pa8VEKIZtKWc1jrek5CCCES6jkFow15qYQQzaUN57DW9Zw05iSEqDMac2ogyq0XIz9Yk/QRPUXTR/QUTV82Rrn1arZonlN+fURPqnN+fURP0fRVYjTPqWaLcuvl10f0pDrn10f0FE1fJUa59WpG2TEn5dYbnj6iJ9U5vz6ip2j6KjHKrddQlFtveOURPanO+csjeopWXiVGufVqtmjMKb8+oifVOb8+oqdo+ioxGnOq2aLcejH0ET2pzvn1ET1F05eNUW69mqF5TkKItqB5Tg2lLXmphBDNpC3nsNY1TmOs8TXrdCHXeO+naIQQIiJNP4fptl6BbusJIeqAbus1EKUvipGCpUn6iJ6i6SN6iqYvG6P0RTVb9Ch5fn1ET6pzfn1ET9H0VWL0KHnNFqUvyq+P6El1zq+P6CmavkqM0hfVDKUvyq+P6El1zq+P6CmavkqM0hc1FKUvGl55RE+qc/7yiJ6ilVeJUfqimi0ac8qvj+hJdc6vj+gpmr5KjMacarYofVEMfURPqnN+fURP0fRlY5S+qGZonpMQoi1onlNDaUvqDyFEM2nLOax1PSchhBAJ9ZyCsfLulSw9fymjx44ysmyE0WNHWXr+0mbMqhZCNJ42nMNa13PSmJMQos5ozKmBKLdejPxgTdJH9BRNH9FTNH3ZGOXWq9mieU759RE9qc759RE9RdNXidE8p5otyq2XXx/Rk+qcXx/RUzR9lRjl1qsZyq2XXx/Rk+qcXx/RUzR9lRjl1msoyq03vPKInlTn/OURPUUrrxKj3Ho1WzTmlF8f0ZPqnF8f0VM0fZUYjTnVbFFuvRj6iJ5U5/z6iJ6i6cvGKLdezdA8JyFEW9A8p4bSlrxUQohm0pZzWOt6TkIIIRLqOQWjDXmphBDNpQ3nsNb1nDTmJISoMxpzaiDKrRcjP1iT9BE9RdNH9BRNXzZGufVqtmieU359RE+qc359RE/R9FViNM+pZoty6+XXR/SkOufXR/QUTV8lRrn1aoZy6+XXR/SkOufXR/QUTV8lRrn1Gopy6w2vPKIn1Tl/eURP0cqrxCi3Xs0WjTnl10f0pDrn10f0FE1fJUZjTjVblFsvhj6iJ9U5vz6ip2j6sjHKrVczNM9JCNEWNM+pobQlL5UQopm05RzWup6TEEKIhHpOwbj4povZ8aQdsWX26LLjSTty8U0X942pmsuqbfqInlTn/PqInqLpNyRmQ85hdaN1Paf3fPs9HH3J0X3XL1u4jKP2Omqdsqr3eNumj+hJdc6vj+gpmn5DYjbkHNaPyD2nKW+czGwf4ARgGnCKu7+/a70V6/cDHgRe7+5Xj7fNso3TxTddzItOf9GEuosOuYgXbvtCIF3R7PTxnXhw9YN99bNmzGL5kuXM22Je6/TQvs+ojXXWZzR4PVT/jDbkHDYekRunKb2tZ2bTgI8C+wLbAweZ2fZdsn2B7YplMfCxQb3/m7/25lK6f/76Pz/6/we//0FWP7J6HDWsfmQ1x11+XCv1ET2pzvn1ET1F029IzIacw+rKlPaczGwP4Bh337t4fSSAux/bofkEcIm7f7Z4fT2w0N1v77fdsj0nW2alvfrR6XMZPXaU+x+6f0L96Eaj3HfEfa3TQ/s+I2hfnfUZDV4P1T+jDTmHjYd6TmvZGril4/WtRVlVDWa22MyuNLMrf/vb3w7c6BgPPPRAJV3b9BE9qc759RE9RdNvaExbmOrGqVez3928l9Hg7ie7+wJ3XzB79uyBmOtF1VxWbdNH9KQ659dH9BRNv6ExbWGqG6dbgad0vH4ycNsGaDaIHWbvUEq34+N3fPT/g3c6mBkjM8bVzxiZwaKdFrVSH9GT6pxfH9FTNP2GxGzIOayuTHXj9ENgOzPb1sxmAn8HnNelOQ84xBK7A/eNN95UhRP3PbGU7oR9Tnj0/8P3OJwZ0ybYeabN4K27v7WV+oieVOf8+oieouk3JGZDzmF1ZUobJ3d/GDgMuBC4Djjb3X9qZkvMbEkhuwC4EVgBfBJYOqj3f+G2L2TZwmXjapYtXLbOI5jztpjHOQeew6wZs9a7wpkxMoNZM2ZxzoHnPPpoaNv0ET2pzvn1ET1F029IzIacw2pL7syzg1jKZiUf46IbL/IdT9pxnUy+O560o19040V9Y1bctcIPPf9QHz121EeWjfjosaN+6PmH9s3+2zZ9RE+qc359RE/R9BsSsyHnsF6grOTDRbn1hBCiOnqUXAghhKiAGichhBDhUOMkhBAiHGqchBBChEONkxBCiHCocRJCCBEONU5CCCHCocZJCCFEOBoxCdfMfgus2sDwrYA7B2inDqjO7UB1bgeTqfM27j68n3WYBI1onCaDmV0ZdYb0sFCd24Hq3A6aWmfd1hNCCBEONU5CCCHCocYJTs5tIAOqcztQndtBI+vc+jEnIYQQ8VDPSQghRDha0ziZ2T5mdr2ZrTCzI3qsNzM7sVi/3Mx2yeFzkJSo82uLui43s++Z2c45fA6SiercodvVzB4xswOm0t8wKFNnM1toZteY2U/N7NtT7XHQlNi3H2dmXzGza4s6vyGHz0FhZqea2R1m9pM+6xt3/sr+a4dTsQDTgJXAU4GZwLXA9l2a/YCvAQbsDvwgt+8pqPPzgM2L//dtQ507dBcDFwAH5PY9Bd/zZsDPgDnF68fn9j0Fdf5X4D+L/2cDdwMzc3ufRJ2fD+wC/KTP+kadv9y9NT2n5wIr3P1Gd38I+Bzw8i7Ny4HTPXE5sJmZPXGqjQ6QCevs7t9z93uKl5cDT55ij4OmzPcM8E/AF4A7ptLckChT59cA57r7LwHcve71LlNnBzY1MwM2ITVOD0+tzcHh7peS6tCPpp2/WtM4bQ3c0vH61qKsqqZOVK3P/yNdedWZCetsZlsDrwA+PoW+hkmZ7/npwOZmdomZXWVmh0yZu+FQps4fAZ4F3Ab8GPhnd18zNfay0LTzF9NzG5girEdZ92OKZTR1onR9zOwFpMZpz6E6Gj5l6nw88E53fyRdVNeeMnWeDswHXgQ8Bvi+mV3u7jcM29yQKFPnvYFrgBcC84Bvmtl33P13wzaXiaadv1rTON0KPKXj9ZNJV1RVNXWiVH3MbCfgFGBfd79rirwNizJ1XgB8rmiYtgL2M7OH3f1LU2Nx4JTdt+90998DvzezS4Gdgbo2TmXq/Abg/Z4GZFaY2U3AM4ErpsbilNO081drbuv9ENjOzLY1s5nA3wHndWnOAw4pnnrZHbjP3W+faqMDZMI6m9kc4FxgUY2vojuZsM7uvq27z3X3ucA5wNIaN0xQbt/+MvCXZjbdzGYBuwHXTbHPQVKmzr8k9RQxsycAzwBunFKXU0vTzl/t6Dm5+8NmdhhwIelJn1Pd/admtqRY/3HSk1v7ASuAB0lXXrWlZJ2PArYETip6Eg97jRNIlqxzoyhTZ3e/zsy+DiwH1gCnuHvPR5LrQMnv+d+B08zsx6RbXu9099pmKzezzwILga3M7FbgaGAGNPP8BcoQIYQQIiBtua0nhBCiRqhxEkIIEQ41TkIIIcKhxkkIIUQ41DgJIYQIhxon0QjMbH8zu7TI3PwHM1tlZl8ys31ye5ssZvZhM/vKALd3gpmdP6jtCTEM9Ci5qD1m9mbgBOBU4EvA70kpa14G3ODu78hob1KY2TzShNnnufuVA9rmE0kTUl/m7hcPYptCDBo1TqL2mNkvgavc/RU91o3UOeGnmX0Y2N3ddx3wds8GNnH3/Qa5XSEGhW7riSawBfDrXis6GyYze72ZuZk9v7jl94CZ3WVmHzWzx3TGmdkyM7vazO4zszvN7OIiLQxdutlmdpKZ3WJmfyr+nmFmG3Vodjaz88zsnuKW43fN7C8nqlSxjYOBz3SVLyzqsb+ZfcLM7i62fZyZTbP0Q4qXmdnvLf3Q3t49Nv85YG8ze0qPdUJkR42TaAJXAK8zs7eb2dNL6M8kpXl5JXAc8A/Ax7o0Wxfr9gdeT/rtp0uLRLkAmNnmwPeAVwMfIqWPeQcprczMQrNLodmieJ9XAXcB/2dm8yfwuTvphwK/02f98aRbmK8m/UTEW4qy00m3OF9J+g2gc81sq67YS0nH/0sm8CBEHnL/2qEWLZNdSL9XtJz0EwEO3Al8Fnhpl+71xfqPd5X/G/AI8PQ+259GykN5PXBCR/l7irjnjOPtItKY0cyu7V0HfGmCer2TlAtvZlf5wqIep3aVX12U79lRtlNR9roe278FODn396dFS69FPSdRezxlVH8OsBfwPtLv+LwCuNDM3tUj5Oyu158j9SKeO1ZgZi82s2+Z2V2kX1BdTWoEn9ER91Lgh+7+o16+iluFewH/C6wpsoJPJyUi/T/ST2+Px5OA33n6tddedP845M+B37v7ZV1lsO7PKYzx2+I9hAiHGifRCNz9EXe/1N3f5e4vBp5K+gXUo4vbb538ps/rreHRW3EXAA+QfoRxd2BX4Fpg4464LUm/o9OPLUi9pHeTGrfO5TDSr9OOdwxuDPxpnPX3dL1+CLi3s6CjYduY9fkD6ccHhQhHK34yQ7QPd7/NzE4hPWK+Hev+yNwTgJ92vQb4VfH3VaTe0ivdffWYqGjkOk/+dzL+T2HfS7ot91HSOFAvn+M9SXgX0N2wDpItSLdDhQiHek6i9ozzxNkzi7/dT/L9bdfrvyM1ImMN2CzSWNKj8yzM7IXAnK64bwDPNbOde725p1+e/Q7pV2evdvcru5dxqgXpltwMM3vyBLrKmNk00q2+6we9bSEGgXpOogn8xMy+BXwRuAkYJT05twQ4291/2aXfz8z+m6JxIf1w2+m+9teAv0568u00M/s0aazp3aztWY1xHPAa0pN37yXdRtwKeDmwxN3vB95GejLuQjP7FHB7odkFmObuR4xTr0uLv89l/NuHG8KOwGM73kOIUKjnJJrAO0n78ntIDc7ngT2AI4BFPfQHkxqcLwKHA58Elo6tdPcLgTcDfwF8Ffh74BDS4+d06O4tNF8s3uvrwAdJtwQfKjRXk8ar7gJOLPydAPw5EzQM7n4zqTf31xN/BJX5K1KP8pIhbFuISaMMEaI1mNnrgU8D27n7ignkISg8nwA80d0fHOB2fwZ8wd3fPahtCjFI1HMSIjZnkG4nLp1IWBYzeznpIZAPDmqbQgwaNU5CBMbdHyHdVhxYr4n0+PjBxW1JIUKi23pCCCHCoZ6TEEKIcKhxEkIIEQ41TkIIIcKhxkkIIUQ41DgJIYQIhxonIYQQ4fj/R1ClCjJZbLEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "nx=nt=26\n",
    "x = np.arange(nx)/(nx-1)\n",
    "t=np.zeros(nt)\n",
    "\n",
    "for it in range(nt):\n",
    "    plt.plot(x,t,'ko')\n",
    "    t+=1/(nt-1)\n",
    "plt.plot(x,np.zeros(nx),'go',ms=10)\n",
    "plt.plot(np.zeros(nt),x,'go',ms=10)\n",
    "plt.plot(np.zeros(nx)+1,x,'go',ms=10)    \n",
    "plt.ylabel('Time (s)',fontsize=16)\n",
    "plt.xlabel('Space (m)',fontsize=16)\n",
    "plt.title('Solution space for 1D Diffusion Problem',fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 1. Example solution space for 1D diffusion problem.  The green points are the locations where the solution is fixed by the initial condition (bottom) or the boundary conditions (left and right sides).  The black dots are the locations where we need to compute the solution.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explicit FTCS solution\n",
    "\n",
    "The first numerical solution that we will investigate is an explicit method: forward in time, centered in space (FTCS). This method uses an $\\mathcal{O}(\\Delta t)$ forward approximation for the temporal first derivative and an $\\mathcal{O}(\\Delta x^2)$ centered approximation for the spatial second derivative. Thus, this numerical approximation can be written as:\n",
    "\n",
    "$$\\frac{\\phi^{n+1}_i -\\phi^{n}_i}{\\Delta t} = \\kappa \\frac{\\left(\\phi^{n}_{i+1}-2\\phi^{n}_i+\\phi^{n}_{i-1}  \\right)}{\\Delta x^2}.\\tag{5}$$\n",
    "\n",
    "Let's now multiply both sides by $\\Delta t$ and then rearrange equation 5 so that we have all terms at the $n+1$ time step on the left-hand side, and all others on the right-hand side:\n",
    "\n",
    "$$\\phi^{n+1}_i = \\alpha \\phi^{n}_{i+1}+(1-2\\alpha)\\phi^{n}_i+\\alpha\\phi^{n}_{i-1}.\\tag{6}$$\n",
    "\n",
    "where $\\alpha = \\frac{\\kappa \\Delta t}{\\Delta x^2}$ is introduced as a dimensionless number called the **diffusion number**.\n",
    "\n",
    "Let's now code up a solver for this system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def diffusion_ftcs_solver_homogeneous_Dirchelet(U,dx,dt,K,LB,RB):\n",
    "    '''Set up FTCS Solver of the Diffusion equations\n",
    "    usage: U = ADI_Solution(U,dx,dt,K,LB,RB)\n",
    "    input: \n",
    "        U : Heat distribution vector (nx) at time t\n",
    "        dx: Spatial sampling\n",
    "        dt: Temporal sampling\n",
    "        K : Thermal conductivity\n",
    "        LB: Left  boundary (Homogeneous Dirchelet)\n",
    "        RB: Right boundary (Homogeneous Dirchelet)\n",
    "    output:\n",
    "        U: Heat distribution vector (nx) at time t+∆t\n",
    "    dependencies:\n",
    "        Setup_Tridiagonal\n",
    "        np.linalg.solve\n",
    "    written by Jeff Shragge, jshragge@mines.edu, 10/2019\n",
    "    '''        \n",
    "    ## . . Get dimensions\n",
    "    nx = len(U)\n",
    "    \n",
    "    ## . . Diffusion number\n",
    "    alpha = K*dt/(dx*dx)\n",
    "    \n",
    "    ## . . Update solution vector\n",
    "    U[1:nx-1] =      alpha *U[2:nx  ]+ \\\n",
    "                (1-2*alpha)*U[1:nx-1]+ \\\n",
    "                     alpha *U[0:nx-2]\n",
    "            \n",
    "    ## . . Enforce boundary conditions\n",
    "    U[0   ] = LB\n",
    "    U[nx-1] = RB\n",
    "    \n",
    "    return U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's also define the *diffusion_1D_animate(i)* animate function once so that we can call it frequently below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def diffusion_1D_animate(i):\n",
    "    global k\n",
    "    ax1.clear()\n",
    "    ax1.plot(xx,c[:,k],'g',xx,d[:,k],'b',linewidth=3)\n",
    "    ax1.legend(['Copper','Iron'])\n",
    "    ax1.grid(True)\n",
    "    ax1.set_xlim([xmin,xmax])\n",
    "    ax1.set_ylim([-0.1*T0,1.1*T0])\n",
    "    ax1.set(xlabel='Distance (m)',ylabel='Temperature ($^\\circ$C)')\n",
    "    ax1.set_title('Heat Flow at %f s'%(k*dt))\n",
    "    k += kskip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Copper and Iron Bars with Fixed Boundaries\n",
    "\n",
    "Let's now create an example that calls our diffusion solver. We're going to discretize a solution domain $x\\in[0,1]$m with $nx=51$ points ($dx=0.02$m) and $t\\in[0,1]$s also with $nt=101$ points ($dt=0.01$s). As the initial distribution, let's use a Gaussian function\n",
    "\n",
    "$$\\phi_0(x) = 100 {\\rm e}^{-(x-x_0)^2/(2\\sigma^2)} \\tag{7}$$\n",
    "\n",
    "where $x_0=0.5$m, $\\sigma=0.05$m and the units of $\\phi$ are in $^\\circ$C. For the thermal diffusivity, we are going to look at two materials:\n",
    "   * Copper: $\\kappa = 1.1 \\times 10^{-4}$ m$^2/s$ (at $25^\\circ$C)\n",
    "   * Iron  : $\\kappa = 2.3 \\times 10^{-5}$ m$^2/s$\n",
    "\n",
    "You'll note that the $\\kappa$ values are defined at a constant temperature. Generally speaking, the thermal diffusivity of materials can depend on the temperature itself, which will make the problem **non-linear** and more difficult to solve.\n",
    "\n",
    "In this scenario, you might recall that the steady-state solution would be for the bars to become uniformly $0^\\circ$C. Let's see if our numerical solution accomplishes that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha (Copper) = 0.4125\n",
      "alpha (Iron  ) = 0.08625\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_1.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=1001,1.5\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 51\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/2,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kc       = 1.1*10**(-4) ## . . thermal diffusivity\n",
    "Ki       = 2.3*10**(-5) ## . . thermal diffusivity\n",
    "LB = RB = 0 ## . . Boundary conditions on left and right\n",
    "\n",
    "print('alpha (Copper) =',Kc*dt/(dx*dx))\n",
    "print('alpha (Iron  ) =',Ki*dt/(dx*dx))\n",
    "\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "c[:,0]=Uc\n",
    "d[:,0]=Ui\n",
    "## . . Iterate over solution\n",
    "for i in range(1,nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = diffusion_ftcs_solver_homogeneous_Dirchelet(Uc,dx,dt,Kc,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = diffusion_ftcs_solver_homogeneous_Dirchelet(Ui,dx,dt,Ki,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "    \n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig1,ax1 = plt.subplots(1)\n",
    "fig1.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig1.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim1 = animation.FuncAnimation(fig1,diffusion_1D_animate,frames=int((nt-kskip)/kskip),interval=100)\n",
    "anim1.save('./movies/Ex7_1.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_1.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 2. Solutions of the 1D heat diffusion equation for an insulated 1m metal bar with the exposed ends held at 0ºC (i.e., zero Dirchelet boundary conditions). Solutions for copper (green) and iron (blue) are presented.**\n",
    "\n",
    "### Playing with Boundary Conditions - Non-zero Dirichelet\n",
    "\n",
    "Let's now look at the same example as above, but let's change the boundary conditions such that $\\phi(x=0,t)=0^\\circ$C and $\\phi(x=1,t)=50^\\circ$C.  Let's rerun the code from above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha (Copper) = 0.4125\n",
      "alpha (Iron  ) = 0.08625\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_2.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=1201,1.5\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 51\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/2,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kc       = 1.1*10**(-4) ## . . thermal diffusivity\n",
    "Ki       = 2.3*10**(-5) ## . . thermal diffusivity\n",
    "LB,RB = 0,50 ## . . Boundary conditions on left and right\n",
    "\n",
    "print('alpha (Copper) =',Kc*dt/(dx*dx))\n",
    "print('alpha (Iron  ) =',Ki*dt/(dx*dx))\n",
    "## . . Iterate over solution\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "c[:,0]=Uc\n",
    "d[:,0]=Ui\n",
    "for i in range(1,nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = diffusion_ftcs_solver_homogeneous_Dirchelet(Uc,dx,dt,Kc,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = diffusion_ftcs_solver_homogeneous_Dirchelet(Ui,dx,dt,Ki,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "    \n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig2,ax1 = plt.subplots(1)\n",
    "fig2.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig2.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim2 = animation.FuncAnimation(fig2,diffusion_1D_animate,frames=int((nt-kskip)/kskip),interval=100)\n",
    "anim2.save('./movies/Ex7_2.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_2.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 3. Solutions of the 1D heat diffusion equation for an insulated 1m metal bar with the exposed ends held at 0ºC on the left and 50ºC on the right. (i.e., Dirchelet boundary conditions). Solutions for copper (green) and iron (blue) are presented.**\n",
    "\n",
    "### Playing with boundary conditions - Neumann boundary conditions\n",
    "\n",
    "The examples above looked at Dirichelet boundary conditions, which is where the boundaries are fixed at a certain value. However, there are important scenarios where you might want to control the heat gradient (i.e., $\\frac{\\partial \\phi}{\\partial x}$) at the boundary.\n",
    "\n",
    "Let's say that we want to modify the above such that we had:\n",
    "\n",
    "$$\\frac{\\partial \\phi(x=0,t)}{\\partial x}=A \\tag{8}$$\n",
    "\n",
    "We can think about applying this by discretizing equation 8 as:\n",
    "\n",
    "$$\\frac{U^n_1-U^n_0}{\\Delta x} = A \\tag{9a}$$\n",
    "\n",
    "or solving for $U^n_0$:\n",
    "\n",
    "$$U^n_0= U^n_1-A\\Delta x. \\tag{9b}$$\n",
    "\n",
    "Similarly, at the right boundary we have\n",
    "\n",
    "$$\\frac{\\partial \\phi(x=L,t)}{\\partial x}=B \\tag{10}$$\n",
    "\n",
    "and thus numerically we have:\n",
    "\n",
    "$$\\frac{U^n_N-U^n_{N-1}}{\\Delta x} = B \\tag{11a}$$\n",
    "\n",
    "or solving for $U^n_N$:\n",
    "\n",
    "$$U^n_N= U^n_{N-1}+B\\Delta x. \\tag{11b}$$\n",
    "\n",
    "where the difference in sign in equations 9b and 11b is due to the opposing directions of normal derivative into the bar.  Let's redo our solver above but use Neumann BC's:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def diffusion_ftcs_solver_homogeneous_Neumann(U,dx,dt,K,LB,RB):\n",
    "    '''Set up FTCS Solver of the Diffusion equations\n",
    "    usage: U = ADI_Solution(U,dx,dt,K,LB,RB)\n",
    "    input: \n",
    "        U : Heat distribution vector (nx) at time t\n",
    "        dx: Spatial sampling\n",
    "        dt: Temporal sampling\n",
    "        K : Thermal conductivity\n",
    "        LB: Left  boundary (Homogeneous Neuman)\n",
    "        RB: Right boundary (Homogeneous Neuman)\n",
    "    output:\n",
    "        U: Heat distribution vector (nx) at time t+∆t\n",
    "    dependencies:\n",
    "        Setup_Tridiagonal\n",
    "        np.linalg.solve\n",
    "    written by Jeff Shragge, jshragge@mines.edu, 10/2019\n",
    "    '''        \n",
    "    ## . . Get dimensions\n",
    "    nx = len(U)\n",
    "    \n",
    "    ## . . Diffusion number\n",
    "    alpha = K*dt/(dx*dx)\n",
    "    \n",
    "    ## . . Update solution vector\n",
    "    U[1:nx-1] =      alpha *U[2:nx  ]+ \\\n",
    "                (1-2*alpha)*U[1:nx-1]+ \\\n",
    "                     alpha *U[0:nx-2]\n",
    "            \n",
    "    ## . . Enforce boundary conditions\n",
    "    U[0   ] = U[1   ]-LB*dx\n",
    "    U[nx-1] = U[nx-2]+RB*dx\n",
    "    \n",
    "    return U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now rerun the example above but with Neumann boundary conditions of $\\frac{\\partial \\phi(x=0,t)}{\\partial x}=1/\\Delta x^\\circ$C and $\\frac{\\partial \\phi(x=1,t)}{\\partial x}=0^\\circ$C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha (Copper) = 0.4125\n",
      "alpha (Iron  ) = 0.08625\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_3.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=400,1.5\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 51\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/2,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kc       = 1.1*10**(-4) ## . . thermal diffusivity\n",
    "Ki       = 2.3*10**(-5) ## . . thermal diffusivity\n",
    "LB,RB = -1/dx,0 ## . . Derivative Boundary conditions on left and right\n",
    "\n",
    "print('alpha (Copper) =',Kc*dt/(dx*dx))\n",
    "print('alpha (Iron  ) =',Ki*dt/(dx*dx))\n",
    "## . . Iterate over solution\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "for i in range(nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = diffusion_ftcs_solver_homogeneous_Neumann(Uc,dx,dt,Kc,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = diffusion_ftcs_solver_homogeneous_Neumann(Ui,dx,dt,Ki,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "    \n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig3,ax1 = plt.subplots(1)\n",
    "fig3.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig3.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim3 = animation.FuncAnimation(fig3,diffusion_1D_animate,frames=int((nt-kskip)/kskip),interval=100)\n",
    "anim3.save('./movies/Ex7_3.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_3.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 3. Solutions of the 1D heat diffusion equation for an insulated 1m metal bar with the exposed ends experiencing a constant inward flux of (1/∆x)ºC/m on the left and 0ºC/m on the right. (i.e., Dirchelet boundary conditions). Solutions for copper (green) and iron (blue) are presented.**\n",
    "\n",
    "## Stability of the FTCS method\n",
    "\n",
    "One important consideration that have not discussed is whether or not the FTCS method is  **unconditionally stable** and, if not, under what conditions is it **conditionally stable**. You may have noticed that I have written out the values of $\\alpha$ when running each example, and that the maximum value has been $\\alpha=0.4125$ m$^2$/s.  Let's now redo the very first example but increase $dt=2$ such that maximum value (for copper) is now $\\alpha=0.55$ m$^2$/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha (Copper) = 0.605\n",
      "alpha (Iron  ) = 0.1265\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_4.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=101,2.2\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 51\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/2,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kc       = 1.1*10**(-4) ## . . thermal diffusivity\n",
    "Ki       = 2.3*10**(-5) ## . . thermal diffusivity\n",
    "LB = RB = 0 ## . . Boundary conditions on left and right\n",
    "\n",
    "print('alpha (Copper) =',Kc*dt/(dx*dx))\n",
    "print('alpha (Iron  ) =',Ki*dt/(dx*dx))\n",
    "\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "## . . Iterate over solution\n",
    "for i in range(nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = diffusion_ftcs_solver_homogeneous_Dirchelet(Uc,dx,dt,Kc,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = diffusion_ftcs_solver_homogeneous_Dirchelet(Ui,dx,dt,Ki,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "    \n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig4,ax1 = plt.subplots(1)\n",
    "fig4.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig4.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim4 = animation.FuncAnimation(fig4,diffusion_1D_animate,frames=int((nt-kskip)/kskip),interval=100)\n",
    "anim4.save('./movies/Ex7_4.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_4.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 4. Solutions of the 1D heat diffusion equation as in Figure 3, except that the value of $\\alpha$ for the copper bar has been increased to 0.605 instead of 0.4125.  Note that the copper bar simulation rapidly becomes unstable, where as the iron bar simulation remains stable ($\\alpha=0.1265 m^2/s$).**\n",
    "\n",
    "So it appears that we have crossed some sort of threshold between when $\\alpha=0.4125$ m$^2$/s and $\\alpha=0.55$ m$^2$/s. To examine this, let's again look at the von Neumann stability analysis that we discussed in Module 7.\n",
    "\n",
    "## Von Neumann Stability Analysis\n",
    "\n",
    "We want to investigate how the error difference $\\epsilon^{n}$ between the true solution to the 1D diffusion PDE, $\\phi_T$ and our numerical solution, $\\phi$, changes from time step $n$ to $n+1$. Recall that this term is called the **amplification factor**\n",
    "\n",
    "$$g^2 = \\overline{g}g=\\left|\\frac{\\epsilon^{n+1}_i}{\\epsilon^n_i}\\right|^2 \\le 1, \\tag{12}$$\n",
    "\n",
    "Again, we can assume that the error term satisfies the discretization itself the amplification factor is given by:\n",
    "\n",
    "$$g^2 = \\left|\\frac{\\epsilon^{n+1}_i}{\\epsilon^n_i}\\right|^2 = \\left|\\frac{\\alpha \\epsilon^n_{i+1}+\\left(1-2\\alpha\\right)\\epsilon^n_i+\\alpha \\epsilon^n_{i-1} }{\\epsilon^n_i}\\right|^2 \\le 1  \\tag{14}$$\n",
    "\n",
    "We now can assume that each term can be represented by \n",
    "\n",
    "$$\\epsilon^n = A_n e^{ik x}, \\tag{15}$$\n",
    "\n",
    "where $A^n$ is an amplitude term at time step $n$ such that \n",
    "\n",
    "$$\\epsilon^n_i     = A_n e^{ik x},            \\tag{16a}$$\n",
    "$$\\epsilon^n_{i+1} = A_n e^{ik (x+\\Delta x)}, \\tag{16b}$$\n",
    "$$\\epsilon^n_{i-1} = A_n e^{ik (x-\\Delta x)}. \\tag{16c}$$\n",
    "\n",
    "Inserting these into equation 14 yields\n",
    "\n",
    "$$g^2 = \n",
    "\\left|\\frac{\\alpha A_n e^{ik (x+\\Delta x)}+\n",
    "\\left(1-2\\alpha\\right)A_n e^{ik x}+\n",
    "\\alpha A_n e^{ik (x-\\Delta x)}}\n",
    "{A_n e^{ik x}}  \\right|^2 \\le 1\n",
    "\\tag{17}$$\n",
    "\n",
    "which simplifies to \n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "g^2 &=& \n",
    "\\left|\\alpha e^{ik\\Delta x}+\n",
    "\\left(1-2\\alpha\\right)+\n",
    "\\alpha e^{-ik\\Delta x} \\right|^2\\le 1\\\\\n",
    "&=& \\left|1+2\\alpha\\,(\\cos\\left(k\\Delta x\\right)-1)\\right|^2 \\le 1\\\\\n",
    "&=& \\left|1-4\\alpha \\sin^2\\left(\\frac{k\\Delta x}{2}\\right)\\right|^2 \\le 1\n",
    "\\end{eqnarray}\n",
    "\\tag{18}$$\n",
    "\n",
    "where in the last step we have used the half-angle formula $2\\sin^2 \\theta =1-\\cos 2\\theta$.  Thus, the **stability criterion** for this numerical scheme is given by\n",
    "\n",
    "$$g^2=\\left|1-4\\alpha \\sin^2\\left(\\frac{k\\Delta x}{2}\\right)\\right|^2\\le 1. \\tag{19}$$\n",
    "\n",
    "Note that the $\\sin^2$ factor ranges between 0 and 1; when it is 0 this criterion is always satisfied.  Thus, we are interested in the maximal case when $\\sin^2=1$.  Accordingly, we rewrite the inequality in equation 18 as\n",
    "\n",
    "$$g^2 = \\left|1-4\\alpha\\right|^2 \\le 1, \\tag{20}$$\n",
    "\n",
    "which is satisfied in the range when $\\alpha\\in[0,1/2]$.  Thus, we have the following restriction on the numerical simulation:\n",
    "\n",
    "$$\\alpha = \\frac{\\kappa \\Delta t}{\\Delta x^2} < \\frac{1}{2} \\tag{21}$$\n",
    "\n",
    "or for a more explicit constraint on the allowed time step $\\Delta t$:\n",
    "\n",
    "$$\\Delta t < \\frac{\\Delta x^2}{2\\kappa}, \\tag{22}$$\n",
    "\n",
    "which is consistent with the numerical simulations above. Thus, while the above FTCS approach is relatively straightforward to implement, it suffers from the severe drawback that it is only **conditionally stable**. This motivates use to look for other numerical methods that are **unconditionally stable**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Implicit Crank-Nicholson method\n",
    "\n",
    "Fortunately, we do not have to look too far to find such a method. The Crank-Nicholson (CN) method represents an extension of the FTCS approach that is $\\mathcal{O}\\left(\\Delta t^2,\\Delta x^2\\right)$ and **unconditionally stable**.  The main idea behind this approach is that unlike the FTCS method that defines a forward-difference formula for **one unknown point** at time step $n+1$, the CN method specifies it for **three unknown points** at time step $n+1$ in the FD stencil. Figure 9-5 below illustrates the stencil associated with this numerical scheme."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAGBCAYAAABmcMZgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3X+YXAV97/H3dwOJGWtvRQLeojtDECVyETBbkdarROqPem/FWpoWYyQ8KCCPvdeqrVr8gcpa7+21Xm+fqk1VVHbVBxV/9QdaBdFSQYNggj+rkFC0Kj9FDSaQfO8f5yxsNrvZc8jMnJnd9+t5zrPMmTPnfOdL5nzm/JhzIjORJGmk6QIkSYPBQJAkAQaCJKlkIEiSAANBklQyECRJgIEgSSoZCJIkwECQJJUMBEkSAAc0XUAdBx98cHY6nabLkKShcs0119yamSvmm26oAqHT6bBp06amy5CkoRIR26pM5y4jSRJgIEiSSgaCJAkwECRJJQNBkgQYCJKkkoEgSQIGIBAi4tKIyIi4oOlaJGkxazQQIuI04Ngma5AkFRoLhIj4NeBtwMt6vazJyUk6nQ4jIyN0Oh0mJyd7vcihZr/qsV/12K96+tqvzGxkADYCny//O4EL5nvN6tWrs66JiYlstVpZLiOBbLVaOTExUXtei4H9qsd+1WO/6ulWv4BNWWW9XGWibg/Ak4BfAo/JHgdCu93eo5lTQ7vdrj2vxcB+1WO/6rFf9XSrX1UDIYpp+yciDgSuAz6ema8pxyUwPvV4xvRnAWcBjI6Ort62rdI1mu4zMjLCbO8xIti9e3f9N7DA2a967Fc99quebvUrIq7JzLF5l1evvK54JbAcGK8ycWZuzMyxzBxbsWLeq7fuZXR0tNb4xc5+1WO/6rFf9fS7X30NhIgYBc4DXgssi4hfKw8uM+3xkm4uc3x8nFartce4VqvF+HilPFp07Fc99qse+1VP3/tVZb9StwbgJGbZHzZjOG6u1z+QYwiZxYGZdrudEZHtdtsDWPOwX/XYr3rsVz3d6BeDeAyh3Bo4bpanLgcmgPeUhf98ttePjY2lN8iRpHqqHkPo6x3TMvNO4Aszx0cEwLbM3Os5SVJ/NH7pCknSYBiIeypnZjRdgyQtdm4hSJIAA0GSVDIQJEmAgSBJKhkIkiTAQJAklQwESRJgIEiSSgaCJAkwECRJJQNBkgQYCJKkkoEgSQIMBElSyUCQJAEGgiSpZCBIkgADQZJUMhAkSYCBIEkqGQiSJMBAkCSVDARJEmAgSJJKBoIkCTAQJEklA0GSBBgIkqSSgSBJAgwESVLJQJAkAQaCJKlkIEiSAANBklQyECRJgIEgSSoZCJIkwECQJJUMBEkSYCBIkkoGgiQJMBAkSSUDQZIEGAiSpJKBIEkCDARJUslAkCQBBoIkqWQgSJIAA0GSVDIQJEmAgSBJKhkIkiTAQJAklQwESRJgIEiSSgaCJAkwEDTdli1w7rmwahUsWwYRxd9Vq4rxW7Y0XeFgsV/12K+BF5nZdA2VjY2N5aZNm5ouY+G54QZYvx6uuw527IBdu/aeZsmS4sN73HFw0UWwcmX/6xwU9qse+9W4iLgmM8fmm84thMXu4ovhmGPg6qth+/bZP6xQjN++vZjumGOK1y1G9qse+zVUDmi6ADXo4othwwa4++7qr5n64G7YUDxeu7YXlQ0m+1WP/Ro6fd1lFBHPAF4JPBZ4KHAL8K/A+Zn5zfle7y6jLrrhhuKb2PbtD3werRZcfz0cfnj36hpU9qse+zVQBnWX0UHANcBLgKcDrwaOBq6KiHafa1nc1q8v9ufujx074PnP7049g85+1WO/hlJfAyEzP5SZf5qZH83MKzLzIuC5wEOAU/tZy6K2eXNxgG+u/blV7doF11678M8OsV/12K+hNQgHlW8r/97TqwVMTk7S6XQYGRmh0+kwOTnZq0UNh3e9a/+/vU3ZubOY30Jmv+qxX13V1/VXZvZ9AJYAS4EjgY8B/wEcMt/rVq9enXVNTExkq9VK4L6h1WrlxMRE7XktGEcdlQkODsMxrFrV9CemMd1afwGbMudfNzfyO4SI2ASsLh9+D3h2Zn5rvtc9kIPKnU6Hbdu27TW+3W6zdevWWvNaMJYtK755ScNg6dLubXEMmW6tv6oeVG4qEFYBvwqsBF4BHAo8KTO3zjLtWcBZAKOjo6tna86+jIyMMNt7jAh2795du/YFIaLpCqTqImCRfla7tf4a1LOMAMjMb2Xm1Zn5IeBk4FeAV80x7cbMHMvMsRUrVtRe1ujoaK3xi8LSpU1XIFV34IFNV9CYfq+/Gj+onJl3Uuw2elQv5j8+Pk6r1dpjXKvVYnx8vBeLGw5eFkDD5Igjmq6gMf1efzUeCBFxKHAU8P1ezH/dunVs3LiRdrtNRNBut9m4cSPr1q3rxeKGw5o1xbVjumHJkuLCZM0feuzd8OIX268m+7VmTXfmNYT6vf7q9y+VPw58DdgM3AU8GvgT4OHAEzLzu/t6vb9U7pLNm+HEE/fvV6RTli+///ozC5X9qsd+DZxBPYZwFfAc4P3APwAvA64AjpsvDNRFj3tccVXJ/f0Wt2QJHH/8wv+w2q967NfQ8vLXi5XXmqnHftVjvwbKoG4haFCsXAkXXlhskj8Qy5cXr18sH1b7VY/9Gkpe/noxm7q08BlnzH3jkpmmbmRy4YWL79LE9qse+zV03EJY7NauLS4edsIJxbeyufb7LllSPH/CCcVm/GL9sNqveuzXUDEQVGzeX3llcTbH2WcX97hdurT4hejSpcXjs88unr/ySjfjZ/ZrNvbrfv77GhoeVJb212yXAhmiz5UWPg8qS5JqMRAkSYCBIEkqGQiSJMBAkCSVDARJEmAgSJJKBoIkCTAQJEklA0GSBBgIkqSSgSBJAgwESVKp8g1yIuJ04DRgFHjQjKczM4/oZmGSpP6qFAgR8VrgDcD1wHXAjl4WJUnqv6pbCGcCb8/MP+llMZKk5lQ9hvAw4NO9LESS1KyqgXAFcGwvC5EkNavqLqOXApdExG3APwK3z5wgM3d3szBJUn9VDYTvln8vnOP5rDEvSdIAqroSfyPFSl+StEBVCoTMPL/HdUiSGuYvlSVJgIEgSSoZCJIkwECQJJUMBEkSYCBIkkoGgiQJqBEIEbE7Iva67HVE3BMR93a3LElSv9W53MQHgNlW/B/ALQ1JGnqVAyEzN8wx/syuVSNJaozf7CVJQM0rlEbEGuBE4DDgB8CXM/PyXhQmSeqvqvdUPgj4CHASxVVP7wAeWjwVXwD+IDP3ukeCJGl4VN1l9P+A3wDWA8szcwWwHHgBMAa8vTflSZL6peouo98FXp2ZH5wakZn3AJPl1sMFvShOktQ/VbcQdgH/Nsdz3ymflyQNsaqB8EngD+d47o+AT3SnHElSU6ruMvo08LaI+AeKg8s/Bg4F1gJHA/8zIp46NXFmXtbtQiVJvVU1ED5a/n0k8DuzPP+x8m9QnIW0ZD/rkiT1WdVAWNPTKiRJjasUCJl5Ra8LkSQ1q+4vlQ8Gngg8DPh0Zt4eEQ8Cdmbm7l4UKEnqj0pnGUXhL4GbgU8B7wU65dOfBM7rSXWSpL6petrpq4GXAG8ETqA4eDzl08B/73JdkqQ+q7rL6IXAGzPzLyJi5hlE3wOO6G5ZkqR+q7qFcBhw1RzP7QQe3J1yJElNqRoIPwD+yxzPHQvc2J1yJElNqRoIHwFeFxG/NW1cRsSjgZcDH+56ZZKkvqoaCOcD3wa+yP0XufsIsKV8/JauVyZJ6quqP0y7OyJOAp4HPIPiQPJtwJuAycy8t2cVSpL6ovI9lTNzV2ZelJnPz8ynZ+Zpmfl+w2AB2bIFzj0XVq2CZcsgovi7alUxfsuWpiscLFP9mo392pv/vgZeZOb8E0XsAk7MzK/M8txq4CuZ2fML2o2NjeWmTZt6vZjF54YbYP16uO462LEDds1ye4slS4oP73HHwUUXwcqV/a9zUNiveuxX4yLimswcm2+6qlsIsY/nllBc4VTD6OKL4Zhj4OqrYfv22T+sUIzfvr2Y7phjitctRvarHvs1VPZ5DCEiRrg/DEbKx9Mtp7gc9q09qE29dvHFsGED3H139ddMfXA3bCger13bi8oGk/2qx34NnTl3GUXE64HXVZzPOzLzj+ddWMSpwGnAGHAIcBNwCfDmzPzZfK93l1EX3XBD8U1s+/YHPo9WC66/Hg4/vHt1DSr7VY/9GihVdxntawvhC1PzogiG91Bc3G66HcA3gb+vWNcrKELgz8t5HU9xSuuaiPhNr5jaR+vXF/tz98eOHfD858OVV3anpkFmv+qxX0NpzkAo74FwBUBEJPB3mfnD/Vze72bmLdMeXxERtwPvB04CvPVmP2zeXBzgm2t/blW7dsG11xZnhxxzTHdqG0T2qx77NbQqHVTOzDdMD4OI+E8RMRYRj6izsBlhMOWr5d/D6syrjsnJSTqdDiMjI3Q6HSYnJ3u1qOHwrnft/7e3KTt3FvNbyOxXPfarq/q6/srMWQeKH6C9ZZbx51HsKtpVDh8EDphrPvMNwDkUZymNzTft6tWrs66JiYlstVpZLiOBbLVaOTExUXteC8ZRR2WCg8NwDKtWNf2JaUy31l/Apsz518f7Oqj88SIv8rnTxj0N+AzFJSveDawCzgb+LDPfWjeMIuIw4Frg65n5tPmmfyAHlTudDtu2bdtrfLvdZuvWrbXmtWAsW1Z885KGwdKl3dviGDLdWn9146Dy8RSXppjuDOCXwDMy80flgqC4pEWtQIiIX6G429q95Xznmu4s4CyA0dHROosA4Kabbqo1flEwDDRM7rmn6Qoa0+/1176OIRwCfH/GuKcB/zIVBqV/AB5dZ6HlfZg/BaykCJeZZy/dJzM3ZuZYZo6tWLGizmKAuUPkgYTLgrF0adMVSNUdeGDTFTSm3+uvfQXCz5h245uIOBJ4GHvfKOcuil8rVxIRBwIfA54APCsze3oBk/HxcVqt1h7jWq0W4+PjvVzsYPOyABomRyzeGzL2e/21r0D4NnDKtMenUBzU+OyM6Q4HflxlYeUvnSeBk4FTMnOuu7B1zbp169i4cSPtdpuIoN1us3HjRtatW9frRQ+uNWuKa8d0w5IlxYXJmj/02LvhxS+2X032a82a7sxrCPV7/bWvg8rPofgV8SUUK/wNFPc+OD6nvag8+ByZ+Zx5FxbxToqzisbZ+8dsN+9r1xH4S+Wu2bwZTjxx/35FOmX58vuvP7NQ2a967NfA2e+L22XmJ4CXAr8BvIBiV9EfzAiDRwBrgH+sWNfvlH/PA748Y3hhxXlofz3uccVVJff3W9ySJXD88Qv/w2q/6rFfQ6vS5a8HhVsIXeS1ZuqxX/XYr4HS7ctfa6FZuRIuvLDYJH8gli8vXr9YPqz2qx77NZQq3UJTC9TUpYXPOGPuG5fMNHUjkwsvXHyXJrZf9divoeMWwmK3dm1x8bATTii+lc2133fJkuL5E04oNuMX64fVftVjv4aKxxB0vy1biguJXX45fP/7xS9EDzywOA98zRo45xwP8E1nv+qxX42pegzBQJCkBc6DypKkWioHQkQcHxGXRMStEXFvRDy+HP/miHhm70qUJPVDpUCIiCdR/HjsKIr7H0x/3W6KXx9LkoZY1S2Et1DcB+Fo4GUznvsa8PhuFiVJ6r+qv0N4PPDczMzy/srT3QrUvy61JGmgVN1C+CXQmuO5/wz8tDvlSJKaUjUQ/gV4aURM/1XJ1JbCmcBlXa1KktR3VXcZvRa4Evg68FGKMDg9Iv4KWE1xRVRJ0hCrtIWQmV8HnkxxX4TzgABeUj79lMz8Tm/KkyT1S+WL22Xm14CTy/shHwTcmZlduAOGJGkQ1L7aaWb+EvhhD2qRJDWociBExCrgVOCRwINmPJ2ZeXo3C5Mk9VelQIiIFwDvpTiY/BNg54xJhucKeZKkWdU5y+iTwJmZeWcP65EkNaRqIDwcOMcwkKSFq+oP064EVvWyEElSs6puIbwEuCQibgM+C9wxc4LM3N3NwiRJ/VU1EG4GrgUm5ng+a8xLkjSAqq7E/w74Q+ATwLfZ+ywjSdKQqxoIpwB/mplv72UxkqTmVD2o/Avgm70sRJLUrKqBcCHwvF4WIklqVtVdRtuA0yLin4FLmf0so/d2szBJUn9VDYR3ln/bwMmzPJ8Ul7aQJA2pqoFweE+rkCQ1rlIgZOa2XhciSWpW1YPKkqQFbs4thIi4Afi9zPx6RNzIvi9xnZl5RNerkyT1zb52GV0B3DXtv73ngSQtYHMGQmaeMe2/N/SlGklSY+Y8hhARN0TEsf0sRpLUnH0dVO4Ay/pUhySpYZ5lJEkC5g8EDyRL0iIx3w/T3hARt1aYT2bm6d0oSJLUjPkC4ThgR4X5uCUhSUNuvkB4TmZ+pS+VSJIa5UFlSRJgIEiSSgaCJAnY96UrDAtJWkRc6UuSAANBklQyECRJgIEgSSoZCJIkwECQJJUMBEkSYCBIkkoGgiQJMBAkSSUDQZIEGAiSpJKBIEkCDARJUslAkCQBDQRCRDwiIv46Ir4cEdsjIiOi0+86JEl7amIL4VHAWuAO4EsNLF+SNIsmAuGLmXloZj4L+EgDy5ckzaLvgZCZu/u9zMnJSTqdDiMjI3Q6HSYnJ/tdwlCxX/XYr3rsVz197VdmNjYALwQS6FSZfvXq1VnXxMREtlqtLJeTQLZarZyYmKg9r8XAftVjv+qxX/V0q1/Apqywjo1i2mZExAuBvwMOz8yt800/NjaWmzZtqrWMTqfDtm3b9hrfbrfZunXeRS469qse+1WP/aqnW/2KiGsyc2ze6QY9ECLiLOAsgNHR0dWzNWdfRkZGmO09RgS7d/d979XAs1/12K967Fc93epX1UAY+N8hZObGzBzLzLEVK1bUfv3o6Git8Yud/arHftVjv+rpd78GPhD21/j4OK1Wa49xrVaL8fHxhioabParHvtVj/2qp+/9qnKgoVcDfTionFkcmGm32xkR2W63PYA1D/tVj/2qx37V041+McgHlSPi1PI/TwbOAc4FbgFuycwr5nrdAzmoLEmLXdVjCAf0o5hZzPxB2jvKv1cAJ/W3FEkSNBQImRlNLFeSNLcFf1BZklSNgSBJAgwESVLJQJAkAQaCJKlkIEiSAANBklQyECRJgIEgSSoZCJIkwECQJJUMBEkSYCBIkkoGgiQJMBAkSSUDQZIEGAiSpJKBIEkCDARJUslAkCQBBoIkqWQgSJIAA0GSVDIQJEmAgSBJKhkIkiTAQJAklQwESRJgIEiSSgaCJAkwECRJJQNBkgQYCJKkkoEgSQIMBElSyUCQJAEGgiSpZCBIkgADQZJUMhAkSYCBIEkqGQiSJMBAkCSVDARJEmAgSJJKBoIkCTAQJEklA0GSBBgIkqSSgSBJAgwESVLJQJAkAQaCJKlkIEiSAANBklQyECRJgIEgSSoZCJIkwECQJJUMBEkSYCBIkkp9D4SIeGREfDQifhoRd0XEJREx2u86JEl76msgREQLuAw4CjgdWA8cCVweEQ/uZy2SpD0d0OflvQhYCTwmM78HEBGbgX8Dzgb+qs/1SJJK/d5l9GzgqqkwAMjMG4ErgVN6tdDJyUk6nQ4jIyN0Oh0mJyd7tagFwX7VY7/qsV/19LVfmdm3AfgR8LezjH8HcMt8r1+9enXWNTExka1WK4H7hlarlRMTE7XntRjYr3rsVz32q55u9QvYlFXW0VUm6tYA7ATeMsv4C4B753v9AwmEdru9RzOnhna7XXtei4H9qsd+1WO/6ulWv6oGQhTT9kdE7ATempmvnjF+HHhlZu51TCMizgLOAhgdHV29bdu2WsscGRlhtvcYEezevbvWvBYD+1WP/arHftXTrX5FxDWZOTbv8uqVt9/uAA6aZfxDy+f2kpkbM3MsM8dWrFhRe4Gjo7Of0TrX+MXOftVjv+qxX/X0u1/9DoRvAEfPMv6xwDd7scDx8XFardYe41qtFuPj471Y3NCzX/XYr3rsVz1971eV/UrdGoCXAvcCK6eN6wD3AC+f7/UP5BhCZnFgpt1uZ0Rku932ANY87Fc99qse+1VPN/rFgB5DeDDwdeBu4DUUB0jeBDwEeFxm/nxfrx8bG8tNmzb1vE5JWkgG8hhCZv4CeCrwXeAiYBK4EXjqfGEgSeqtfv9Smcy8Cfj9fi9XkrRvXu1UkgQYCJKkkoEgSQIMBElSyUCQJAEGgiSpZCBIkgADQZJU6uulK/ZXRNwC1Lv+9Z4OBm7tUjmLgf2qx37VY7/q2Z9+tTNz3stFD1Ug7K+I2FTleh4q2K967Fc99quefvTLXUaSJMBAkCSVFlsgbGy6gCFjv+qxX/XYr3p63q9FdQxBkjS3xbaFIEmaw4IOhIh4ZER8NCJ+GhF3RcQlEeHdvOcQEY+IiL+OiC9HxPaIyIjoNF3XoIqIUyPiYxGxLSLujojvRMRfRMRDmq5tEEXEMyLisoj4UUTsiIibI+LiiHhs07UNi4i4tPxcXtCL+S/YQIiIFnAZcBRwOrAeOBK4vLyVp/b2KGAtcAfwpYZrGQavAHYBfw48E3gn8GLgnyNiwX629sNBwDXAS4CnA68Gjgauioh2k4UNg4g4DTi2l8vo+x3T+uhFwErgMZn5PYCI2Az8G3A28FcN1jaovpiZhwJExAspPrSa2+9m5i3THl8REbcD7wdOovhColJmfgj40PRxEfEV4NvAqcBbm6hrGETErwFvA/4E+GCvlrOQv8U8G7hqKgwAMvNG4ErglMaqGmCZubvpGobJjDCY8tXy72H9rGWI3Vb+vafRKgbf/wa+UYZqzyzkQDgauH6W8d8A3GepXnlK+fdbjVYxwCJiSUQsjYgjgb8FfgR8uOGyBlZEPAl4AXBur5e1kHcZHUSxL3ym24GH9rkWLQIRcRjwRuBzmbmp6XoG2NXA6vK/vwc8NTN/0mA9AysiDqQIzf+Tmd/p9fIW8hYCwGw/soi+V6EFLyJ+BfgkcC9wRsPlDLr1wBOB5wF3URyE7zRZ0AB7JbAcGO/HwhbyFsIdFFsJMz2U2bccpAckIh4EfIriJIanZObNDZc00DJzanfa1RHxT8BW4FXAOY0VNYDKU+TPA14ILIuIZdOeXlYeaP5ZZu7q1jIX8hbCNyiOI8z0WOCbfa5FC1S5Sf8x4AnAszJzS8MlDZXMvJNit9Gjmq5lAK0EHgRMUHyJnRqgOOX5DuCYbi5wIQfCp4AnRsTKqRHlZulvlc9J+6X8rcEkcDJwSmZe1XBJQyciDqX4rdD3m65lAF0HrJllgCIk1lCEadcs2GsZlT8++zpwN/AaiuMJbwIeAjwuM3/eYHkDKyJOLf/zZIpN+HOBW4BbMvOKxgobQBHxTooejQN/P+Ppm911tKeI+DjwNWAzxbGDR1OcV/9w4AmZ+d0GyxsaEZHAeGa+puvzXqiBAPftg3sb8DSKg8mfB16amVubrGuQlf/YZnNFZp7Uz1oGXURsBeb6he0bMvP8/lUz+CLilRS/hD8CWAr8O/AF4C/8TFZnIEiSem4hH0OQJNVgIEiSAANBklQyECRJgIEgSSoZCJIkwEBQwyLiORHxxYj4SXkbym0R8YmIeGbTte2PKu8rIk6KiPMH6e5qZT2ei75IDcw/RC0+EfE/gI9T3MXuTOC/AVP3in1qU3Xtrxrv6yTg9QzW5/DdwIlNF6Fm+MM0NSYibgKuyczfm+W5kWG9g1vV9xUR51MEwoGZeW9/q5T2NkjfTLT4HERxt6y9TA+DiNgQERkRTy53u/w8Im6LiL+JiOXTXxcRb4iIr0XETyPi1oi4LCKeOHP+EbEiIt4REf8eETvKvxdNv8RwRBwbEZ+KiDvK3T5XRsR/7cb7mhYGAPeU7+++b2cR0YqI/xURN0bEzvLvedN3L5W7nDIifj8i3lfWeVdETEbEwyrUuRd3GS1uC/l+CBp8XwFOj4gbgE9WuLjZBHAx8A6Ky02/DngwsGHaNIdRXL/q5vK55wNfjIixzNwMEBEPBf6VYsV9AcXF1g6huNf2UmBHRDwe+BJwLfAiYDvFhew+FxG/mZnX7Of7ejfwCIpdSk8C7rumfUQcAHyG4lLtbwK2UNxQ5rVlzS+fMa//C3wOOA04Engz8Ovcf2VMIuJ9wOmZ6Q2iNLfMdHBoZKC42uVmiivRJnAr8CHg6TOm21A+/64Z48+jWJE+eo75L6H40vMd4O3Txr+xfN3x+6jt8xT3RV46Y37fAj7Rpfd1fvn8ATPGry/HP3mW97sTOKR8fFI53aUzpltXjj952rj3APdW+H9yfrFaaP7fh0P/B3cZqTFZfHM+nuLG9OMU13//PeAzETHblRwvnvH4wxS7PZ8wNSIifjsiLo+I2yhuZ3kPxQr6MdNe93Tgq5l57Wx1lbuhngJ8BNgdEQeU39qD4pv4k7v8vmZ6JrAN+NepZZfL/yxwIMXWwnQz+/IRYDfTDg5n5pmZ6R4B7ZP/QNSoLG7/98VyICJ+HbgUeH1E/E1mTr/d6Y9nvHzq8WHlax8P/CPF7pYzgf+g2BJ4N8Wdp6Y8jOJeGXM5iGJr4LXlsJf5DnrXfF8zHUJxWe175nh+5vGBPfqSmTsj4g7KvkhVGQgaKJn5w4h4N/B2iv3hX5n29KEUt0ad/hjgB+Xf36fYKnhuZt63Mi2PGdw57XW3su+V5Z0U37D/BvjAHHXWOgNqnvc1023AjRT3DpjN1hmPD53+ICKWUtw7/AdINRgIakxEPDIz/32Wp44q/848U2ctcNm0x39EseKeWrm2KLYIpp+t81RglGIFO+WzwGsi4tjM3GtLITN/ERFfAo4FvlZ35V/jfe0o/y4HfjZtukspwu3nmfntCotcC7x32uM/oNiV9uXKRUsYCGrW9RFxOcWPuG4EfhV4FsXZPBdn5k0zpn9WRPwlxQr9CRSnbX4g7z+L51LgpcD7IuJCimMHr2Xvb8pvA55HccbQBRRn8RxMcZbROZn5M+BlFLt7PhMR76HY/XQw8HhgSWa+qgvv65vl35dHxD8BuzJzE8V9ms8APh8Rb6XYvbWU4k5jzwaek5nbpy3v6PL9frh8z+MUd7j7/NQE5Xs43eMI2qemj2o7LN6BYgX5KYoxEhcMAAAA2UlEQVQDqL8EfkFxmuefsefZPRsoz7oBPgn8HLidYpfO8hnz/GOKlfDdwFeB36a4TeMXZkx3CLCRYkW/k+J2ju8Hlk2bZhXFSvYnFN/mby7rfVaX3teS8j38hGJLJ6c99yCKM36+XS779vL9nE95VhL3n2X0XOB9FLu6fgZ8EDh4Rk3vo8LZQ3iW0aIe/KWyBl5EbAAuBI7MzO81XM7AiIiTgMuBp2Xm5xouRwuAp51KkgADQZJUcpeRJAlwC0GSVDIQJEmAgSBJKhkIkiTAQJAklQwESRIA/x8+HVii8L3F+AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(6,6))\n",
    "nx=ny=5\n",
    "x = np.arange(nx)\n",
    "y=np.zeros(ny)\n",
    "\n",
    "for iy in range(nx):\n",
    "    plt.plot(x,y,'ko')\n",
    "    y+=1\n",
    "plt.plot(1,2,'ro',ms=20)\n",
    "plt.plot(2,2,'ro',ms=20)\n",
    "plt.plot(3,2,'ro',ms=20)\n",
    "plt.plot(1,3,'ro',ms=20)\n",
    "plt.plot(2,3,'ro',ms=20)\n",
    "plt.plot(3,3,'ro',ms=20)\n",
    "plt.plot([1,3],[2,2],'r-',linewidth=5)\n",
    "plt.plot([1,3],[3,3],'r-',linewidth=5)\n",
    "plt.plot([2,2],[2,3],'r-',linewidth=5)\n",
    "\n",
    "plt.ylabel('Time Step: n',fontsize=16)\n",
    "plt.yticks(range(5),fontsize=16)\n",
    "plt.xlabel('Space Step: i',fontsize=16)\n",
    "plt.xticks(range(5),fontsize=16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 5. Illustration of the Crank-Nicholson stencil where known values at time step $n=2$ are used to calculate solutions at time step $n=3$.**\n",
    "\n",
    "As we see in the above figure, we now have a system where we have three unknown at time step $n+1$ for each location of the stencil.  However, if you shift the stencil one step to the left and right, you will also have constraint equations that can be used to help solve for the unknown forward value at $[i,n]=[2,3]$.\n",
    "\n",
    "### CN Numerical Scheme\n",
    "\n",
    "The CN numerical scheme is based on taking an equally weighted average of the points at time step $n$ and $n+1$ for approximating the second spatial derivative. We can write the CN numerical scheme in the following way:\n",
    "\n",
    "$$\\frac{\\phi^{n+1}_i -\\phi^{n}_i}{\\Delta t} = \n",
    "\\frac{\\kappa}{2}\n",
    "\\left( \n",
    "\\frac{\n",
    "\\phi^{n+1}_{i+1}-2\\phi^{n+1}_i+\\phi^{n+1}_{i-1} + \n",
    "\\phi^{n}_{i+1}-2\\phi^{n}_i+\\phi^{n}_{i-1}\n",
    "}{\\Delta x^2}\n",
    "\\right)\n",
    ".\\tag{23}$$\n",
    "\n",
    "Again, by defining $\\alpha=\\frac{\\kappa \\Delta t}{\\Delta x^2}$, we can rearrange equation 23 as\n",
    "\n",
    "$$2\\left(\\phi^{n+1}_i -\\phi^{n}_i\\right) = \n",
    "\\alpha\n",
    "\\left( \n",
    "\\phi^{n+1}_{i+1}-2\\phi^{n+1}_i+\\phi^{n+1}_{i-1} + \n",
    "\\phi^{n}_{i+1}-2\\phi^{n}_i+\\phi^{n}_{i-1}\n",
    "\\right),\n",
    "\\tag{24}$$\n",
    "\n",
    "and then move the $n+1$ and $n$ terms to the left- and right-hand sides, respectively:\n",
    "\n",
    "$$-\\alpha\\phi^{n+1}_{i+1}+2(1+\\alpha)\\phi^{n+1}_{i}-\\alpha\\phi^{n+1}_{i-1}=\n",
    "\\alpha\\phi^{n}_{i+1}+2(1-\\alpha)\\phi^{n}_{i}+\\alpha\\phi^{n}_{i-1}. \\tag{25}$$\n",
    "\n",
    "While this seems quite a bit more complicated, we have actually a similar system in Module 8!  To illustrate this, let's write this in a linear algebra framework in the standard form of ${\\bf Ax}={\\bf b}$:\n",
    "\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{ccccccc}\n",
    "2(1+\\alpha) & -\\alpha & 0 & 0 & ... & 0&0 &0 &0 \\\\\n",
    "-\\alpha & 2(1+\\alpha) & -\\alpha & 0 &... & 0&0  &0 &0 \\\\\n",
    "0 & -\\alpha & 2(1+\\alpha) & -\\alpha & ... & 0&0 &0 &0 \\\\\n",
    "\\vdots & \\vdots & \\vdots& \\vdots & \\ddots & \\vdots &\\vdots &\\vdots & \\vdots\\\\\n",
    "0 &0&0&0&...&-\\alpha&2(1+\\alpha)&-\\alpha&0\\\\\n",
    "0 &0&0&0&...&0&-\\alpha&2(1+\\alpha)&-\\alpha\\\\\n",
    "0 &0&0&0&...&0&0&-\\alpha&2(1+\\alpha)\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "u^{n+1}_{1} \\\\\n",
    "u^{n+1}_{2} \\\\\n",
    "u^{n+1}_{3} \\\\\n",
    "\\vdots\\\\\n",
    "u^{n+1}_{I-3} \\\\\n",
    "u^{n+1}_{I-2} \\\\\n",
    "u^{n+1}_{I-1} \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "\\alpha\\phi^{n}_{2}+2(1-\\alpha)\\phi^{n}_{1}+\\alpha(\\phi^{n}_{0}+\\phi^{n+1}_{0}) \\\\\n",
    "\\alpha\\phi^{n}_{3}+2(1-\\alpha)\\phi^{n}_{2}+\\alpha\\phi^{n}_{1} \\\\\n",
    "\\alpha\\phi^{n}_{4}+2(1-\\alpha)\\phi^{n}_{3}+\\alpha\\phi^{n}_{2} \\\\\n",
    "\\vdots\\\\\n",
    "\\alpha\\phi^{n}_{I-2}+2(1-\\alpha)\\phi^{n}_{I-3}+\\alpha\\phi^{n}_{I-4} \\\\\n",
    "\\alpha\\phi^{n}_{I-1}+2(1-\\alpha)\\phi^{n}_{I-2}+\\alpha\\phi^{n}_{I-3} \\\\\n",
    "\\alpha(\\phi^{n}_{I  }+\\phi^{n+1}_{I})+2(1-\\alpha)\\phi^{n}_{I-1}+\\alpha\\phi^{n}_{I-2} \\\\\n",
    "\\end{array}\n",
    "\\right]. \\tag{26}\n",
    "$$\n",
    "\n",
    "Here we see that modeling matrix ${\\bf A}$ is a **banded matrix** on the main, super- and subdiagonal. Thus, this system of equations again represents a **tridiagonal system** that can be solved using standard methods of linear algebra!  Also, note the first and last entries in known solution vector ${\\bf b}$ contain the terms $\\phi^n_0$, and $\\phi^{n+1}_0$ (first row) and $\\phi^n_{I}$ and $\\phi^{n+1}_{I}$ (last row). This is because these values are known and are used to represent the **boundary conditions**.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solution by Thomas Algorithm\n",
    "\n",
    "Previously, we have solved systems of equations like those in equation 26 using a \"brute force\" method of *np.linalg.solve(Ax,bx)*, which computes the “exact” solution, ${\\bf x}$, of the well-determined, i.e., full rank, linear matrix equation ${\\bf Ax} = {\\bf b}$. However, it should be pretty clear that the FD matrix ${\\bf A}$ is actually very **sparse** and has **banded structure**.  Thus, it is judicious to ask whether there are any more efficient numerical methods that can take advantage of this matrix structure.  \n",
    "\n",
    "Since you have seen the title of this subsection, you have probably surmised that the answer is yes.  Let's investigate the [Thomas algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm), which can be used to solve certain tridiagonal systems (i.e., a [diagonally dominant matrix](https://en.wikipedia.org/wiki/Diagonally_dominant_matrix)) defined by:\n",
    "\n",
    "$$a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i}, \\tag{27a}$$\n",
    "\n",
    "where ${a_{1}=0}$ and $c_{n}=0$ and written in matrix form as\n",
    "\n",
    "$${\\begin{bmatrix}{b_{1}}&{c_{1}}&{}&{}&{0}\\\\{a_{2}}&{b_{2}}&{c_{2}}&{}&{}\\\\{}&{a_{3}}&{b_{3}}&\\ddots &{}\\\\{}&{}&\\ddots &\\ddots &{c_{n-1}}\\\\{0}&{}&{}&{a_{n}}&{b_{n}}\\\\\\end{bmatrix}}{\\begin{bmatrix}{x_{1}}\\\\{x_{2}}\\\\{x_{3}}\\\\\\vdots \\\\{x_{n}}\\\\\\end{bmatrix}}={\\begin{bmatrix}{d_{1}}\\\\{d_{2}}\\\\{d_{3}}\\\\\\vdots \\\\{d_{n}}\\\\\\end{bmatrix}}.\\tag{28b}$$\n",
    "\n",
    "This algorithm uses a two-step approach involving a forward and then a backward \"sweep\" calculation. The forward sweep consists of modifying the coefficients as follows, denoting the new coefficients with primes:\n",
    "\n",
    "$${\\displaystyle c'_{i}={\\begin{cases}{\\begin{array}{lcl}{\\cfrac {c_{i}}{b_{i}}}&;&i=1\\\\{\\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}}&;&i=2,3,\\dots ,n-1\\\\\\end{array}}\\end{cases}}\\,} \\tag{29a}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$d'_{i}=\n",
    "{\\begin{cases}\n",
    "{\\begin{array}{lcl}{\\cfrac {d_{i}}{b_{i}}}&;&i=1\\\\{\\cfrac {d_{i}-a_{i}d'_{i-1}}{b_{i}-a_{i}c'_{i-1}}}&;&i=2,3,\\dots ,n.\\\\\n",
    "\\end{array}}\n",
    "\\end{cases}} \\tag{29b}\n",
    "$$\n",
    "The solution is then obtained by back substitution:\n",
    "\n",
    "$$ x_{n}=d'_{n} \\tag{30a}$$\n",
    "\n",
    "and\n",
    "\n",
    "$$ x_{i}=d'_{i}-c'_{i}x_{i+1}\\quad \\ i=n-1,n-2,\\ldots ,1. \\tag{30b}$$\n",
    "\n",
    "Let's now implement the Thomas algorithm using two different functions.  The first function *TDMAsolver(a,b,c,d)* implements the Thomas Algorithm solution, while the second function *Solve_Tridiagonal_homogeneous_Thomas(a,b,c,d)* sets up the coefficients to be passed to the *TDMAsolver*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def TDMAsolver(a, b, c, d):\n",
    "    '''\n",
    "    TDMA solver, a b c d can be NumPy array type or Python list type.\n",
    "    refer to http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm\n",
    "    '''\n",
    "    nf = len(d) # number of equations\n",
    "    ac, bc, cc, dc = map(np.array, (a, b, c, d)) # copy arrays\n",
    "\n",
    "    ## . . Forward sweep\n",
    "    for it in range(1, nf):\n",
    "        mc = ac[it-1]/bc[it-1]\n",
    "        bc[it] = bc[it] - mc*cc[it-1] \n",
    "        dc[it] = dc[it] - mc*dc[it-1]\n",
    "    \n",
    "    ## . . Backsubstitution sweep\n",
    "    xc = bc\n",
    "    xc[-1] = dc[-1]/bc[-1]\n",
    "\n",
    "    for il in range(nf-2, -1, -1):\n",
    "        xc[il] = (dc[il]-cc[il]*xc[il+1])/bc[il]\n",
    "\n",
    "    return xc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Solve_Diffusion_homogeneous_Thomas(U,K,dx,dt,LB,RB):\n",
    "    '''Set up tridiagonal matrix and call Thomas Algorithm \n",
    "    usage: x = Solve_Diffusion_homogeneous_Thomas(U,K,dx,dt,LB,RB):\n",
    "    input: \n",
    "        U: heat solution at time step n   (nx)\n",
    "        K : thermal conductivity\n",
    "        dx: spatial sampling\n",
    "        dt: temporal sampling\n",
    "        LB: Left  boundary condition (Dirchelet)\n",
    "        RB: Right boundary condition (Dirchelet)\n",
    "    output:\n",
    "        U: heat solution at time step n+1 (nx)\n",
    "    depends on:\n",
    "        TDMAsolver\n",
    "    \n",
    "    written by Jeff Shragge, jshragge@mines.edu, 10/2019\n",
    "    '''\n",
    "    n = len(U)\n",
    "    \n",
    "    alpha = K*dt/(dx*dx)\n",
    "    \n",
    "    ## . . Set up coefficients\n",
    "    aa =     -alpha *np.ones(n)\n",
    "    bb = 2*(1+alpha)*np.ones(n)\n",
    "    cc =     -alpha *np.ones(n)\n",
    "    \n",
    "    ## . . Treat ends\n",
    "    aa[0  ]=0\n",
    "    cc[n-1]=0\n",
    "    \n",
    "    ## . . Set up known right-hand side solution matrix\n",
    "    dd = np.zeros((nx))\n",
    "    \n",
    "    dd[1:nx-2] =     alpha *U[2:nx-1]+\\\n",
    "                2*(1-alpha)*U[1:nx-2]+\\\n",
    "                     alpha *U[0:nx-3]\n",
    "    dd[1   ]+=alpha*U[0   ]\n",
    "    dd[nx-2]+=alpha*U[nx-1]\n",
    "    \n",
    "    ## . . Call Solver to update solution\n",
    "    U = TDMAsolver(aa, bb, cc, dd)\n",
    "    \n",
    "    return U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now redo the scenario from above that proved to be unstable with $dt=2$ and $\\alpha=0.55$m$^2$/s, but go back to simple Dirchelet boundary conditions where $\\phi(x=0,t)=\\phi(x=1,t)=0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha (Copper) = 0.55\n",
      "alpha (Iron  ) = 0.115\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_5.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=400,2\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 51\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/2,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kc       = 1.1*10**(-4) ## . . thermal diffusivity\n",
    "Ki       = 2.3*10**(-5) ## . . thermal diffusivity\n",
    "LB,RB = 0,0 ## . . Derivative Boundary conditions on left and right\n",
    "\n",
    "print('alpha (Copper) =',Kc*dt/(dx*dx))\n",
    "print('alpha (Iron  ) =',Ki*dt/(dx*dx))\n",
    "## . . Iterate over solution\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "for i in range(nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = Solve_Diffusion_homogeneous_Thomas(Uc,Kc,dx,dt,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = Solve_Diffusion_homogeneous_Thomas(Ui,Ki,dx,dt,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "\n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig5,ax1 = plt.subplots(1)\n",
    "fig5.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig5.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim5 = animation.FuncAnimation(fig5,diffusion_1D_animate,frames=int((nt-kskip)/kskip),interval=100)\n",
    "anim5.save('./movies/Ex7_5.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_5.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 6. Movie showing the same example as in Figure 4, but using a Crank-Nicholson solver with dt=2.0s.  We see that the CN approach does not suffer from any (apparent) numerical instabilities.**\n",
    "\n",
    "Let's now rerun the solution above, but increase the time step to, say, $\\Delta t=20$s:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha (Copper) = 5.5\n",
      "alpha (Iron  ) = 1.15\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_6.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=400,20\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 51\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/2,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kc       = 1.1*10**(-4) ## . . thermal diffusivity\n",
    "Ki       = 2.3*10**(-5) ## . . thermal diffusivity\n",
    "LB,RB = 0,0 ## . . Derivative Boundary conditions on left and right\n",
    "\n",
    "print('alpha (Copper) =',Kc*dt/(dx*dx))\n",
    "print('alpha (Iron  ) =',Ki*dt/(dx*dx))\n",
    "\n",
    "## . . Iterate over solution\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "for i in range(nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = Solve_Diffusion_homogeneous_Thomas(Uc,Kc,dx,dt,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = Solve_Diffusion_homogeneous_Thomas(Ui,Ki,dx,dt,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "    \n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig6,ax1 = plt.subplots(1)\n",
    "fig6.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig6.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim6 = animation.FuncAnimation(fig6,diffusion_1D_animate,frames=int((nt-kskip)/kskip),interval=100)\n",
    "anim6.save('./movies/Ex7_6.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_6.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 7. Movie showing the same example as in Figure 4, but using a Crank-Nicholson solver with ∆t=20.0s.  We see that the CN approach does not suffer from any (apparent) numerical instabilities.**\n",
    "\n",
    "### Importance of CN method\n",
    "\n",
    "There are a number of different factors that have come to light in this series of investigations:\n",
    "\n",
    "   * **Unconditional Stability** - The CN method allows you to run simulations that are guaranteed to be stable, which can be shown to be the following through a von Neumann staility analysis for all values of $\\alpha$:\n",
    "\n",
    "$$g^2 = \\left| \\frac{1-2\\alpha \\sin^2(\\gamma/2)}{1+2\\alpha \\sin^2(\\gamma/2)} \\right|^2 \\le 1\\tag{31}$$\n",
    "\n",
    "   * **Large time-stepping** - We have seen that we can use large $\\Delta t$ time steps, which means that it is much quicker to forward simulate the numerical solution to the steady-state solution (if there is one!).\n",
    "   * **Reduced computational complexity** - The computational complexity of the Thomas Algorithm is $\\mathcal{O}(nx)$ compared to $O(nx^3)$ for Gaussian elimination, which means that this approach scales much better as $nx$ gets very large."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the 1D Convection-Diffusion equation\n",
    "\n",
    "The 1D convection-diffusion equation for a convective velocity $v_x$ and diffusivity $\\kappa$ is given by:\n",
    "\n",
    "$$\\frac{\\partial \\phi}{\\partial t} + v_x \\frac{\\partial \\phi}{\\partial x} = \\kappa \\frac{\\partial^2 \\phi}{\\partial x^2} \\tag{32}$$\n",
    "\n",
    "which is subject to some initial condition\n",
    "\n",
    "$$\\phi(x,t=0)= \\phi_0(x). \\tag{33}$$\n",
    "\n",
    "Here, both the convection velocity $v_x$ and diffusivity $\\kappa$ are assumed to be constant.\n",
    "\n",
    "## Crank-Nicholson Numerical Approach\n",
    "\n",
    "Let's take the Crank-Nicholson approach that we learned in the section above and apply it to the the 1D convection-diffusion equation.  This leads to the following discretized equation\n",
    "\n",
    "$$\\frac{\\phi^{n+1}_i-\\phi^n_i}{\\Delta t} + \\frac{v_x}{2} \\left(\\frac{\\phi^{n+1}_{i+1}-\\phi^{n+1}_{i-1}}{2\\Delta x}+\\frac{\\phi^{n}_{i+1}-\\phi^{n}_{i-1}}{2\\Delta x} \\right) \n",
    "=\n",
    "\\frac{\\kappa}{2}\n",
    "\\left(\n",
    " \\frac{\n",
    "\\phi^{n+1}_{i+1}-2\\phi^{n+1}_i+\\phi^{n+1}_{i-1} + \n",
    "\\phi^{n}_{i+1}-2\\phi^{n}_i+\\phi^{n}_{i-1}\n",
    "}{\\Delta x^2}\n",
    "\\right)\n",
    "\\tag{34}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now multiply through by $\\Delta t$ \n",
    "\n",
    "$$\\phi^{n+1}_i - \\phi^n_i + \\frac{v_x \\Delta t}{4\\Delta x} \\left(\\phi^{n+1}_{i+1}-\\phi^{n+1}_{i-1}+\\phi^{n}_{i+1}-\\phi^{n}_{i-1} \\right) \n",
    "=\n",
    "\\frac{\\kappa \\Delta t}{2\\Delta x^2}\n",
    "\\left(\n",
    "\\phi^{n+1}_{i+1}-2\\phi^{n+1}_i+\\phi^{n+1}_{i-1} + \n",
    "\\phi^{n}_{i+1}-2\\phi^{n}_i+\\phi^{n}_{i-1}\n",
    "\\right),\n",
    "\\tag{35a}$$\n",
    "\n",
    "and use the following subsitutions of the Courant number $C = \\frac{v_x \\Delta t}{\\Delta x}$ and the diffusion number $\\alpha = \\frac{\\kappa \\Delta t}{\\Delta x^2}$ and multiply through by 4 such that \n",
    "\n",
    "$$4\\phi^{n+1}_i - 4\\phi^n_i + C \\left(\\phi^{n+1}_{i+1}-\\phi^{n+1}_{i-1}+\\phi^{n}_{i+1}-\\phi^{n}_{i-1} \\right) \n",
    "=\n",
    "2\\alpha\n",
    "\\left(\n",
    "\\phi^{n+1}_{i+1}-2\\phi^{n+1}_i+\\phi^{n+1}_{i-1} + \n",
    "\\phi^{n}_{i+1}-2\\phi^{n}_i+\\phi^{n}_{i-1}\n",
    "\\right).\n",
    "\\tag{35b}$$\n",
    "\n",
    "Let's now bring all of the $n+1$ terms to the left and all of the $n$ terms to the right:\n",
    "\n",
    "$$4\\phi^{n+1}_i + C \\phi^{n+1}_{i+1}-C\\phi^{n+1}_{i-1} \n",
    "-2\\alpha\\phi^{n+1}_{i+1}+4\\alpha\\phi^{n+1}_i-2\\alpha\\phi^{n+1}_{i-1} \n",
    "=\n",
    "4\\phi^n_i\n",
    "-C\\phi^{n}_{i+1}\n",
    "+C\\phi^{n}_{i-1}\n",
    "+2\\alpha\\phi^{n}_{i+1}-4\\alpha\\phi^{n}_i+2\\alpha\\phi^{n}_{i-1}.\n",
    "\\tag{35c}$$\n",
    "\n",
    "Now let's collect like terms to yield\n",
    "\n",
    "$$(C-2\\alpha) \\, \\phi^{n+1}_{i+1}\n",
    "+4(1+\\alpha) \\, \\phi^{n+1}_i \n",
    "-(C+2\\alpha)\\,\\phi^{n+1}_{i-1} \n",
    "=\n",
    "- (C-2\\alpha)\\,\\phi^{n}_{i+1}\n",
    "+4(1-\\alpha)\\,\\phi^n_i\n",
    "+(C+2\\alpha)\\,\\phi^{n}_{i-1}.\n",
    "\\tag{35d}$$\n",
    "\n",
    "Again, we see that even though the resulting discretized system is more complex, it is still represented by a tridiagonal system!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's look at the resulting system in a linear algebra framework\n",
    "\n",
    "$$\n",
    "\\left[\n",
    "\\begin{array}{ccccccc}\n",
    "4(1+\\alpha) & (C-2\\alpha) & 0 & 0 & ... & 0&0 &0 &0 \\\\\n",
    "-(C+2\\alpha) & 4(1+\\alpha) & (C-2\\alpha) & 0 &... & 0&0  &0 &0 \\\\\n",
    "0 &-(C+2\\alpha) & 4(1+\\alpha) & (C-2\\alpha)& ... & 0&0 &0 &0 \\\\\n",
    "\\vdots & \\vdots & \\vdots& \\vdots & \\ddots & \\vdots &\\vdots &\\vdots & \\vdots\\\\\n",
    "0 &0&0&0&...&-(C+2\\alpha)&4(1+\\alpha)&(C-2\\alpha)&0\\\\\n",
    "0 &0&0&0&...&0&-(C+2\\alpha)&4(1+\\alpha)&(C-2\\alpha)\\\\\n",
    "0 &0&0&0&...&0&0&-(C+2\\alpha)&4(1+\\alpha)\\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "u^{n+1}_{1} \\\\\n",
    "u^{n+1}_{2} \\\\\n",
    "u^{n+1}_{3} \\\\\n",
    "\\vdots\\\\\n",
    "u^{n+1}_{I-3} \\\\\n",
    "u^{n+1}_{I-2} \\\\\n",
    "u^{n+1}_{I-1} \\\\\n",
    "\\end{array}\n",
    "\\right]\n",
    "=\n",
    "\\left[\n",
    "\\begin{array}{c}\n",
    "- (C-2\\alpha)\\,\\phi^{n}_{2}+4(1-\\alpha)\\,\\phi^n_1+(C+2\\alpha)\\,(\\phi^{n}_{0}+ \\phi^{n+1}_{0})\\\\\n",
    "- (C-2\\alpha)\\,\\phi^{n}_{3}+4(1-\\alpha)\\,\\phi^n_2+(C+2\\alpha)\\,\\phi^{n}_{1} \\\\\n",
    "- (C-2\\alpha)\\,\\phi^{n}_{4}+4(1-\\alpha)\\,\\phi^n_3+(C+2\\alpha)\\,\\phi^{n}_{2} \\\\\n",
    "\\vdots\\\\\n",
    "- (C-2\\alpha)\\,\\phi^{n}_{I-2}+4(1-\\alpha)\\,\\phi^n_{I-3}+(C+2\\alpha)\\,\\phi^{n}_{I-4} \\\\\n",
    "- (C-2\\alpha)\\,\\phi^{n}_{I-1}+4(1-\\alpha)\\,\\phi^n_{I-2}+(C+2\\alpha)\\,\\phi^{n}_{I-3} \\\\\n",
    "- (C-2\\alpha)\\,(\\phi^{n}_{I  }+\\phi^{n+1}_{I})+4(1-\\alpha)\\,\\phi^n_{I-1}+(C+2\\alpha)\\,\\phi^{n}_{I-2} \\\\\n",
    "\\end{array}\n",
    "\\right]. \\tag{36}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Solve_Convection_Diffusion_homogeneous_Thomas(U,K,vx,dx,dt,LB,RB):\n",
    "    '''Set up tridiagonal matrix and call Thomas Algorithm \n",
    "    usage: x = Solve_Tridiagonal(a,b,c,u):\n",
    "    input: \n",
    "        u: heat solution at time step n   (nx)\n",
    "        K : thermal diffusivity\n",
    "        vx: convection velocity\n",
    "        dx: spatial sampling\n",
    "        dt: temporal sampling\n",
    "        LB: Left  boundary condition (Dirchelet)\n",
    "        RB: Right boundary condition (Dirchelet)\n",
    "    output:\n",
    "        u: heat solution at time step n+1 (nx)\n",
    "    depends on:\n",
    "        TDMAsolver\n",
    "    \n",
    "    written by Jeff Shragge, jshragge@mines.edu, 10/2019\n",
    "    '''\n",
    "    n = len(U)\n",
    "    \n",
    "    AA = K *dt/(dx*dx) ## Diffusivity: alpha\n",
    "    CC = vx*dt/ dx     ## Courant :C\n",
    "    \n",
    "    ## . . Set up coefficients\n",
    "    aa =-(CC+2*AA)*np.ones(n)\n",
    "    bb =  4*(1+AA)*np.ones(n)\n",
    "    cc = (CC-2*AA)*np.ones(n)\n",
    "    \n",
    "    ## . . Treat ends\n",
    "    aa[0  ]=0\n",
    "    cc[n-1]=0\n",
    "    \n",
    "    ## . . Set up known solution matrrix\n",
    "    dd = np.zeros((nx))\n",
    "    dd[1:nx-2] = -(CC-2*AA)*U[2:nx-1]+\\\n",
    "                   4*(1-AA)*U[1:nx-2]+\\\n",
    "                 +(CC+2*AA)*U[0:nx-3]\n",
    "    dd[1   ]+=+(CC+2*AA)*U[0   ]\n",
    "    dd[nx-1]+=-(CC-2*AA)*U[nx-1]\n",
    "    U = TDMAsolver(aa, bb, cc, dd)\n",
    "    \n",
    "    return U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def convection_diffusion_1D_animate(i):\n",
    "    global k\n",
    "    ax1.clear()\n",
    "    ax1.plot(xx,c[:,k],'g',xx,d[:,k],'b',linewidth=3)\n",
    "    ax1.legend(['Fast Diff','Slow Diff'])\n",
    "    ax1.grid(True)\n",
    "    ax1.set_xlim([xmin,xmax])\n",
    "    ax1.set_ylim([-0.1*T0,1.1*T0])\n",
    "    ax1.set(xlabel='Distance (m)',ylabel='Temperature ($^\\circ$C)')\n",
    "    ax1.set_title('Heat Flow at %f s'%(k*dt))\n",
    "    k += kskip    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "diffusion #: alpha (fast) = 100.0\n",
      "diffusion #: alpha (slow) = 1.0\n",
      "Courant   #: C = 0.2\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_7.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Time axis parameters\n",
    "nt,dt=600,1\n",
    "\n",
    "## . . Space axis parameters\n",
    "xmin,xmax=0,1\n",
    "nx = 101\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "dx = xx[1]-xx[0]\n",
    "\n",
    "## . . Initial Condition for copper and iron bars\n",
    "sigma,midpoint,T0=0.05,(xmax-xmin)/4,100\n",
    "Uc = np.zeros((nx))\n",
    "Ui = np.zeros((nx))\n",
    "Uc = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "Ui = T0*np.exp(-(xx-midpoint)**2/(2*sigma**2))\n",
    "\n",
    "## . . FD parameters\n",
    "Kf = 0.01 ## . . thermal fast\n",
    "Ks = 0.0001 ## . . thermal slow\n",
    "vx = 0.002\n",
    "LB,RB = 0,0 ## . . Derivative Boundary conditions on left and right\n",
    "\n",
    "print('diffusion #: alpha (fast) =',Kf*dt/(dx*dx))\n",
    "print('diffusion #: alpha (slow) =',Ks*dt/(dx*dx))\n",
    "print('Courant   #: C =',vx*dt/dx)\n",
    "\n",
    "## . . Iterate over solution\n",
    "c = np.zeros((nx,nt))\n",
    "d = np.zeros((nx,nt))\n",
    "\n",
    "for i in range(nt):\n",
    "    ## . . run for copper bar\n",
    "    Uc = Solve_Convection_Diffusion_homogeneous_Thomas(Uc,Kc,vx,dx,dt,LB,RB)\n",
    "    c[:,i]=Uc\n",
    "    ## . . run for iron bar\n",
    "    Ui = Solve_Convection_Diffusion_homogeneous_Thomas(Ui,Ki,vx,dx,dt,LB,RB)\n",
    "    d[:,i]=Ui\n",
    "    \n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig7,ax1 = plt.subplots(1)\n",
    "fig7.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig7.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim7 = animation.FuncAnimation(fig7,convection_diffusion_1D_animate,frames=int((nt-2*kskip)/kskip),interval=100)\n",
    "anim7.save('./movies/Ex7_7.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_7.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 8. Solution of the 1D convection-diffusion equation for both slow (blue) and fast (green) diffusion coefficients.**\n",
    "\n",
    "## Solving the 2D Convection-Diffusion Equation\n",
    "\n",
    "Let's now look at expanding our dimensionality to cover the 2D convection-diffusion equation.  This equation is given in general by the following:\n",
    "\n",
    "$$\\frac{\\partial \\phi}{\\partial t} + {\\bf v} \\cdot \\nabla \\phi = \\kappa \\nabla^2 \\phi, \\tag{37a}$$\n",
    "\n",
    "and with specific Cartesian partial derivatives\n",
    "\n",
    "$$\\frac{\\partial \\phi}{\\partial t} + v_x \\frac{\\partial \\phi}{\\partial x} + v_y \\frac{\\partial \\phi}{\\partial y}= \\kappa \\left(\\frac{\\partial^2 \\phi}{\\partial x^2}+\\frac{\\partial^2 \\phi}{\\partial y^2} \\right).\\tag{37b}$$\n",
    "\n",
    "Given that we (1) know how to solve the 1D convection-diffusion equation from above and (2) used an alternating direction implict (ADI) solver in Module 8, we might expect that their combination presents us with a two-step procedure to solve the 2D convection-diffusion equation\n",
    "\n",
    "### Step 1 - First half-step in x direction\n",
    "\n",
    "The first half step provides the following update in the $x$ direction:\n",
    "\n",
    "$$\\frac{\\phi^{n+1/2}_{i,j}-\\phi^n_{i,j}}{\\Delta t/2} + v_x \\left(\\frac{\\phi^{n+1/2}_{i+1,j}-\\phi^{n+1/2}_{i-1,j}}{2\\Delta x}\\right) + v_y \\left(\\frac{\\phi^{n}_{i,j+1}-\\phi^{n}_{i,j-1}}{2\\Delta y} \\right) \n",
    "=\n",
    "\\kappa\n",
    "\\left(\n",
    "\\frac{\\phi^{n+1/2}_{i+1,j}-2\\phi^{n+1/2}_{i,j}+\\phi^{n+1/2}_{i-1,j}}{\\Delta x^2} + \n",
    "\\frac{\\phi^{n}_{i,j+1}-2\\phi^{n}_{i,j}+\\phi^{n}_{i,j-1}}{\\Delta y^2}\n",
    "\\right)\n",
    "\\tag{38}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplying through by $2\\Delta t$ and then using $C_x = \\frac{v_x\\Delta t}{\\Delta x}$, $C_y = \\frac{v_y\\Delta t}{\\Delta y}$, $\\alpha_x = \\frac{\\kappa \\Delta t}{\\Delta x^2}$, and $\\alpha_y = \\frac{\\kappa \\Delta t}{\\Delta y^2}$ leads to\n",
    "\n",
    "$$4\\phi^{n+1/2}_{i,j}-4\\phi^n_{i,j} \n",
    "+ C_x\\left(\\phi^{n+1/2}_{i+1,j} - \\phi^{n+1/2}_{i-1,j}\\right) \n",
    "+C_y\\left(\\phi^{n}_{i,j+1}-\\phi^{n}_{i,j-1} \\right) \n",
    "=\n",
    "2\\alpha_x \\left(\\phi^{n+1/2}_{i+1,j}-2\\phi^{n+1/2}_{i,j}+\\phi^{n+1/2}_{i-1,j}\\right) +\n",
    "2\\alpha_y \\left(\\phi^{n}_{i,j+1}-2\\phi^{n}_{i,j}+\\phi^{n}_{i,j-1}\\right)\n",
    "\\tag{39}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collecting like terms and moving those dependent on time step $n+1/2$ to the left-hand side and $n$ to the right-hand side results in\n",
    "\n",
    "$$\n",
    "(C_x-2\\alpha_x)\\phi^{n+1/2}_{i+1,j}\n",
    "+4(1+\\alpha_x)\\phi^{n+1/2}_{i,j} \n",
    "-(C_x+2\\alpha_x)\\phi^{n+1/2}_{i-1,j} \n",
    "=\n",
    "-(C_y-2\\alpha_y)\\phi^{n}_{i,j+1}\n",
    "+4(1-\\alpha_y)\\phi^{n}_{i,j}\n",
    "+(C_y+2\\alpha_y)\\phi^{n}_{i,j-1} \n",
    "\\tag{40}$$\n",
    "\n",
    "which is a tridiagonal matrix as is expected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2 - Second half-step in y direction\n",
    "\n",
    "The second half step provides the following update in the y direction:\n",
    "\n",
    "$$\\frac{\\phi^{n+1}_{i,j}-\\phi^{n+1/2}_{i,j}}{\\Delta t/2} \n",
    "+ v_x \\left(\\frac{\\phi^{n+1/2}_{i+1,j}-\\phi^{n+1/2}_{i-1,j}}{2\\Delta x}\\right) \n",
    "+ v_y \\left(\\frac{\\phi^{n+1  }_{i,j+1}-\\phi^{n+1  }_{i,j-1}}{2\\Delta y}\\right) \n",
    "=\n",
    "\\kappa\n",
    "\\left(\n",
    "\\frac{\\phi^{n+1/2}_{i+1,j}-2\\phi^{n+1/2}_{i,j}+\\phi^{n+1/2}_{i-1,j}}{\\Delta x^2} + \n",
    "\\frac{\\phi^{n+1  }_{i,j+1}-2\\phi^{n+1  }_{i,j}+\\phi^{n+1  }_{i,j-1}}{\\Delta y^2}\n",
    "\\right)\n",
    "\\tag{41}$$\n",
    "\n",
    "Multiplying through by $2\\Delta t$ and again using $C_x = \\frac{v_x\\Delta t}{\\Delta x}$, $C_y = \\frac{v_y\\Delta t}{\\Delta y}$, $\\alpha_x = \\frac{\\kappa \\Delta t}{\\Delta x^2}$, and $\\alpha_y = \\frac{\\kappa \\Delta t}{\\Delta y^2}$ leads to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$4\\phi^{n+1}_{i,j}-4\\phi^{n+1/2}_{i,j} \n",
    "+ C_x\\left(\\phi^{n+1/2}_{i+1,j} - \\phi^{n+1/2}_{i-1,j}\\right) \n",
    "+C_y\\left(\\phi^{n+1}_{i,j+1}-\\phi^{n+1}_{i,j-1} \\right) \n",
    "=\n",
    "2\\alpha_x \\left(\\phi^{n+1/2}_{i+1,j}-2\\phi^{n+1/2}_{i,j}+\\phi^{n+1/2}_{i-1,j}\\right) +\n",
    "2\\alpha_y \\left(\\phi^{n+1  }_{i,j+1}-2\\phi^{n+1  }_{i,j}+\\phi^{n+1  }_{i,j-1}\\right)\n",
    "\\tag{42}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Collecting like terms and moving those dependent on time step $n$ to the left-hand side and $n+1/2$ to the right-hand side results in\n",
    "\n",
    "$$\n",
    "(C_y-2\\alpha_y)\\phi^{n+1}_{i,j+1}\n",
    "+4(1+\\alpha_y)\\phi^{n+1}_{i,j} \n",
    "-(C_y+2\\alpha_y)\\phi^{n+1}_{i,j-1} \n",
    "=\n",
    "-(C_x-2\\alpha_x)\\phi^{n+1/2}_{i+1,j}\n",
    "+4(1-\\alpha_x)\\phi^{n+1/2}_{i,j}\n",
    "+(C_x+2\\alpha_x)\\phi^{n+1/2}_{i-1,j} \n",
    "\\tag{43}$$\n",
    "\n",
    "which, again, is a tridiagonal matrix as is expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Solve_2D_Convection_Diffusion_homogeneous_Thomas(U,K,vx,vy,dx,dy,dt,LB,RB,BB,TB):\n",
    "    '''Set up tridiagonal matrix and call Thomas Algorithm \n",
    "    usage: x = Solve_2D_Convection_Diffusion_homogeneous_Thomas(U,K,vx,vy,dx,dy,dt,LB,RB,BB,TB):\n",
    "    input: \n",
    "        U: heat solution at time step n   (nx)\n",
    "        K : thermal diffusivity\n",
    "        vx: convection velocity in x-direction\n",
    "        vy: convection velocity in y-direction\n",
    "        dx: spatial sampling in x-direction\n",
    "        dy: spatial sampling in y-direction\n",
    "        dt: temporal sampling\n",
    "        LB: Left   boundary condition (Dirchelet)\n",
    "        RB: Right  boundary condition (Dirchelet)\n",
    "        BB: Bottom boundary condition (Dirchelet)\n",
    "        TB: Top    boundary condition (Dirchelet)        \n",
    "    output:\n",
    "        u: heat solution at time step n+1 (nx,ny)\n",
    "    depends on:\n",
    "        TDMAsolver\n",
    "    \n",
    "    written by Jeff Shragge, jshragge@mines.edu, 10/2019\n",
    "    '''\n",
    "    ## . . Get dimensions\n",
    "    nx,ny = np.size(U,0),np.size(U,1)\n",
    "    U1 = np.zeros((nx,ny))\n",
    "    \n",
    "    ## . . Define diffusivity and Courant numbers\n",
    "    AX = K *dt/(dx*dx) ## Diffusivity x: alpha\n",
    "    AY = K *dt/(dy*dy) ## Diffusivity y: alpha\n",
    "    CX = vx*dt/ dx     ## Courant x : C\n",
    "    CY = vy*dt/ dy     ## Courant y : C\n",
    "    \n",
    "    ## . . Set up coefficients in X\n",
    "    ax =-(CX+2*AX)*np.ones(nx)\n",
    "    bx =  4*(1+AX)*np.ones(nx)\n",
    "    cx = (CX-2*AX)*np.ones(nx)\n",
    "    \n",
    "    ## . . Set up coefficients in X\n",
    "    ay =-(CY+2*AY)*np.ones(ny)\n",
    "    by =  4*(1+AY)*np.ones(ny)\n",
    "    cy = (CY-2*AY)*np.ones(ny)    \n",
    "    \n",
    "    ## . . Treat ends\n",
    "    ax[0]=cx[nx-1]=ay[0]=cy[ny-1]=0\n",
    "    \n",
    "    ## . . Set up known Solution matrix\n",
    "    dx = np.zeros((nx))\n",
    "    dy = np.zeros((ny))\n",
    "    \n",
    "    ## . . Solve first for update in x direction\n",
    "    for iy in range(1,ny-1):\n",
    "        dx[1:nx-2] = -(CY-2*AY)*U[1:nx-2,iy+1]+\\\n",
    "                       4*(1-AY)*U[1:nx-2,iy  ]+\\\n",
    "                     +(CY+2*AY)*U[1:nx-2,iy-1]\n",
    "        \n",
    "        d[0   ]+=+(CY+2*AY)*U[0   ,iy]\n",
    "        d[nx-1]+=-(CY-2*AY)*U[nx-1,iy]\n",
    "        \n",
    "        U1[:,iy] = TDMAsolver(ax, bx, cx, dx)\n",
    "    \n",
    "    ## . . Solve second for update in yu direction\n",
    "    for ix in range(1,nx-1):\n",
    "        dy[1:ny-2] = -(CX-2*AX)*U1[ix+1,1:ny-2]+\\\n",
    "                       4*(1-AX)*U1[ix  ,1:ny-2]+\\\n",
    "                     +(CX+2*AX)*U1[ix-1,1:ny-2]\n",
    "        \n",
    "        d[0   ]+=+(CX+2*AX)*U[ix,0   ]\n",
    "        d[ny-1]+=-(CX-2*AX)*U[ix,ny-1]\n",
    "        \n",
    "        U[ix,:] = TDMAsolver(ay, by, cy, dy)\n",
    "        \n",
    "    return U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAAEgCAYAAACAZQdyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnXmUZVV97z9fm7EFGW2GBqQRRNEEoghqjLbhJYC+iEkwTkviiKwEh6yYgD6ToCaKMb6gQSWEEMKKCb4XJ3xByGAUjYK0pmUSsIEWWsYGmkHGht/7Y+9ddWrXubfurbrnVlX397NWrV3nnL33+Z1z7/2d7/ntSRGBMcZ0xZPm2wBjzKaNnYwxplPsZIwxnWInY4zpFDsZY0yn2MkYYzrFTmYTRNJaSWuHLHOKpJC0co7nPifXs+9c6lno9Lpfed83hqxrk75nm7WTkbSLpLdJ+pKkNZIeknSvpG9LequkafdH0sr8hWj+PSjpVkkXS/q4pF+YhS3lS3vKSC5uev1vyvW/qYv654qkJ0t6j6SvS7pD0qOSNkj6nqQ/k7TffNs4W0blwBcrW8y3AfPMq4HPArcC/wncBOwG/AZwFnC0pFdHe4/FnwDn5P+3Ap4KPBd4L/BeSf8IvCMiHuj0Cto5YhZlTgfOI92DsSLpBcA/A8uBdcAFwC3Ak4FfAE4C/kDSCyLiB+O2b0ieBTw4ZJn3AacCPx29OfPP5u5krgNeCfxLRDxRdkp6P/A94DdJDucLLWXXRsQp9U5JhwDnAq8HdgaOHr3Z/YmI62dRZj2wvgNz+iLpmcBFwHbAycAnImJjlWcF8DHgKeO2b1gi4ppZlLmV9KDbNIkI/7X8Ae8HAvirav/KvP8bfcouA+7I+V414PlOyflPqfafk/fvC7wDuAJ4GLgdOBPYoaWutSQnWLa/keto+9u3Ov/Kqq5XAf9Acsg/Ax4Avg+8C3hSy7nPadY7wHX/W87/kQHybl1tH0By6D8FHiWpn3OBA/rc35XAsaSHyIPA3SQFt7zHOZ8HXAjcD9wH/Dvwwj73a8p3I38Wrfd+kHsG/BZwMXAv8FD+/N9X34vm5w4sBT5OUqWPAGtIalDz8Vva3JVMPx7L6ca+uVqIiDsk/TXwAeANwJdHYM+fA0cCXwX+FXgZ8HZgf+CXZyh7DrABOAb4CrC6cWzDDGVPBZ4ALiX9mHfI5/sk8HzgjUNcwxSyQvkfJKf55zPlj4hHGmWfT/rBbw+cD1wNPJN0v4+RdERErGqp5ndI6vV84JvA4cBrgIMlHVKd40X5HFsBXyT9WA8hOe2vD3iZp5Ec9UuBvyc5gYGQ9BGSQ1kP/CPJwR8NfAQ4UtKvRMRjVbEtSd+PPYGvkb6/ryJ9jtsAHxz0/CNjPjzbQv8jvUZeQXq6HFkdW8kMSibnOyLn+8mA5zyF/krmJmCfysaL87HDqjJraSiZvO9NOe+bZjj/ymr/01vyPon0gwng8B727jvANb8x5/32kJ+PgB/lsm+ojr0m77+GhtJqXN99wM9VZf4xH/ut6hzX5P3HVPnfzaQiqe/XtO9Gr3vb756R1FL53HevPvev5mPvb/ncgxTT2raxfxnpYbIB2HLcv6fNunWpD6cCzwEuiIiLZllHCeI9dTQm8aGImAjKRopb/F3ePGxE55hGtMR3IsWvPpk3j5xD9XvkdN2Q5V5EUi3fjYjPVbZ9Hvg2cCDw4payn4qIK6p9f5PT5n18Ua7j4oj4SpX/dGDouNeQvCWnfxoRt5Wd+XP/fZK6fFuPsu+KiIcaZe4gKdgdSNc0VuxkKiS9i/QhXsMcXgVIT0JIT5ZR0Cb9b87pTiM6xzRyM/+pki6X9EBptifFZSC1CM26+pwOe4+em9Neryxlf1tXgkHvYznHN+vMEfE4yZF1Sc9rjIjrSI55haQdq8P3RsSalvo6/670wjGZBpJ+l/SEvho4IiLunkN1e+b0zjkblmiLnZR40ZIRnWMK+Qt8GbCCFCg9lxQo3QjsSHpt2HoOp7glp3sNWW6HnPZqkSn76x8gDH4fyzlu73GO23rsHxWDXOM+OV/zmnrF2Dr9rvTDTiYj6T3AXwJXkhzMHXOs8mU5vXSO9cwnbyM5mA9G1Vwv6YUkJzMXiho4VNIOEXHvgOVKvt17HN+jyjcbStndehzvde5R0bzGtlezUVzjWPDrEiDpJJKDWQ28bK4ORtIyUnMzwOf65R0jj+d0mCfZ/jlt6yf00rmZAxFxI6n1ZhvgD2bKL6mopv/O6coeWcv+uXTcK2WnXaekJbTHe3oxm3vf8xol7U9SfzdGxEytg/POZu9kJP0RKdD7fZKCmVOHNEkHk/p+7EoKHJ8/dytHwl053WeIMmtzurK5Mw+beN/cTQLgnaQWn/dJ+n1J09S1pH0knUdqcQH4L+Ba4MWSjq3yHgu8hNSvZy5xk+/kc7xE0jHVsROBpw9R12zu/dk5/YCkicaD7OD+gvTb/dsh6ps3NuvXJUm/DXyI9KT5FvAuSXW2tRFxTkvxfRvjjLYkOZXn5T9ICuYdLeXmi++SOp+9R9LOTMYa/qrPa8q5JIVxmqSXAT8mdYD7n6R+I6+Zq1ERcY2kI0lq6S+Ad0v6DyaHFRwM/CIpOPyxXCbyZ/dvwOclfYUUqD+Q1CfkfuC4aPTinoVdIemt+RxfkFT6yRxM6ttzIXDUgNX9J6k16KOSngPck8/xp33O/x1Jfw78IXClpH8mdYY8mtTy+W1Sh7sFz2btZEjxBkgy9j098nyTyTFKTZ4G/En+/2FSwO3HpB/K5yJidUuZeSMi7pH0mySb30z6AUPqzdvqZCLiFkm/RFJ6LyY1V19D6tD274zAyeTzXCLpQFLnwlcCryC1gjxI+mF/Ajgzv16VMpfmDnkfIP3of43Uae2fgA9HxLUjsOu/8vX/GZPDQy4lKbsjGdDJRMSPslN8L+nebZMP9XQyudxJkv6bpJyOIz3Mridd8yci4tGhLmieUO6sY4wxnbDZx2SMMd0yVicj6ew8V8iVPY5L0qfy3C6XS3puWz5jzOJh3ErmHPq/xx5NCiweABxPmuvFGLOIGauTiYiLST1Ge3EMcG4kLgF2lLRHn/zGmAXOQmtdWs7kGAtI4zOW09K1WtLxJLUDbPm81IJsjOmOW9dHxNADfheak5nWSYUeg+ci4kzSpE1Ie8aEvzHGdMQHfzKbUgutdWkdsHdjey8mB9EZYxYhC83JnA8cl1uZXkAatr7pzn1qzGbAWF+XJP0TqbfkrpLWkXqfbgkQEWeQZvR6OamX54OknqnGmEXMWJ1MRLxuhuMB/O6YzDHGjIGF9rpkjNnEsJMxxnSKnYwxplPsZIwxnWInY4zpFDsZY0yn2MkYYzrFTsYY0yl2MsaYTrGTMcZ0ip2MMaZT7GSMMZ1iJ2OM6RQ7GWNMp9jJGGM6xU7GGNMpdjLGmE6xkzHGdIqdjDGmU+xkjDGdYidjjOkUOxljTKfYyRhjOsVOxhjTKXYyxphOsZMxxnSKnYwxplPsZIwxnWInY4zpFDsZY0yn2MkYYzrFTsYY0yl2MsaYTrGTMcZ0ip2MMaZTxu5kJB0l6VpJaySd3HJ8B0lflfRDSVdJevO4bTTGjI6xOhlJS4BPA0cDBwGvk3RQle13gasj4mBgJfAJSVuN005jzOgYt5I5DFgTETdExKPAecAxVZ4AtpckYDvgbmDjeM00xoyKcTuZ5cDNje11eV+T04FnAbcAVwDvjogn6ookHS9plaRV8GBX9hpj5si4nYxa9kW1fSSwGtgTOAQ4XdJTphWKODMiDo2IQ2Hp6C01xoyEcTuZdcDeje29SIqlyZuBL0ZiDXAj8Mwx2WeMGTHjdjKXAQdIWpGDua8Fzq/y3AQcASBpN+BA4IaxWmmMGRlbjPNkEbFR0onARcAS4OyIuErSCfn4GcCHgXMkXUF6vTopItaP005jzOgYq5MBiIgLgAuqfWc0/r8F+NVx22WM6Qb3+DXGdIqdjDGmU+xkjDGdYidjjOkUOxljTKfYyRhjOsVOxhjTKXYyxphOGXtnPLOpsmXH9T/Wcf2mK6xkjDGdYiVjhmQuiqXX122QOcl6ndcKZ6FjJWOM6RQ7GWNMp/h1yfRhkFejYb5CowgO169WbXX6FWohYSVjjOkUKxnToJfS6Pc1mU2ZQfLOZoGKUqa2ycpmPrGSMcZ0ipWMYfqTv/5atKmVmfLM9atVK5miRrbocbxfHcU2K5r5wErGGNMpVjKbJYPGUbZs2d+2r21/v5akQdTIY1We+nwbq3xtZevzOFYzH1jJGGM6xUpms2GYuEqtRvopmbK9bZ8y/WyAdkVRK5WHqu02hVPv66WU2pSNVU1XWMkYYzrFSmaTZy4KZttqf3Pf0upYrWS2ZZJBv2ZN5VErl4eqPA9W++vybfRroXILVFdYyRhjOsVOxhjTKX5d2mQZpHNc/Zq0bZWW409plKnzlHT7qkzzdalfMBjam6Pr16X7q/3bVtsA91XneaglT7POtqCxX5tGjZWMMaZTrGQ2efopmDqwW6uTp1QpTCqW+ljZr+mn26Y6TR13LenDLaYSOS1K5r4qvb+Rd4vqWE2voQptRplRYSVjjOkUK5lNjl4d3pr7a0nRS7kUdbJzo+xTqn1ZuezI1HS7RpHZKJkHcroh17+htqlWXc2Ky7XWiuahKl+TOgbj2MyosJIxxnSKlcwmQ6/WpLaWnZliLzv3SGGiE96utKe1ooFJVbMN7RQF80Bj34YqXV/SrGzW75J3tLViDTrNZ93q1MStTaPCSsYY0yljVzKSjgI+CSwBzoqIU1vyrAROIz1G1kfES8dq5KKm/kjbpl7opWRq5bI8V9FoMdorp7v32C5pU8mU/+vYTB2L2cAk5f/bqnRdTos6Wrd0oggbm6qmjVqNtLUkzWZyLNOPsToZSUuATwO/Qvq6XCbp/Ii4upFnR+AzwFERcZOkZeO00RgzWsatZA4D1kTEDQCSzgOOAa5u5Hk98MWIuAkgIu4Ys42LjJlak/r1wO3VipQVzDZZwRS1ArDvDGmtbIAn7fgzAJZulwY1LtnicQAe37gEgAcfSGrkiQ1PnixUK5e1Oa1br5rf4HXZ3oeXM5XHqrSokzbl02syLMdmZsu4YzLLgZsb2+uY+EZP8AxgJ0nfkPR9Sce1VSTpeEmrJK2aHJFrjFlojFvJtHQHnejSWdgCeB5wBOlR811Jl0TEdVMKRZwJnAkg7VnXsRnSqzWpTmHG1qQtKgWzf6PoM6t9E9vpCb/709IzZDcmBeiOOcCyND8MlpCVzNZZyTw5KZkNu00Gcm4/ML0l3/aTvdOO3bcslSXalExhbbZ/Y7muWsG0qZW2+9TEsZnZMm4nsw7Yu7G9F3BLS571EfEz4GeSLgYOBq7DGLPoGNjJSNoaeCHwAmBPkspYD1wLXFziLDNwGXCApBXAT4HXkmIwTb4CnC5pC2Ar4HDgLwe10xizsJjRyUjaH3gP8AZgB+AJ4F5ST6adSQ2TIen7pFahcyPiiba6ImKjpBOBi0hN2GdHxFWSTsjHz4iIH0m6ELg8n+usiLhyjte5iTFIZ7N+gd96eoby2pSbg8tr0r45La9EAM+Zmu5wSIrQPn3rNQCsyBHaPRoCdVfuymdLgxknXpdIr0v3ZzvWs8tEmVvZE4Abn5aMuH739H5273Y5otzWsa9+G1pbmrfL9dVTPzTvSV241+uR5wUelr5ORtLpwNuB1cCHgIuBH0bExkae3Ujq5hXA/wZOkvSmiLi0rc6IuAC4oNp3RrX9ceDjQ1+NMWbBMZOS2Qs4PCJW98oQEbeTXnG+IumdwDtIMZRWJ2O6plfAt9FprVcnvDI0oDQ/75vTZuA3K5inHn4TAAfl3gcl3Z81uejaiSLLuB2AnXIAuFYy9+Ro7h3sNlFmbT75shxA3n7rNObg6sMPAuBO9kkZ2wZV1un6EgAuCqZMD9G8J/XgyZJarcyVvk4mIl41TGUR8QjwqTlZZIzZpPAAyUXPTCs29hsgWZRMbvItSqZPTKbEYIpy+QX+G4BDWD1l/9OzogHY+eYsN+7KOx7J6dY5zaGYu/eeDLRcn+XTLrnQ1hOFEqsP2QqAex9o9PqrB1VODK4sPSfK9d6d07b1pOpBk/V9dVP2sAzlZHKLzwtJzdDTQm8RcfaI7DLGbCIM04T9XOBLpOdcr051djKdMpdWpWZLyvZT09LBrWdMZjIuUVqRimIpCuZQVgHw7DtyT4ZmFO/GnN6e01rJ5FDMzismAyw7H5IaFJcum9qb+5Fc6P6tU2+8H+w/2SLFbfma60GV5fo2VNc9ZVKremLymVqZwEMNBmMYJXMGKZT2KuAa4NFOLDLGbFIM42QOAn4rN0GbBU8dq2mqoPJ/j6kzq2kcylABmOwHU1qRiqKZUDD/lTOuapzuxzktIw1+ltMyHrKMsz+gUSbnefYvpnofXJZagu7KAZw7cqFbnrbnRJHb1u431f76uspUnn3vicOUo2aYAZLXMfm1MMaYgRjGbb8f+JikS8s0DGYhUffv6LVwW/0/k4MNeyia5mDH0pO39IOZaEUqMZhVVQrwo5TclTsBP/R4tiJ1k2GXIkbuZTr5sfb0X03nuTkPfbsxB4yatt22+37t19Gc1DyduUph5gGmjrvMloGdTERcmGes+7Gk64B7pmfxDHbGmKkM07p0MvCHwJ2ksPzjXRlljNl0GOZ16T3AXwMnRoQdzKKjpeNZ6elU0uq1qcxot2Nj8t0y2LEMFZjoaFeaqUuQ90eTZ7s6x41/krdLo/TS/C0qceWDmubukNM8eqCcZ9net0+xo2lbsfeJHZ889Xrq63zYQd5xMkzgdynwf+1gjDHDMIwr/xqpt+/XO7LFDExbd/hetB3fYmo1Ja2e+GVO3qWN6U3LdA1lsOPEUIHS0S7HYe9qTEVWFMzanJYucM0VtgF2a5TZpcRzS735PDvtvWGKHU3bir0PbPPkqdczbdK7fuuD96LtuIcYDMIwTuY04BxJABcyPfDLgBNXGWM2I4ZxMqWb1YdJc8u0sWRu5pjRMcBH22ta27xdVhVY0ojxl/8n9pUhAiXNnegearxUF61RFEw9BLEcb5aZ6LBX1V+ff4pt2d6eS08N9G13nGbUDHNH38L0Sb+NMaYvw/STOadDO8zIGSBe0GsMYN4u6yKVVQVgcqKpkk4McixpDods29C0pRWpjsFUk35OKTPRt7yqvz7/4w3xXOztdT2DhVAcZxk1XgvbGNMpM83x+0ngoxFxW798VZnfALaKiPPmapzpRdvTttdH2dYdfuPUauo1qXNaVnYs6yLB5KTfZcrMibm/y8yZebDjLpPjFnlaczk/Gv1kynGml5kYNFnqzecp5y12PNiYQrPYW1/HdCXTbw3sXtRrN5lBmUnJ7AfcIOnzkl4paZc6g6QnSTpE0h9Juhb4LJNTjxljNnNmmuP31yS9BHgv8AXgSZJuIQ0teATYiclZ8m4F/gb4y4i4r0eVZt5oPoHzU7l+4peJt3MXmLI2dXNlx/UTUy0kiVGmzJyYcKpM19AY7Fh68u420wDJZzVMLPWsYMp5ynmLHRuYtG1iLe3SCbhcT32dViVjZcbAb0RcDFwsaQ/gV0mLre1Jciw3AJ8jLZXyrV7rLRljNl+GaV26Ffj7/GeMMQPhnkebDL0Wky/bzS5wVXe46jVpIs3h/tsPXDaRtazsWNZFKqsKlDl5JzrRNcmDHXcZZma8Q3N6CFPOU85b7LidSdsm5vStr+MBKuqVJKH3ffMr1VxxE7YxplOsZDZZ6idzs4m2/J87cJe5bysFw7q8+ZO9J0qWtanLyo5lXaSyqkCZk3fKRK15uoaZVisoQV5gQsFctSzNdnd1Dh+vyYqmzIzXtK3Y21PRTHRY73dPrFxGjZWMMaZTrGQWFeVp229agpKnfLRt8Yf7p6Ybcgf/9Xl3UQJrc7r75Pmu3z0pibI2db2yY1lVoMzJC3NbQbIomNVZ2pTt6x/JC3SvadyLYm+xv1zPhJKprrs1TtWmcnrheX8Hoa+SkfTOcRlijNk0mUnJnCbp1cBbImLNDHnNvFDHEOqPtC3WUJ7apc9kXlGxrBldYhv1rP/AvdulJQyuPnzKZJkTKzuWdZHKqgIwOWVmmXCqTM9QBjeWoQJ3TARnJluR1lSKpqT3rs5LKVzTMGItU+0vSmYiFlNPNtHSQXEajtXMlZliMkcAewA/lPRe5RmrjDFmUGYaVvANST8HfBD4CHCspLdExNVjsc7Mgrp/R0mba0qX9Ybum7q9PgdHygTcdQoTU1remZuMVh+yFTC5NnVZ2fHGiYW0Jyf9LlNm1kqmDHZcz+TQuNIPptRTYjATCiZ3y6Gpr9fmtI7JTAylq5VM8570um9mrgwyrOBh4CRJnwfOAn4g6ctMXws7IuK3O7DRGLOIGaZ1aQ1pncBDgF+ixcmMyigzE834Qa+Wpl6tTM3/SytLtaLiuqVTi7Z9S3KD0b0PJGXxg/2TCilrUzdXdizLlpRJv2slU6ZraA52LD15J/rBlFakEoMpCqZfTGbaxJ91q1LzngzaquQWpWEZqJ+MpFcCVwOvAt4eEcsjYkX1t9+AdR0l6VpJa/KCcb3yPV/S45KOHehKjDELkpmasJ+aX5O+BPwAeE5E/O1sTyZpCfBp4GjSDACvk3RQj3wfAy6a7bmMMQuDmV6XrgGeAI6LiM+N4HyHAWvK0imSzgOOIamkJu8kzV/z/BGcczOhV6CyLZBZXh/q96H8SrKxvDa1NCaWaupBlbelsretTYL2tt0nhW1Z2bGsi1RWFShz8pYZ7Sbmg4FpQxsmXoV6pc28G8ub+91Vel+VNu/JTANMzWyZ6XXp68CzR+RgAJYDzckY1+V9E0haDvw6cEa/iiQdL2mVpFVTWwmMMQuJmZqwXz3i87X1s6kDxqcBJ0XE4/265UTEmcCZANKem3HQuddQgzrw26TsK0/0Las083D2/2sbn0NPJZPT3MLc7MBX1qaeWNmxnL6eW3hyWevp9dbKpt6GhoL5aU6LgikB33K95d4ME/h1wHe2jHvs0jpodAWFvYBbqjyHAudlB7Mr8HJJGyPiy+Mx0RgzSsbtZC4DDpC0gvS4eS3w+maGiJgY8C/pHOD/2cEMQ5EHtVxoUq/hOMOUzBt3nvx/bW7e7qVkWoYiTHTm24Z26jmGm/WWdH2PdMqrch2D6dUJrx4MCb0XazJzZaxOJiI2SjqR1Gq0BDg7Iq6SdEI+3jcOY4xZfIx9qoeIuAC4oNrX6lwi4k3jsGnToI7N9HsS10pmpjqbT/wyLURWN2VQZa1gWoYiTGvM6rXuE/SeEpS65aipwnrFYHopmbYBkvV9cyxmrnjSKmNMp3jSqk2OmVqbmrRNMt7M22+aiJJmZbMhTxexoaVFcDZKZoKiXHqpk/sbeQdVLnVfmDasYEaFlYwxplOsZDZ56tYmmP6U7vU1aOstXA+uLHGQrGQmFNS2EyV4uOyb6TxNu+oWoF6DG5sqrFc/mFrJ9FN1ZtRYyRhjOsVKZpOlX2vTFlWeQv3E37YlX9m3tKqr7j3cUDIDf83aFFOtSupJuNoWaKvL9lIwberFsZhRYyVjjOkUOxljTKf4dWmTp61Juw4Gz/SK0Pya1K8k9etRvyDvTLP4NalfcQZ5Beo1TYNfk+YTKxljTKdYyWw2tM0L3KvZtlYCbUqm1wqV1QRYrfQbvNnLhl5qpJ+Sqff3OofpEisZY0ynWMlsltRP8FrZ9IvVtKmb5n56HB+WmQYqDhJX8WDHhYCVjDGmU6xkDMNNEzFomVErmfq8oypjusZKxhjTKVYypsFMsZq51NWsb5C8vRhGZZmFgJWMMaZTrGRMHwZRI6OYImEudVi1LHSsZIwxnWInY4zpFL8umSHp9XrSbxhBwa9FmyNWMsaYTrGSMSPCSsO0YyVjjOkUOxljTKfYyRhjOsVOxhjTKXYyxphOsZMxxnSKnYwxplPsZIwxnWInY4zpFDsZY0ynjN3JSDpK0rWS1kg6ueX4GyRdnv++I+ngcdtojBkdY3UykpYAnwaOBg4CXifpoCrbjcBLI+LngQ8DZ47TRmPMaBm3kjkMWBMRN0TEo8B5wDHNDBHxnYi4J29eAuw1ZhuNMSNk3E5mOXBzY3td3teLtwJfazsg6XhJqyStggdHaKIxZpSMe6oHteyL1ozSy0hO5sVtxyPiTPKrlLRnax3GmPln3E5mHbB3Y3sv4JY6k6SfB84Cjo6Iu8ZkmzGmA8b9unQZcICkFZK2Al4LnN/MIGkf4IvAGyPiujHbZ4wZMWNVMhGxUdKJwEXAEuDsiLhK0gn5+BnAHwO7AJ+RBLAxIg4dp53GmNGhiMUfzkgxmePn2wxjNnE++P3ZPPDd49cY0yl2MsaYTrGTMcZ0ip2MMaZT7GSMMZ1iJ2OM6RQ7GWNMp9jJGGM6xU7GGNMpdjLGmE6xkzHGdIqdjDGmU+xkjDGdYidjjOkUOxljTKfYyRhjOsVOxhjTKXYyxphOsZMxxnSKnYwxplPsZIwxnWInY4zpFDsZY0yn2MkYYzrFTsYY0yl2MsaYTrGTMcZ0ip2MMaZT7GSMMZ1iJ2OM6RQ7GWNMp9jJGGM6xU7GGNMpdjLGmE6xkzHGdMrYnYykoyRdK2mNpJNbjkvSp/LxyyU9d9w2GmNGx1idjKQlwKeBo4GDgNdJOqjKdjRwQP47HvjsOG00xoyWcSuZw4A1EXFDRDwKnAccU+U5Bjg3EpcAO0raY8x2GmNGxBZjPt9y4ObG9jrg8AHyLAdubWaSdDxJ6QA8Ah+8crSmdsquwPr5NmJAFpOtsLjsXUy2Ahw4m0LjdjJq2RezyENEnAmcCSBpVUQcOnfzxsNisncx2QqLy97FZCske2dTbtyvS+uAvRvbewG3zCKPMWaRMG4ncxlwgKQVkrYCXgucX+U5HzgutzK9ALg3Im6tKzLGLA7G+roUERslnQhcBCwBzo6IqySdkI+fAVwAvBxYAzwIvHmAqs/syOSuWEz2LiZbYXHZu5hshVnaq4hp4Q5jjBkZ7vFrjOkUOxljTKcsKiezmIYkDGDrG7KNl0v6jqSD58POhj197W3ke76kxyUdO077KhtmtFWHX18lAAAFL0lEQVTSSkmrJV0l6ZvjtrGyZabvwg6Svirph9neQeKQnSDpbEl3SGrtdzar31hELIo/UqD4emA/YCvgh8BBVZ6XA18j9bV5AXDpArb1RcBO+f+j58vWQe1t5Ps6KTh/7EK1FdgRuBrYJ28vW8j3Fng/8LH8/1OBu4Gt5snelwDPBa7scXzo39hiUjKLaUjCjLZGxHci4p68eQmpP9B8Mci9BXgn8AXgjnEaVzGIra8HvhgRNwFExEK3N4DtJQnYjuRkNo7XzGxIxMX5/L0Y+je2mJxMr+EGw+YZB8Pa8VbS02G+mNFeScuBXwfOGKNdbQxyb58B7CTpG5K+L+m4sVk3nUHsPR14FqnT6RXAuyPiifGYNzRD/8bGPaxgLoxsSMIYGNgOSS8jOZkXd2pRfwax9zTgpIh4PD1w541BbN0CeB5wBLAt8F1Jl0TEdV0b18Ig9h4JrAZ+GXg68G+SvhUR93Vt3CwY+je2mJzMYhqSMJAdkn4eOAs4OiLuGpNtbQxi76HAednB7Aq8XNLGiPjyeEycYNDvwfqI+BnwM0kXAwcD8+FkBrH3zcCpkYIeayTdCDwT+N54TByK4X9j8xUQm0VAagvgBmAFkwG0Z1d5XsHUoNT3FrCt+5B6Nb9oMdzbKv85zF/gd5B7+yzgP3LepcCVwHMWsL2fBU7J/+8G/BTYdR6/D/vSO/A79G9s0SiZ6G5IwnzZ+sfALsBnsjrYGPM0IndAexcEg9gaET+SdCFwOfAEcFZEzMtUIAPe2w8D50i6gvTjPSki5mUKCEn/BKwEdpW0DvgTYMuGrUP/xjyswBjTKYupdckYswixkzHGdIqdjDGmU+xkjDGdYidjjOkUOxkzI5L+j6S7Je1e7V8iaZWkH0vadoB6fj+P3B1Jl2FJv5fr8/d4AeMPxwzCiaSu45+p9r+XNGL3bRHxUL8KJO1IGm38oRhdv4kzgGXAb4+oPtMBdjJmRiKNYv494NclvRpA0jOAU4C/johB5mt5K/AY8KUR2vUQcC7J2ZkFip2MGYiIOBe4EDhd0lOBvwXuBE4asIq3AZ+PiMfLDkn7SgpJJ0j6qKTbJN0v6R8kLZW0v6SLJD2QJ0lqUyznAQdJetEcL9F0hJ2MGYZ3kMYCXUIaNX5CDDBSWNI+pAF/3+qR5X3AnqTXnj8GXkN6FfoS8C+kKSYuB/5O0rOrsquB+4Cjhr0YMx4WzdglM/9ExE2STgdOJk0KdcGARV+Y0x/2OH59RBSVcpGkXwLeCLwxIv4BJlYvfCVwLHBVw6YnJF1OGqxnFiBWMmZgJD2F9OMP4PmSth+w6J45vbPH8XrCrmtyelHZEWkWwTuYOs1A4c7GOcwCw07GDMPHgZ1Iw/2XAR8dsNw2OX2kx/F7qu1H++zfhuk8RJqcyixA7GTMQEh6KfB24AMR8TXgT4HfGTDgWibk2qkj83YG5mVqBDMzdjJmRnJHu7NIa5l/Mu/+GGkyqLOU1jXvR3n92a8bC1kBXNtR3WaO2MmYQfgQ8DRSp7snACLiMVLflwOB/zVD+e+RXpUOG7VhuZPfM4CLR123GQ12MqYvkg4ldcQ7NSKuaB6LiKJsTm5pWm7mexj4CvBrHZj4ClKsZmSd/Mxo8cx4ZixIWklaGG7fyOshjajer5EmDX/jqOo0o8VOxowNSf8KXBcRJ46ovkNIHQOfExFrRlGnGT1+XTLj5F3AulGNwgZ2B95sB7OwsZIxxnSKlYwxplPsZIwxnWInY4zpFDsZY0yn2MkYYzrl/wN0DaaNzgs10gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## . . Space axis parameters\n",
    "xmin,xmax,ymin,ymax=0,1,0,1\n",
    "nx,ny=101,101\n",
    "\n",
    "xx = np.linspace(xmin,xmax,nx)\n",
    "yy = np.linspace(ymin,ymax,ny)\n",
    "dx = xx[1]-xx[0]\n",
    "dy = yy[1]-yy[0]\n",
    "\n",
    "## . . Initial Condition\n",
    "sigma=0.05\n",
    "mx = (xmax-xmin)/2\n",
    "my = (ymax-ymin)/2\n",
    "U = np.zeros((nx,ny))\n",
    "\n",
    "xxx, yyy = np.meshgrid(xx, yy, sparse=True)\n",
    "\n",
    "U = np.exp(-(xxx-mx)**2/((2*sigma)**2)-(yyy-my)**2/((2*sigma)**2))\n",
    "\n",
    "plt.figure()\n",
    "plt.imshow(U,cmap='jet',extent=[xmin,xmax,ymin,ymax])\n",
    "plt.xlabel('X (m)',fontsize=16);plt.ylabel('Y (m)',fontsize=16);\n",
    "plt.title('2D Initial Condition',fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Figure 9. 2D initial condition for the solution of the 2D diffusion-advection equation below.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
       "  <source src=\"./movies/Ex7_8.mp4\" type=\"video/mp4\">\n",
       "</video>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "## . . Velocity components parameters\n",
    "vx,vy = -0.1,-0.1 # velocity components (m/s)\n",
    "K = 0.01 ## . . diffusivity (m^2/s)\n",
    "\n",
    "# Boundary conditions on left, right, bottom and top\n",
    "LB=RB=BB=TB=0 \n",
    "\n",
    "## . . Time axis parameters\n",
    "nt,dt=501,0.0125\n",
    "\n",
    "## . . Solution space\n",
    "p2d = np.zeros((nx,ny,nt))\n",
    "\n",
    "## . . Loop over time\n",
    "for i in range(nt):\n",
    "    U = Solve_2D_Convection_Diffusion_homogeneous_Thomas(U,K,vx,vy,dx,dy,dt,LB,RB,BB,TB)\n",
    "    p2d[:,:,i]=U\n",
    "\n",
    "def convection_diffusion_2d_animate(i):\n",
    "    global k\n",
    "    ax1.clear()\n",
    "    ax1.imshow(p2d[:,:,k],cmap='jet',extent=[xmin,xmax,ymin,ymax])\n",
    "    ax1.grid(True)\n",
    "    ax1.set(xlabel='X (m)',ylabel='Y (m)')\n",
    "    ax1.set_title('2D Advection Solution',fontsize=20)\n",
    "    k += kskip\n",
    "\n",
    "## . . Animate solution\n",
    "k = 0\n",
    "kskip = 5\n",
    "\n",
    "## . . Set up movie\n",
    "fig8,ax1 = plt.subplots(1)\n",
    "fig8.subplots_adjust(0.1,0.1,0.9,0.9)\n",
    "fig8.set_dpi(100)\n",
    "\n",
    "## . . Call the animator.  \n",
    "anim8 = animation.FuncAnimation(fig8,convection_diffusion_2d_animate,frames=int((nt-2*kskip)/kskip),interval=100)\n",
    "anim8.save('./movies/Ex7_8.mp4')\n",
    "plt.close()\n",
    "\n",
    "HTML(\"\"\"\n",
    "<video width=\"500\" controls style=\"display:block; margin:0 auto;\">\n",
    "  <source src=\"./movies/Ex7_8.mp4\" type=\"video/mp4\">\n",
    "</video>\n",
    "\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "**Figure 10. Solution of the 2D advection-diffusion equation for the initial condition shown in Figure 9.**"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}