{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "###### Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 L.A. Barba, G.F. Forsyth, C. Cooper. Based on [CFDPython](https://github.com/barbagroup/CFDPython), (c)2013 L.A. Barba, also under CC-BY license." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Space & Time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Stability and the CFL condition" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Welcome back! This is the second Jupyter Notebook of the series *Space and Time — Introduction to Finite-difference solutions of PDEs*, the second module of [\"Practical Numerical Methods with Python\"](https://openedx.seas.gwu.edu/courses/course-v1:MAE+MAE6286+2017/about).\n", "\n", "In the first lesson of this series, we studied the numerical solution of the linear and non-linear convection equations, using the finite-difference method. Did you experiment there using different parameter choices? If you did, you probably ran into some unexpected behavior. Did your solution ever blow up (sometimes in a cool way!)? \n", "\n", "In this Jupyter Notebook, we will explore why changing the discretization parameters can affect your solution in such a drastic way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the solution parameters we initially suggested, the spatial grid had 41 points and the time-step size was 0.25. Now, we're going to experiment with the number of points in the grid. The code below corresponds to the linear convection case, but written into a function so that we can easily examine what happens as we adjust just one variable: **the grid size**." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy\n", "from matplotlib import pyplot\n", "%matplotlib inline" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Set the font family and size to use for Matplotlib figures.\n", "pyplot.rcParams['font.family'] = 'serif'\n", "pyplot.rcParams['font.size'] = 16" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def linear_convection(nx, L=2.0, c=1.0, dt=0.025, nt=20):\n", " \"\"\"\n", " Solves the 1D linear convection equation\n", " with constant speed c in the domain [0, L]\n", " and plots the solution (along with the initial conditions).\n", "\n", " Parameters\n", " ----------\n", " nx : integer\n", " Number of grid points to discretize the domain.\n", " L : float, optional\n", " Length of the domain; default: 2.0.\n", " c : float, optional\n", " Convection speed; default: 1.0.\n", " dt : float, optional\n", " Time-step size; default: 0.025.\n", " nt : integer, optional\n", " Number of time steps to compute; default: 20.\n", " \"\"\"\n", " # Discretize spatial grid.\n", " dx = L / (nx - 1)\n", " x = numpy.linspace(0.0, L, num=nx)\n", " # Set initial conditions.\n", " u0 = numpy.ones(nx)\n", " mask = numpy.where(numpy.logical_and(x >= 0.5, x <= 1.0))\n", " u0[mask] = 2.0\n", " # Integrate the solution in time.\n", " u = u0.copy()\n", " for n in range(1, nt):\n", " u[1:] = u[1:] - c * dt / dx * (u[1:] - u[:-1])\n", " # Plot the solution along with the initial conditions.\n", " pyplot.figure(figsize=(4.0, 4.0))\n", " pyplot.xlabel('x')\n", " pyplot.ylabel('u')\n", " pyplot.grid()\n", " pyplot.plot(x, u0, label='Initial',\n", " color='C0', linestyle='--', linewidth=2)\n", " pyplot.plot(x, u, label='nt = {}'.format(nt),\n", " color='C1', linestyle='-', linewidth=2)\n", " pyplot.legend()\n", " pyplot.xlim(0.0, L)\n", " pyplot.ylim(0.0, 2.5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's examine the results of the linear convection problem with an increasingly fine mesh. We'll try 41, 61 and 71 points ... then we'll shoot for 85. See what happens:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEbCAYAAACLNQJjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl4VNX5wPHvm40khABJCLIKKCCLIQqKO0EFFyoFpcXdQhVRq1K1FbVa19Zq3ZdaqhWr/qrWDTeUWgggIggSFtlkEwhC2LOQPef3x50JyTCTTJKZ3Lk37+d55pmZO/feeW9m5s05555zrhhjUEopp4myOwCllGoMTV5KKUfS5KWUciRNXkopR9LkpZRyJE1eSilH0uSllHKkGDveVEQygZuAEz0xxAJfAg8ZY3bXs+0W4ICfl+4wxnwZ4lCVUhFK7OikKiJrge+Bq40xRSLSBfgfVklwkDGmuI5ttxhjejRPpEqpSGVntfFOY0wRgDEmF3gc6A1caGNMSimHsKXaCGQYY8p8lu3w3Ldv7mCUUs5jS8nLT+IC6AMYYF5924vIYyKyRETWi8gsERkd8iCVUhHNrpJXLSISDUwEXjHGrK9n9TxgGXAvUAVMAmaIyM3GmOfreI9JnnWJj48f3L1795DEHomqqqqIinLniWQ3Hxu4//jWr1+/xxjTIRT7sqXB/oggRO4HLgKGGWMKG7H9p8CZQLoxpqS+9fv27WvWrVvX4DidIjs7m6ysLLvDCAs3Hxu4//hEZKkxZkgo9mV7iheRCcAvgQsak7g8FgFtgAEhC0wpFdFsTV4ichVwO3C2MSYviPUTRCTJz0uVnvvoUManlIpctiUvEbkSuBM41xiz07PsZ562Ke86HUWkZozjgSf87G4wUAqsDmPISqkIYkvyEpErgH8A04FzReRKTzK7COjsWed0rO4TL/hsfpmInFRjX+OBMcBjTah2KqUcxq6zjc8B8VgdU3094LkvBA4CP9V4baZnmxdFJBZoB+wHJhtjpoUvXKVUpLEleRljUoJYZzmQ4rNsF/CQ56aUasFsP9uolFKNoclLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOZFvyEpFMEfmHiCwVkeUislpEnhWRDkFsGysiD4nIWhFZJSJfi8gZzRG3Uioy2FnyegtIAc4yxgwCRgAjgQUiklDPts8B44EzjTEDgX8C/xWRzHAGrJSKHHZXG+80xhQBGGNygceB3sCFgTYQkb7AJOBRY8xuz7YvA5uAR8IesVIqIsTY+N4Zxpgyn2U7PPft69huLCDAHJ/ls4HJIpJkjCkMUYyOtLOoio27/f8JUlvH0S4xDoCCknLyCkoD7qdnamuiogSA3APFlJRX+l0vqVUMHZPjASirqGLb/kMB99m5bQIJcdEA7Cks5WBxud/1YqOi6J6aWP180+5CDHCgpCrgvlXLYlvy8pO4APoABphXx6YZQBWw1Wf5Zqzj6Q8sDkWMTnX/18WUzJ/r97V7LuzHdWf1AmD22jxufSsn4H5WPXAeSa2sr8htb+ewaPM+v+uNHtSZZy87AYDt+w9xzhP+3xvg39edwqnHpALwUvZGXv5qs9/1eqa1Zs4dWdXPL3ruK4rKrOSZ0ms3Z/Wpt2m0YQ7tg9KCI5fHtYbWaaF9LxUSdpa8ahGRaGAi8IoxZn0dq6YBh4wxvsWAfM99aoD9T8KqbtKhQweys7ObFnAES4s3VBj/LQI/bd1IdraV9zflVXBUogTcz1fz5xMfY70eVVoScN3Sg3nVf8/dh6rq3OfqlTmUbrNKXgfyygKu29qU1PqM0uINGCgqh4++WkbVjriA7xGshEO5dNi9kLQ935Bc8EPA9QqSerEn7RT2pJ1CUevuIIGPr6kKCwtd/d0MJTHG2B0DACJyP3ARMKyuap+IzAJONca08Vl+HTANuNAYM7Ou9+rbt69Zt25d04OOUNnZ2WRlZdkdRsi9MGcDj3+xjsnDjmHqBcc1bif7NsOyN2DtJ7B77eHlMfHQOv3I9Q/tgfIa1eCUXtDvIsi8Ejr0aVwMdXDrZ+clIkuNMUNCsa+IKHmJyATgl0BWEO1Ve4BEEYn2KX15k9necMToBN9s2svWvYeoKnJnu1ByQixAwHayeq35BD64Hso8X7H4ttDnfCsZHXO2VUX0VV4Cm7Jh7cewbibs2wQLnoFFf4eLnoFBlzYuFtVkticvEbkKuB042xiTF8QmK4DLgG7AlhrLewIVwJpQx+gU73+3nXeWbGfCgKZXqSJRW0/yyi9pYPKqqoLsP8O8x6znfUfByddCjzMhOrbubWPjoe/51q2yArYuhKXTYdW7ViL8aTmMeAiibf8ptTi2dpUQkSuBO4FzjTE7Pct+5mmf8q7TUURqxvkBVqN+ls/uhgOzjDF+Wl1bBm+JJDE2fG0ydkptHUdqvFQnsaCUHIS3LrMSl0TBuQ/ApW9aJa36Epev6BjoeSZc8jKMehKiYuGbF+H1MVC0p2H7Uk1mZw/7K4B/ANOBc0XkSk8yuwjo7FnndKzuEy94tzPGrMNq27pLRNI8600EjgHuac5jiDT5xRUAtHZp8jr92DSeyErkT2OPD26D3evhH+fA+s8hvh1c8S6cMaXpDe4icNKv4VefWO1kW+bDtCyrFKaajZ1l3eeAeKyOqb4e8NwXAgeBn3xevxn4I1Zv/HKgABhpjAl83r8F8FanErUGA7lL4bWfQ1kBpA+AS9+wGttDqfspcP1cePsqyF0Cr4yEy9+BXsNC+z7KLzv7eaUEsc5yrCFEvsvLgT94bsrD7dVGL2MMUlfpqbwY3r/eSlz9LoKxf/ffGB8KyZ1hwmfwyW8h50348Aa4caF1MkCFld3Dg1QI5XuTV4w7k1d5ZRU3/6+Ifvd9Tp1dfOY8Ant/gA7HwcUvhy9xecW0gouehS6DIT8XvmjRrRfNRpOXS1RVGcoqrS4SCS6tNsZGR1FaCSXlVRQHGKrEtsXw9fNW4/zPX7TOFjaH6BgY8zeIbgXLXocfvmye923BNHm5RFSUsPahC1j70PlER7mz5AWHq8R++3qVF1vVNgycfit0Hdy8wXXoC8Pvth5/fIt1plOFjSYvl4mPjbY7hLBK9PRu8J5ZrWX2w7B3g1VdzLqreQPzOu1m6DLEU328254YWghNXspRvO15R5S8ti6ChS+ARMOYF612KDtERdeoPr4BP/zXnjhaAE1eLrH0x32c99Q87v/oe7tDCStvtTG/ZvIqO1S7utilmauLvjr0gbM9jfYf3QLFB+yNx6U0ebnE7oJS1u0qYMeBYrtDCavWnpMRtYYIzXkE9m2EDv0ga6o9gfk69TfQ9SQo2KHVxzDR5OUS3jag5IYMnXGgUzrHcN/P+pPR1dOP6sBW+OZv9lcXfUVFW2c7o1tZ/b92ubtEbAdNXi7hbQNKjnd38hrUIYaJZ/Tk2HTPJCKL/g6mEgZeAl1OtDc4Xx36wOBrrMcLX7Q3FhfS5OUS3mpUgwYtO11JPix9zXp86o32xhLI0MmAwMp3oGCX3dG4iiYvl6guebm1h6rH3uIqPli2na9+2GN1Bi0rgKPPgM4n2B2af6nHwHGjoLIMvn3Z7mhcRZOXS3jPvrm95LX5YBW/fXs5r3+9Eb55yVp46k32BlUfb3zfvmx1pFUhocnLJc7o3YHLh3and3qb+ld2MO90P8ftz4aDW62ZIvqcb29Q9el+qlUyLN4Hy9+yOxrX0OTlEuMGd+VPY4/n+K7uns3A28P+vIL3rQen3AhREf41FrG6ToB1ZrTKndN0N7cI/9SVqi0xRjhR1tO/cq01wWDm5XaHFJz+P4fkLrBnHWz8n93RuIImL5dYtGkvK7YfoKLS3f/VE2OFiTGei0MNmRj+6W5CJToWhl5vPV74vL2xuIQmL5e45tXFjH5+AaUV7k5e7St2cUHUYspMNJUnXWd3OA1z4jUQ29q6GtHOVXZH43iavFygtKKSkvIqYqKExDh3zyrRLfdTosXwSdWpFMQ67ErWCe3gxKusx99op9Wm0uTlAjWHBtU5PbLTlRyk00/WLA1jJj9Cu0QHXuLN22l1xTtQsNPuaBxNk5cLtJje9cveIKayGHqcSVSXTLujaZyUnlan1apy+PYVu6NxNE1eLnB4XKOLe9cbY82PBYcbvp1q6GTrPuf/tNtEE2jycoH86qFBLi557VwJeaspj2nDM1t78bPn5jNnXTAXWI9AR58ObbtD/nb48Su7o3EsTV4ucLAlJK8VbwOQl34GuQUVrMrNZ+fBEpuDaqSoKMj4pfV4+dv2xuJgmrxcYET/jsy+fRhTzz/O7lDCo7ICVv4HgF0dh1e37eX7uwiHUwy61LpfPcOaCVY1mCYvF0iMi6FXhyS6pSTaHUp4bM6Gwl2Qcgz5yX2q5yzzewUhp0jrDZ1PtGbFWPeZ3dE4kiYvFfm8VatBl4IIbT0DHGtNBe1Egy6z7ldo1bExNHm5wPQFm7n1rWUs3rzP7lBCr7QA1nxsPfa0E3lLXn4vf+YkAy+BqBjY8D8odOjJBxtp8nKBxVv2MSNnB7vyHdqAXZc1H0NFsTWtTPsewOEJFx1dbQRonQrHjrCmsV75rt3ROI4mLxdw9cU3vPNfZYyvXtQ9pTW/HNKVrL4dbAoqhAZ5jmuFzvPVULb2ahSRTsCrwHnGGBePawkv1/awP5gLm+dBdBwMGFO9+Nj0JB4bN8jGwEKozwXQqi38tBzy1todjaPYVvISkbHAQuCYRmybLSKrRSTH53Z16CONfK7tYb/yP4CxZkpNaG93NOERGw8Dfm491tJXg9hZbZwKjAAWNHL7C40xmT63f4UwPsdw5fz1xhyuMnrPylW/ZFjzUz6LNu21IbAwqD7r+A4YHS4ULDuT1+nGmB9sfH9XMMaQX2K1ebVx0zUbd66A3WsgIQWOPfeIl3/+/ALGT/uGkvJKG4ILsW6nQLvukJ9LuwM6z1ewbEtexhiHn+eODGWVVZzVO41Te6USF+Oi8y/evl0DL4GY2lPfiEj1GUdH97L3ioqqPiHRcdccm4NxDid/228TkcUislZE5onIBLsDskOrmGhenXAy/550it2hhE6N4UDVw2h8eM+sOr67hFeGdZwddn+tw4WC5NQW3gPABuAuoAQYC7whIgOMMXf420BEJgGTADp06EB2dnYzhdr8CgsLHX18KXu/I6Moj0MJnVn8QwFsyK5+rfrYyqzrH2Z/vZjc9u6YPfbENr1JLviB1R/8lbyOZ9kdTsRzZPIyxozxWfSuiAwHfisizxpjtvrZZhowDaBv374mKysr/IE2g5LySgpLK0iOj62uNmZnZ+Po43vPmrcr8ZSJZA0bXusl77G9umkxmw7u5ph+A8k6rqMdUYZe4vXw2R30L19B/6z77I4m4jm52uhrEdbxnGR3IM1p3vrdDHn4S258c6ndoYRGaQGs+cR67J02xo/DM0u4qOl0wMVUSbR1abSCXXZHE/Ecl7xEJE5E/F1Z1XvayR11iCB5zzQmu+VM4+qPPMOBToP2RwdczTVDhGpqncq+lMFWd4lVOlyoPhGfvEQkVURqnm46DXjHz6qDPffLwh9V5HDdRITejpqDxte52uRhx/DFlLMYe2KXZgiq+ezqmGU9WK4dVusT0clLRHoCucAMn5fOEZFRNdbLAq4HXm9pfcdcNQX0we2weT5Et4L+vs2atXVtn0jfo9q4p8TpsTf1JIhva/Vz27Xa7nAimp3Dgx4XkRxgtOe5d4hPzVJWMbAP2FFj2XfA74G7RWS5iGwAXgQeBiY2T/SRw1VDg1a8Axjoe4F1jcMWqCo6DgaMtZ7ocKE62faNN8b8Loh1dgKdfZblA096bi2eawZlG3N4Ur4AfbtqWr0jn1e+2syx6UnckNXg4bGRLeNSWDodVvwHzvkjRLWoZtygRXS1UdXPNdPh/LQcdq+FxFS/w4F87Ssq473vtjP/h93NEFwz634KtDsaCnbAlvl2RxOxNHk53C3nHMuLV5xIZjeHV7O8pa6B4yC6/kTsyrONXiKHS596daGANHk5XEbXdlx4fCc6JsfbHUrj1RoOVPdZRq/qfl5On8c+EO/ki6tnQFmRvbFEKE1eyn4bZ0PRbkj1XFEnCK6Zxz6Q1GOg60lQXgRrP7U7moikycvhnpy1jr9lb6S80sHzQNXs2yXBTajbxnN2Nb+knKoqE67I7OUtfWmfL780eTlYeWUVz87ewONfrCU6yB99xCnJP1yyyAiuyggQEx1FUqsYjIHCMpeWvgZeAlGxsGkOFOy0O5qIo8nLwQpKDp9pjIpyaPJaPQMqSuDoM6wJ+Rrg5J4pnHFsGuUVDi511iUxBfqcZw0X8rYJqmou6NnYch3uoOrgbhLLrBkkgm2or+mfv2oBY/AzxsPaT6y/06m/Cbpa3RJoycvBHD93/falsO0b6+o53l7lqrY+50PSUVYfuI2z7Y4momjycrDDg7IdWoD+5gXrfvA10KpNgzevqjIcPFROcZkL5rEPJCYOhk6yHi98wd5YIowmLwfz9nFyZLXxwDb4/kOQaBh6faN28fv3VjDowVl8vHxH/Ss72eAJEJtozfOlg7WrhTR5iciLodyfqpsgdGob78wOqov/bl3mfsBYaNu1Ubuo7uvl1o6qXokpkHm59fgb/Yl5Nai+EcRFXS9sQiyqgUZldGJURie7w2i40gJY+pr1+NSbGr0bVw8R8jX0Bvj2FWvmjXPug6R0uyOyXUMbS6bX8ZpLewqqkFv2BpTmW7OldgmuR70/h3vZt4DklXasNVXQus+sJDb8Lrsjsl1Dq41rgJ4+t4HAOOADoP7pAFTLVlUJ3/zNetyEUhfUHN/o0k6qvrx/r29fhvISe2OJAA1NXrcYY370ua02xrwPXAVMCUOMKoDb3slhyMNf8uVqB12sYe2ncOBHaN/TKkk0geuu3Vifo0+HToPg0B5Y6W8m9JalQcnLGPO/Ol4rBo5rckQqaHsKy9hTWEp0tIM6LnpP959yY5Mn2Tt8BaEWkrxErI6qYP0dTctuqWlog72/K2EK0B4Yg3UBWNVMHNfDfvsSq1NqfNvDZ8+a4Nj0JJ65NJNObRNCEJxD9B8D//2jp9Pq/4KauNGtGtpgn43/hnkBtgNXNjUgFbwCp/Ww95a6Bk+AVklN3l1K6zh+numuqwfVy9tp9cv7rb+nJq+gbQSu9VlWCeQBG40xLu7qHHkc1cN+/xZrEHZUDJw8ye5onG3wr2DuY9ZwoZ0r4ajj7Y7IFg1tsH/BGDPX5/aVMWa9Jq7mZYxxTg97Y+CT26xOqQPHQdvQlZZeXbCZJ2eto8ytM0v4k9AeTrzGevzJb60zuC1QQxvsnw5XIKphissrKa80tIqJIj42wq8u892/rPaZhPYw4sGQ7vq52Rt4dvaGlnPG0StrKrTpDNu/hYXP2x2NLXRso0NFifCHUf2Ycm4fu0Op24Ft8MU91uMLHoc2HUO6e9fPZR9IQjsY/az1ePYjsHudvfHYQJOXQ8XHRnPtmb0i+5qFxsBHN0NZARz3Mzh+XMjfwnux3RbTXaKm3iMg80qoLIUPb2xx1UdNXip8vvuXNYVxQnsY9WRYJtJrcR1VfZ33iFV9zF3S4qqPmrwcatu+Q8zIyWXF9gN2h+JfzerihX8NeXXRK7mlDRHy1YKrj5q8HGrJj/u49a0cXvlqs92hHMm3ujjwkrC9lfdMa4steYFVfTyh5VUfNXk51MFDEdxNYun0sFcXvdKS4khLisNBA6TCY2SN6uPXz9odTbNwQO9G5Y+3mhRRveuNsWaMmPUH63kYq4tet4/sy+0j+4b1PRzBW318c5zV+76yHM68A6LcWz5x75G5XMT1ri8vhg9vgC/usjqjnnlHWKuLyo/eI2DEQ4DAnEfgP1dbEz+6lK3JS0Q6icjnItKyh8c3Qn4kDco+mAuvXgDL/23NtT7uVTjnXr1Mlx1OvwUuf8e6ItOaj+HlEbB3o91RhYVtyUtExgILgUZ1VBKRKSKyWkRWiMh3IjImtBFGNm+nTNurjT8uhGnDYMcy66Kxv54FAy9utrdfsmUfp/75f1z3ryXN9p4Rr89IuG42pPWF3WvgH8NhQ8DZrBzLzpLXVGAEsKChG4rIVOAPwEXGmAzgTuA/ItK02e0cJL/48NWym11FGaybCe/+Gl67CIp2Q8+zYNLcZh8kHB0l/HSwhF35OhtTLWnHwrVfQt8LoeSg1Rb29lXW4PjyYrujCwk7G0xON8ZUSAOrFiLSDrgXeMIYsxHAGPNfEZkF/BWYGfJII9Ab1w6lsKSC+Lhm+v9TVQmb58Gq92DNR9YPwuuUG622lujm/zq1uAkJGyI+Gca/CXMfhXmPW5/bmo8grg0cN8pqkzxmOERHQNNDI9iWvIwxje1VeD6QCMzxWT4b+KuIHGeMWVvXDqoqStn8/SK/rx3VNoEEz0DnvYWlAcfMxURH0a19YvXzzXuKCHQNkpTWrap/ZAUlFewpDFxKODqlNVFRVkLfcaCY0gr/fXZat4olvU2rI18whtaFW2DnKt8XwFR5Zt801r0xUFUOFSVWaaqyFCpKrUbe/B1wcDvkb/fc77DW8+o40PryD7wY2vcIeDzh5i157j9Uztqd+SS1iqGr53Mpq6hi057CgNt2bZ9IUivrJ5BXUMK+ojK/68VERXFs+uH5x37YVUBlgFlM05JakZZkfS75JeXsOBC4lHNshyRioq1/Ptv2HaKorIJtBVWs3Zlfa70mH1O/3xDT7RLabPyY5A0zSNi9Ala8Zd2i4yC5M4cSOlGe1JnypM5UJHWmMjYJE90KEx1HclISbZOSIDqOwvJKdheWAeJp04yq1bbZrX0CMZ4znDvzSykuD1/n4Qg5VdUgGZ57396Zm2u8Xmfyalu8jZ7/GVnvG6V6bsHoGeR6bTy3YHQOcj1fJwGEowko5RhPwroE0iNjxu+2NYYHnf/0fM4fcBQvXTUYgF35JZz/9PyA2/5r4smc1acDANMXbOHFbP8N213aJbBg6tnVz3/x94UcOOT/n9rvzuvLTcOPBWDBD3u44c3vAr7/sntH0L51HAB3f7CS+T/ssV5YUDvm0B1TP6AfR8tOrmi9hEntl1ltYvu3kLh/S8B91pTkuQXjqCDXaywnJq80z73vOWDvvyu/+UZEJgGTAAZ2imcj3fzuPCVeiPXMCV9QVkWA7ygxUZCacLjKlldUFfDab23ihMRYa5/FFYb80sAnVzsmSvV/sn3FVZQHmKYqIQaSW/mvMlZVVRHlt39PFEbAmvjWeo+qqBiqouKoioqlKioWI7FURsdT2iqV0lZplMR3oLRVGqWt0qiM8Uy3vHqndbNBYWEh2dnZtZad3yOWVXus//BVhXuqX99fUkXXpMDNEuu+X0HVDquUfWBnecB120WX1XrPDnGVJEX5X3f39s1kZ28HYOOeyjrf/+uvF9Da872IKi6la5JQWVVFtM9nF+pjqqQTH7X6OX0GjCeqspRWpbt5b9kO2lfuJp29pJt9JFBCHOXEUsFRrSpIiatATAUl5Yb9JZUIxvMtqv0FTU+Mwvun2V9SRdkRFYfvA8beYMYYW29Y14I0DVh/Glb9LNVn+QjP8hvq20efPn2Mm82ZM8fuEMLGzcdmjPuPD1hiQpQ7nNhJ1VO2PqL25X2+txljUUrZxInJa4XnvofP8p4+ryulXCzik5eIpIpIXI1FnwOHgCyfVYcDq009ZxqVUu4Q0clLRHoCucAM7zJjzAHgIeAmEenlWe9c4DzgDjviVEo1P9vONorI41iN7N09z3M8L51sjPF2uCkG9gE7am5rjHlUREqAT0SkAuvya78wxrSIDqpKKXs7qf4uiHV2EqC7k7GuZKRXM1KqhYroaqNSSgWiyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSE68bmOzyM/PJy8vj/Jy511Gvm3btqxZs8buMMKiMccWGxtLeno6ycnJYYpK2UGTlx/5+fns2rWLLl26kJCQgEjgi3xGooKCAtq0Cfa63M7S0GMzxlBcXExubi6AJjAX0WqjH3l5eXTp0oXExETHJS5Vm4iQmJhIly5dyMvLszscFUKavPwoLy8nISHB7jBUCCUkJDiyCUAFpskrAC1xuYt+nu6jyUsp5UiavJRSjqTJS0WMl156if79+yMiTJ8+3e5wVITT5NUC5OXlkZmZSUpKCiJCZmYmL7/8coP2ceeddzJkyJBay3Jycrj//vs5cOBAreWlpaX07NmTZ555pkHvMXnyZD777LMGbaNaLk1eLUB6ejo5OTmMHj0asJLOtdde2+B9dO/evdaynJwcHnjggSOSV0xMDN27dyc1NbVpgStVB+2kqoJy++23c/vttwe1bnR0NHPnzg1zRKql05JXC+atSvbo0YOZM2dy9tln07VrV0aMGMH27dur17vpppvo3r07IsKWLVsAuOeee7jvvvsAuPDCC8nMzGTcuHHs3buXzMxMkpKSyMrKqt5HSUkJU6dOZfDgwZx44olkZGQwefLkI0ptSgVLk1cL5q1K7t+/n4ULFzJ79mzWrl3L1q1b+f3vf1+93gsvvMCDDz5Ya9tHHnmketlnn31GTk4O7777LqmpqeTk5BzRPnbgwAFeffVVPvzwQ7777jsWL17Mvn37uPrqq8N/oMqVtNrYAD2mfhrwtT+NPZ7Lh1ptQv+3aCt3f7Ay4LpbHh1V/fhnz81nVW6+3/UuO7kbf744A4CV2w9yfNe2jQm7XgUFBUyZMgWApKQkRowYwfvvvx/S90hLS+Prr7+mW7duAMTHxzNx4kQuuOACdu3aRceOHUP6fsr9bEteIpIOPAV4/0WvBKYYY7YH3qp62y2Av/rGHcaYL0MWZAuRlpY88iVJAAAPfklEQVRGSkpK9fOUlBR27doV0veIiYlh/fr13HjjjeTm5hITE0NhYSEAmzZt0uSlGsyW5CUiccB/gfXAAMAA/wTmiMgJxpjC+vZhjMkMb5RHqlliqsvlQ7tXl8Lq88nNZwa1XrhKXQCJiYm1nkdFRVFVVRXS95g5cyajRo3iySef5NZbb0VEyM7OZvjw4ZSWlob0vVTLYFeb1zVABnCnMabCGFMJ3An0Am6wKSYVRq+//jpJSUlMmTJFxxmqkLAreV0CbDXGbPIuMMbsBFZ7XlMOEBsbC1hzZgF88cUX7Nu3z++6paWlREXV/rrt3LkzvAEqV7MreWUAm/0s3wwcH8wOROQxEVkiIutFZJaIjA5phKpePXv2BGD79u0UFBQwduxYCgoK/K47atQoDh48WN2zv6CggKeffrrZYlXuY1eDfRqw1M/yfCBRRBKMMcV1bJ8HLAPuBaqAScAMEbnZGPO8vw1EZJJnPTp06EB2dnbAnbdt2zbgj9AJKisra8W/e/duxowZw7Zt2wDIyMhg0qRJvPPOO6xcuZKioiIyMjL46KOPePTRR/nggw+q13vyySd5++23+fzzzwE4//zz+c1vfsM111zD8ccfz9VXX82VV15JQkICU6ZMwRhDRkYGmzZtqt7Hhx9+yLhx49i8eTMPPPAATz31FEcddRRZWVksWrSIiRMnMnnyZOLi4njppZcAuPfee5k/f/4RCc732BqipKSkzs89EhQWFkZ8jJFCvEX+Zn1TkTLgC2PMRT7L3wQuBxLrSV7+9vkpcCaQbowpqWvdvn37mnXr1gV8fc2aNfTr168hbx9RdBpo/5zwuWZnZ9fq3Os2IrLUGDOk/jXrZ1e1cQ/g7xvYBjjU0MTlsciz/YCmBKaUcga7ktcKoIef5T2x+nsFJCIJIpLk56VKz31000JTSjmBXcnrfeBoEenhXSAiHYF+wHs1VxSRjiJSM87xwBN+9jkYKMU6Y6mUcjm7ktd0rBLWX0QkxpOcHsU62/g370oicjqwA3jBZ/vLROSkGuuNB8YAjwXTwVUp5Xy2nG00xpSJyAis4UGrsXrYrwLO9kk+hcBB4Kcay2YCjwMvikgs0A7YD0w2xkxrjviVUvazbWyjMWYX1pnFutZZDqT4LNsFPOS5KaVaKJ0SRynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJS9Xy9NNP8+GHHzbre5aXl/Pvf/+bkSNHMnjwYAYMGMCQIUN47bXX8DdxwNKlSxk2bBgDBw6kb9++3HHHHZSU1DkWX7mQJi9Vix3Ja+nSpVxxxRVMmDCBpUuX8v333zN16lR+9atfVV9ezWvDhg0MHz6ciy++mFWrVrFo0SK++OILJkyY0KwxK/tp8lIR4bTTTuOyyy6rfj5u3DjOOOMMnnnmmVqlrz//+c+kpKRwyy23ANCuXTvuu+8+3nrrLb799ttmj1vZR5OXy3300UdkZmYiItx7773ceeedDB48mK5du3LPPfdUr7dt2zYyMzPZsWNH9TaZmZl8+WX4L8Y0dOhQ5syZc8Tyzp07U1RURHl5OQAVFRV89tlnDBs2rNY8+GeffTYA77333hH7UO6l1210udGjRzN69GhEhNdff5333nuPv/zlL8yaNYvzzjuPYcOGMXLkSLp160ZOTg49evQgKyuL6dOn17vvJUuWcO2119a73pAhQ6qnf/ZHRKrnw69p/fr1nHrqqcTFxQHWJdKKioqqp5/2Sk1NpU2bNqxYsaLeWJR7aPIK1v3hu/RYg9x/sNGbZmZmMnjwYABGjhxJUlIS2dnZjBw5slH7GzJkCDk5OY2Opy6LFy9mxYoVtUpke/bsAfA7k2pycjJ79+4NSywqMmm1sQXp06dPreft27cP+cVlQ6GwsJBf//rXPPzww5x11llBbWPHdObKXlryClYTSjyRwt/FZSsrKwOsbY+ysjIuueQSRo4cyV133VXrtbS0NAC/F+AoKCggNTW1WWJUkUGTl2q0ULV5eZWVlXHxxRfTv39/nnjiyMlye/XqRevWrdmyZUut5Xv37qWgoICMjIygY1fOp8lL1RIbG1tdBfvxxx/Jzc3ltNNO87tuKNu8vCWu3r1789RTT1Uvv/7667n//vvp1KkTMTExXHDBBcydOxdjTPUZR2+72MUXXxySWJQzaPJStfTs2ZPt27cDMG3aNH766aeAyStUysrKGDduHJs2bWL8+PG88cYb1a/NmzeP0tLS6ud33303w4YN4/nnn+fmm2/m4MGDPPjgg1x66aWcfPLJYY1TRRhjTIu79enTx9Rl9erVdb4e6fLz86sfz5s3zwwaNMgApmPHjubGG280+/fvN4MGDTKxsbGmffv25pxzzqlef8GCBea4444zAwcONEOHDjXr1q0Le7wzZswwWFOB+71t3ry51rF9++235qyzzjL9+/c3vXv3NrfddpspLi6u932c8LnOmTPH7hDCClhiQvQ71pKXy5155pl+q3aBqnunnXYaa9asCXdYtYwePbpBZwuHDBnC3LlzwxiRcgLtKqGUciRNXkopR9LkpZRyJE1eSilH0uSllHIkTV5KKUfS5BVAQ07dq8inn6f7aPLyIzY2luLiYrvDUCFUXFzsd84w5VyavPxIT08nNzeXQ4cO6X9shzPGcOjQIXJzc0lPT7c7HBVC2sPej+TkZAB27NhRPQWxk5SUlBAfH293GGHRmGOLjY2lY8eO1Z+rcgdNXgEkJyc79suenZ3NCSecYHcYYeHmY1MNo9VGpZQj2Za8RCRdRN4UkXWe27si0jXIbWNF5CERWSsiq0TkaxE5I9wxK6Uihy3JS0TigP8CccAAoD9QBMwRkaQgdvEcMB440xgzEPgn8F8RyQxTyEqpCGNXyesaIAO40xhTYYypBO4EegE31LWhiPQFJgGPGmN2AxhjXgY2AY+ENWqlVMSwK3ldAmw1xmzyLjDG7ARWe16ry1hAAN+rlM4GRgZZclNKOZxdySsD2Oxn+Wbg+CC2rQK2+tk2BqsKqpRyObu6SqQBS/0szwcSRSTBGBOoi3sacMhT1fTdFsDv9a9EZBJWdROgVERWNTBmJ0kD9tgdRJi4+djA/cfXN1Q7irR+XhKubY0x04BpACKyxBgzpAnvFdHcfHxuPjZoGccXqn3ZVW3cAxx5zXZr2aE6Sl3ebRNFJNrPtgB6zXelWgC7ktcKoIef5T2BlUFsGwV087NtBdC8V49QStnCruT1PnC0iPTwLhCRjkA/4L2aK4pIRxGpGecHWJfEyvLZ53BgljHmyGvBH2law0N2FDcfn5uPDfT4giZ2zJrg6aS6BKuUdAXW2cNXgDOAE4wxhZ71TgfmAdOMMTfU2P4lrGR1ujFmj4hMBF4ATjXGhOYSzkqpiGZLycsYUwaMACqx+natAZKBs72Jy6MQOAj85LOLm4H/AAs8Zw2vA0Zq4lKq5bCl5KWUPyLSCXgVOM8Y05QzzxHHzcdmF9fMKuH2gd5NPL4tIpLj53ZuuOMOloiMBRYCxzRy+ykislpEVojIdyIyJrQRNl5Tjk1Esj3H5fvZXR36SBtORDJF5B8islRElntifVZEOgSxbdN+d8YYx9+wBngvx6pKxgDRwGvAD0BSENu/BKwHOnieXwsUA5l2H1uIjm+L3ccQRIyLgN7AdOtr2aBtp2J1oTnG83wEUA5cYPdxheDYsoEedh9DHfGtxTrJ1trzvItn2XogoZ5tm/S7s/3gQ/QHvA7rDGSvGsuOwmpT+1092/bFOmEw0Wf598Cndh9bU4/Ps+4Wu48hiBhjPPcN+oED7bBmJHnQZ/mnwPd2H1dTjs2zjROS17E+y37t+b5eUsd2Tf7duaXa6PaB3k05PkcwxlQ0ctPzgUT8f379ReS4JgUWAk04NifIMMZs8Fm2w3Pfvo7tmvy7c0vycvtA76YcHwAi8piILBGR9SIyS0RGhzRC+2R47n3/Ppt9Xney20RksadtaJ6ITLA7IC9j9Rzw1Qer5DWvjk2b/LtzS/JKA/x1Tq0e6F3Ptg0e6N3MmnJ8AHnAMuB0rMkfZwAzROQ3IY3SHmmee9+/TyR9fk1xANiA1a9xAPAs8DcR+autUQXgGbY3EXjFGLO+jlWb/LuLtIHZoRa2gd4RIqgYjTEn+yx6QUQuBP4kIi8bY0pCH5rtnPD51csY43vW9F0RGQ78VkSeNcb4llzsdi/WML3fNnL7oD83t5S83D7QuynHF8giz/YDmhJYBPBOH+P794mkzy/UFmH9dk+yO5CaPNXZX2Kd5S2sZ/Um/+7ckrzcPtC70ccnIgkBGj+9xXXfL4/TrPDc9/BZ3tPndccRkTgRaevnpYj77ETkKuB2rFEyeUFs0uTfnVuSl90DvcOtKcc3HnjCzz4HA6VYZywdQ0RSPWNjvT4HDuH/81ttjFnbXLE1lZ9jOw14x8+qgz33y8IfVf1E5Eqsa1Cc6zkLjoj8zDMBqHed0P/u7O4nEqK+JnFYmfxtrHa8KKyhGLU6cWI1WFcCf/PZ/iVgHZDmeT6RyOuk2qjjA36F1Qh6Uo1l47HO9DzYXMfQgGOdToC+UFj/lUuAmT7LpwK78fSDA84lgjqpNvbYPD/sCmCUz7Ii4F92H48nnis8v5U7gCtr3P4O3G8CfC89y5v0u3NFg70xpkxERgBPYZUkDLCKhg30/iPWQO9yrDNXETPQu4nHNxN4HHhRRGKxOnXuByYba3bZiCAij2P1jO/uee79259sDp+OLwb2cbgfEQDGmEdFpAT4REQqsH4ovzDGzGyW4OvRhGP7Dvg9cLeI/AloDZQBD2N9ppHgOSAe//E84LkPy+9OB2YrpRzJLW1eSqkWRpOXUsqRNHkppRxJk5dSypE0eSmlHEmTl1LKkTR5KaUcSZOXUsqRNHkppRxJk5dSypE0eSmlHEmTl4ponmtV5otIlYh86Vn2gojsF5HNInKt3TEqe+jAbBXxROQXWPNaXWeMeVlEjgY+Bk4z9c/YqVxKk5dyBBF5H2uerkzgn8CfjDGz7I1K2UmTl3IEETkKay6zSuATY0zEXP5L2UPbvJQjGGt64QewLpnle6FS1QJpyUs5gmf+82wgAWtG0gHGmD11bqRcTUteyiluxbrk1xisaYefsTccZTcteamIJyLHAO9inV0sFpHrsS7ecJEx5hN7o1N20ZKXimgi8gjwFXAU1tVlAG703L8pIu/aEpiynZa8lFKOpCUvpZQjafJSSjmSJi+llCNp8lJKOZImL6WUI2nyUko5kiYvpZQjafJSSjmSJi+llCNp8lJKOdL/A/k2gNPa2E41AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "linear_convection(41) # solve using 41 spatial grid points" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEbCAYAAACLNQJjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNXd+PHPNxtZCIEkJCBrRBYRQyyoFUUWgbryuNBSq7WVKqVWW1p9itXHVq22ttatVuWh2tKfbbWtax830ELAurAaUNlEQCAsIQRC9m3O7487MyTDTDKTTHLn3nzfr9e8ZubOPXe+NzPzzbnnnnuOGGNQSimnibM7AKWUag9NXkopR9LkpZRyJE1eSilH0uSllHIkTV5KKUfS5KWUcqQEO95URAqA7wNf8saQCLwD/MIYc6iNsruAo0Feus0Y806UQ1VKxSixo5OqiGwBPgWuM8ZUicgA4N9YNcGxxpiaVsruMsYM7ZpIlVKxys7DxgXGmCoAY0wx8CAwHLjYxpiUUg5hy2EjkG+MqQ9Yts9736erg1FKOY8tNa8giQtgBGCAlW2VF5HfiMhaEdkmIktFZGbUg1RKxTS7al4tiEg8MAd4xhizrY3VS4CPgLsADzAXeFVEbjHG/L6V95jrXZfk5ORxgwcPjkrsscjj8RAX584TyW7eN3D//m3btq3UGNM3GtuypcH+hCBE7gYuAyYZYyrbUf51YCKQY4ypbWv9kSNHmq1bt0Ycp1MUFhYyefJku8PoFG7eN3D//onIOmPM+Ghsy/YULyLXA18DLmpP4vJaBaQDp0UtMKVUTLM1eYnIN4FbganGmJIw1k8RkZ5BXmry3sdHMz6lVOyyLXmJyLXAAmCaMeaAd9ml3rYp3zq5ItI8xtnAQ0E2Nw6oAzZ1YshKqRhiS/ISkWuAPwCLgWkicq03mV0GnORd51ys7hNPBBS/WkTObLat2cDlwG86cNiplHIYu842Pg4kY3VMDXSP974SKAf2N3vtTW+ZJ0UkEegNHAHmGWMWdV64SqlYY0vyMsZkhrHOBiAzYNlB4Bfem1KqG7P9bKNSSrWHJi+llCNp8lJKOZImL6WUI2nyUko5kiYvpZQjafJSSjmSJi+llCNp8lJKOZImL6WUI2nyUko5kiYvpZQjafJSSjmSJi+llCNp8lJKOZImL6WUI2nyUko5kiYvpZQjafJSSjmSJi+llCNp8lJKOZImL6WUI2nyUko5kiYvpZQjafJSSjmSJi+llCNp8lJKOZImL6WUI2nyUko5kiYvpZQjafJSSjmSbclLRApE5A8isk5ENojIJhH5nYj0DaNsooj8QkS2iMgnIvK+iJzXFXErpWKDnTWv54FM4HxjzFhgOjADeE9EUtoo+zgwG5hojBkD/BF4W0QKOjNgpVTssPuwcYExpgrAGFMMPAgMBy4OVUBERgJzgQeMMYe8ZZ8GdgD3d3rESqmYkGDje+cbY+oDlu3z3vdppdwVgADLA5YvA+aJSE9jTGWUYnScQxV17K/08Pmhln+CtKQE+mUkA9DQ5GF3WXXIbfTrlUxaD+urUVZVz5HqwI/JkhAnDMlK8z/fVVpFkzFB1+2TmkRmWhIAlXWNHDxWG/L9h2SmkhBv/V/dX15DdX2T/7WjdZ6Q5ULyeKCmDGqOhFlAQHy3OJB4iE+ExBRITIN4O382yse2TyFI4gIYARhgZStF8wEPsDtg+U6s/RkNrI5GjE6z7osyvrrwAzwG+M+KFq9NOzWXp781HoDSyjoueGhFkC1Y/vTtM5kyKsd6/N5OHl+2Peh6/Xol8+EdF/ifz1r4PqWVwRPdj6eP4AcXDAfg/e2lzH12Xcj3X3PnNPqm9wDgpy99TOHWQy1ez8w7xPkj2mga3bYU3nsUjnwBlQfA09j6+pGIS4T0/pB9CmSPgOzhMPIS6NU/eu+h2hQz/0JEJB6YAzxjjNnWyqrZQLUxpilg+THvfVaI7c/FOtykb9++FBYWdizgGLRsd4OVuIB+qdLitabKw/59PlrnOeH15rZu+hg5EA9A6b6GkOv2iq9v8XfsndBIQoh1DxXvorCwGIDth5taff8PP3if9CTrdU9VnX/digZDVQO8/l4Rnn2JQcsm1ZVxyvanyTn0XovlDQk9aUhMx6q0t02Mp9m9hzhPI3GeOuKb6hBPA5Tvtm6fL7PifGMBB3Mns2fQFVSnDQzrPYKprKx05XezM4gJUc3vaiJyN3AZMKm1wz4RWQqcY4xJD1h+I7AIuNgY82Zr7zVy5EizdevWjgcdY8qrG9hXXsMnH63lqxdPtTucqHt46VZ+t2w786cNZ/60ESeusP5ZWHIn1JVDYipMuQNOvQx65lqHfNFgDDTWQnkxlG6zbntWw9Y3sA4aBE69FC59DNKC/h9tVWFhIZMnT45OrDFIRNYZY8ZHY1sxUfMSkeuBrwGTw2ivKgVSRSQ+oPblS2aHOyNGJ8hITSQjNZGDW+0+D9M5UpKsr2tNfWClG9j5LvzrZuvxiAvh4geh9+DoByFiJcLsU6yb79zS4c/h/ceh6G+w+f+g4gB86/+ilzTVCWz/lovIN4FbganGmJIwimzEintQwPI8oBHYHN0IVaz4xtmDeWhSir/tzM8YeOdu6/F5P4arn++cxNWarGFw2aNwyzrIGAR718CLN4AnSKJVUWFr8hKRa4EFwDRjzAHvsku97VO+dXJFpHmcL2PVzycHbG4KsNQYU9G5Uceuf6zZw23/3MCWMnf+YDJSEslKifOfCfXb8joUr4W0vjDxVqt2ZJfeg+Caf0KPDNjyGiy5w0quKurs7GF/DfAHYDEwTUSu9Sazy4CTvOuci9V94glfOWPMVqy2rZ+KSLZ3vTnAMODOrtyHWLN6VxkvrNtLSXU7uhM4lacJ/n2v9fj8n0CPnvbGA5BzKnz9rxCfBKsWwgdPtF1GRczOmtfjQDJWx9Rnm93mNlunEigH9geUvQX4J1Zv/E+AG4EZxpiizg46lvnagnrE21jz6ESf7ivnsfW1/HZJs5MtG56D0q3WYeK4b9sW2wnyJsLlT1mP374LynbaG48L2dnPKzOMdTZgXUIUuLwB+B/vTXlV11t9mXrE2xxIJzlW08hHJU0kppVZCxpqYfmvrMdT/gcSkuwLLpjTZ8Fnb8PG52H1IrjwV3ZH5Cq2N9ir6Kl2ec0rNcnKyrUN3ja9tc/Asb2QOwZO/6qNkbXinJus+/XPQu2x1tdVEdHk5SI1Db7kZXMgncSXvKrrm6xG8A+etF6YehfExehXuf9YGHIe1FdA0V/tjsZVYvQTV+3h9javFG/yqqlvgkNbrFpXWg4Mn2FzZG348ves+w+f0q4TUaTJy0VG9e/F2IEZpAa/csbxUhJ9Na9G2P6OtfCUC2K31uUz8iLoMxSOfgFbW734Q0Ugxj91FYnHrz6DV28+jz7J7vxYU7097Kvrm5olr2k2RhSmuHg4e571+MOn7I3FRdz5LVeulJwYx5jseC4cno754n1AYJhDruE841ro0Qu++A/s32B3NK6gycsljDFU1DbQ5HFvb24R4bbxyTz25UqkqR4GjIPUNnvcxIYe6XDGN63HHy60NxaX0OTlEtX1TZx+91LG/HyJ3aF0PicdMjZ35nes+61vaMN9FGjycglfHy9fdwK3qmowNG5723ritOSVNcxquK89Cvu69cUgUaHJyyV83SRSXJ68/t+HO0ko30VTjwwY8CW7w4mcr43OO4ihaj9NXi5R3WBdGuT2mte5YjV2l5800TqL5zQnT7HudwROwaAipcnLJar9Na+YGF+y03zZYyWv0tyJNkfSTnnnW5N67FkFdd129Kao0OTlEr7DxtREB9ZGwtVQy1jPpwDsy55gczDtlNIbBoy3JgTZ9V7b66uQNHm5RLdosN/9Pj2o51PPEMoTIh8fPmYM8x46artXh2jycokxA3rx6OwC5pyXZ3concf7Y1/pyW8xl6PjaKN9VGjycon+GSlcfsYAzj0l2+5QOs/etQCs8pzq7OQ1YBwkpcPhz+DoHrujcSxNXsoZmhr9faNumH0VXzkt1+aAOiA+0Wq4Bz3r2AGavFxi7a4ynn53B0V7jtodSuc4tBkaa6hJ7sd5Y0cxsE+q3RF1jL/dS5NXe2nyconlW0u47/XNvLvtkN2hdA7vIeOxXkEmm3UiX7vXjkLwdKMJU6JIk5dLVLu9h33xOgB2Jw3joaVbeW3jPpsD6qDMkyFjMNSUwQEdZaI9NHm5hL+fl1s7qRavB2CznMLjy7az5NODNgfUQSJ66NhBmrxcwtX9vOoqrTYviedoqtUVpMbJZxt9hpxr3XtrlSoymrxcwtWHjfs3gPFA7mnEJyYDUOO9ltPRBoyz7jV5tYsmL5eocfOF2b4f94Bx/pmRHN3PyyfzZEjOgIr9cMzhbXg20OTlEglxcfRIiHN98vLtnisOG+Pi4CTvsD7es6kqfJq8XOLPc85i630XMW6IQ4ZFjkSLmpc1rZsral6gh44doMlLxbaKg1C+B5J6Qt+RJCcIOek9yO6ZZHdk0aHJq91cel5ducY+q4sEJ50BcfFk9BBW3+mw4Z9b40te+4p0XPsIac3LJWY8soJpD6+wJmR1E/8howOHfA5Hei5kDIL6Cij9zO5oHMXW5CUi/UXkLRFx73xdXcAYw/aSSraXVJIU77L/R83au1zLl5j10DEitn3TReQK4ANgWDvKForIJhEpCrhdF/1IY19dowePgaT4OBLclLyMCZq8pj+8gjE/X0JZVb1NgUWZv91LzzhGws42r9uB6cCdwCntKH+xMWZXVCNyKNfOHFS2A2rLoWcu9BrgX1xV10hlXSPV9Y1kprmg4b55o/3I/7I3Fgex89/0ucYYPciPguoGl14a1LzWJeJf7EvSrujrBdC/wJqU4+CnxDXV2R2NY9iWvIwxLmtZtk+Nt5HedTUvX8fNk1o21vsuPndNX68ePaHvKPA00rNyp93ROIaTG0h+LCKrRWSLiKwUkevtDsgurr0o21fzGtiysd6XpF2TvMDfaJ9eoQcj4XJqP6+jwHbgp0AtcAXwFxE5zRhzW7ACIjIXmAvQt29fCgsLuyjUzne4xsPFeYlk9KilsLCQyspKx++feBqYuK+IOOA/O6pp3FMIQGVlJbUVtQCsWvcRdXuc+hVuqX9VOiOB1LJNjv/suoojP3ljzOUBi14QkSnAj0Tkd8aY3UHKLAIWAYwcOdJMnjy58wPtQlc1e1xYWIjj9694PaxshKzhnDf9Uv/iwsJCBvZPZ2PpAU4ZNZrJ+SfZGGQU7c+EbU+RWb3D+Z9dF3HyYWOgVVj7c6bdgagoaKV/18yxA1hw4ShG9Uvv4qA6Uc6pkJBCSu0BqC6zOxpHcFzNS0SSgBRjTHnAS74GEJc1/LSt+GgNOw9VMbBPCkOz0+wOJzq8I6cGS14XjunXxcF0gfhE6D8W9nxo7ftwF10C1UlivuYlIlnehOUzAfhHkFV93/KPOj+q2PLvzQe59plVPP2fHXaHEj3doWd9IL1IOyIxnbxEJA8oBl4NeOkCEbmk2XqTge8Cz3bHvmPVbhu/vrYcSrdBXCL0G3PCyztLq3ht4z4+KQ6sfDucXiYUETsvD3pQRIqAmd7nvkt8mteyaoAyoPkwk+uBnwB3iMgGEdkOPAncB8zpmuhji38I6ESXHDHvKwIM9DsdEnqc8PLbmw5w898+4pWPirs+ts7UvOZl9HLfttj2r9oY899hrHMAOClg2THgYe9NcbyTqmv6efn7d40P+nKKr5Nqg4v6eQH0GUpDQjqJ1aVwdDf0GWJ3RDEtpg8bVXhc10m1jfauVG8Ns9ZNnVQBRI5PqquHjm3S5OUCxy/MdkmbV1vJy4097L00eYVPk5cLuKrmdWyfNZtOjwzIDD5aUrIvebntsBGoSB9uPfB1FVEhueRfdff2qytP5/aLRpHlhnHd/bWuM6zZdYLwHTbWuG3UWJolr/1F0NQI8foTDUVrXi7QJy2JodlppCcn2h1Kx4XRv8t1o0o005DUC/oMhYZqa5ZwFZKmdRVb/Mkr+JlGgFH901l/13R3HCYHM2AcHNll/S36nW53NDFLa14ucOfLH3Pz39ZTcqzW7lA6xtMExd4LJFqZcCMxPo7MtCSS3dKvLZD2tA+LJi8XWLalhNc27qfB4/COjfuKrFl0+gyFdBdevxguf/LSRvvWaPJyAf/ZRqfXRHYss+5PntLqao1NHr75zCq+9r8fdEFQNuiXDxIPJZugvsruaGKWJi8XcM0EHJ8vt+6HTW11tfg44f3PD7N6ZxkNTZ4uCKyLJaVC7mgwHti/we5oYpYmL4drbPJQ3+QhTqBHgoM/zroK2LPKmogi7/xWVxURfy3TjWccAW33CkNUv+0i8mQ0t6fadnzmoASk2Qw7jrPrPfA0Wj/alN5tru66GYQC+c62avIKKaKuEmFM6npxB2JR7eCeQ0Zve1cbh4w+xy8Rcl9HVeB4zWuvTkQbSqT9vBa38prDT3U5kwhMGdmXXikO76C6w9ve1UZjvU+y2w8b+46E5N5QvgcOfw5ZEU8s73qRJq/NnFi7SgNGAd8AnohGUCp8OenJ/On6s+wOo2PK91qDDyalhxwGJ5Cv5lXjwusbAYiLt2qhn74E2/+tySuISNu8fmCM+SLgtskY8xLwTWB+J8So3M53ljFvojWWexgmDu/LfxWcRIbTa5ytOcU7jv32d+yNI0ZFVPMyxvy7lddqRGRUx0NSkahtaKKitpGePRKc2+4VYXsXwI+mj+ikYGLIKRdY97vehYZaSEy2N54YE2mDfbBz2AL0AS7HmgBWdaHCrYeY95d1TB+dyx+uC++QK6Z4PLCj0HocQfLqFtL7Qe7pcPBj2P0BDAuvPbC7iLTNq5DgDfMC7AWu7WhAKjI1DQ4fAvrABqgpg4zBkHly2MWOVNVzuKqOPqlJZPU8cZx71zjlAit5bX9Hk1eASNu8PgemBtwmYTXY5xljVkY3PNUWxw9E6O9VP8U6dRqmx5dtZ9rDK3lpvcsm4Qg0fLp1r+1eJ4i05vWEMWZFp0Si2sXfzyvRoaMbbXnNuo+wVuHmoaBbGHiWdRb20BY4ugd6D7I7opgRUc3LGPNoZwWi2sfRNa89a6we5Mm9YfiMiIqm+IeCdmknVZ+EJDh5kvX485Dny7olB18Mp6DZnI1OTF4feq8mG389JKVFVNSXrF03g1AwvrOOeujYgiYvh3PsnI1H98CmVyEuAc68MeLi3eawEY7399qxApoa7I0lhji0oUT5fOPsIZwzLJtR/dLtDiUya/4ApglOmwUZAyIu7r88yK097JvrPRiyR0LpVti7BoZMsDuimKA1L4cb2S+dC8f0Y2h2ZIddtqqrhHWLrcdfvqldm/BNwuHaUSUC+Wpf25bYG0cM0eSlut6G56C2HAadDQNDzxLUmnFD+vC3G8/mjou7yUUdoy6x7tcthpqjtoYSKzR5Odzzq3fz+2Wfse9ojd2hhMfjgVULrcdf/l67N5OZlsSEYdmckuOww+X2GjIBhk6E2qPw/u/sjiYmaPJyuOfW7OG3S7dxwCkzB216BQ5vh4xBMOoyu6NxDhG44OfW4w+fgooD9sYTAzR5OZyjzjbu+whevdl6POEHHZoNurymgfte28SDS7ZEKTgHGHQmjLrUmpB25YN2R2M7W5OXiPQXkbdERAcybKfjMwfF+InjI7vgr1+DhirI/zqcFXn3iObqGz08/Z+dPLd6T3Tic4qpd1nj/K9bbA1S2I3ZlrxE5ArgA6Bdo6yJyHwR2SQiG0VkvYhcHt0IncERw0BXl8FfZkFVCeRNgpmPR3QdYzCuHwY6lJxRMPYb1nj/y39pdzS2srPmdTswHXgv0oIicjvwP8Blxph8YAHwTxG5KLohxr6Yvzzo4Cb422w4/BnkjIbZz1qXvHRQirefV22DB4/TJ9uN1OTbIb4HfPICFD0Hppvtv5edyetcY8xnkRYSkd7AXcCTxpjPAYwxbwNLgd9GN8TY5vEY/zDIybE04awx1mxAf/0qPHUO7F0N6f3hmn9CckZU3iIuTkhOtL6+rh0KOpTeg2CCt+3wlXnw7BVQtsPemGxgW0OJMaa99f0LgVRgecDyZcBvRWSUMabVVlxPYx07P10V9LV+GSn+/+qHK+s4Vhv8coyE+DgG9Un1P99ZWkWoOUgy03r4hyuuqG2ktDL0mcEhmWnExVmHVPuO1lDXGPyHmdYjkYyURM5LP0iTMcQf2nT8tcpdcPDTkO8RNv9/dON9bMDTZD02TdBQYzUeN1RDzREo3W7VsA5ttSaOAEhIgS99E86d366e9K1JTUqgtqGeozUNpPWwvsrFR2uoCPGZpSYmMDjL+syaPIbPSipCbrt/Ror/MztcWcehyrqg68WJMCL3eHeNzw9VhpwINzM1iZxe1mioVXWN7DlSfcI6eyo8bDlwjLzsNHokxIfep9E/JCO+Pzkf3k/8juXw5Dl4xn2bksSBNKbl0Jiai0lIaVEkOz2ZdO/f6Wh1PUeq64PGKSIMzTre6XnPkWoaQ+xTr+RE/3hq1fVNHDwWusvOgGa/l2iI8VbeoPK99zsDlu9s9nqrySujZg95/2x7FIMs7y0ceWGul+69heOkMNb5i+/BU8eXnQlg94xZKX3grLnWLS27U94iNSmesirYW1bNgN7WD/U3b23h1aJ9Qdc/Ky+Tf3z3HMAaPvvCR98Nue1HZxdw+RlWsn35o2Lue31zyBg23Xuh//mNf17LjtKqoOt+57w87rp0NAAfF5fz9UUfBn/z997l37dOYljfnm3s00lcMPgpnhnxCmx8nrhVC+kXco9a6u29hSPcQXhSCf93EA1OTF6+X0Lgv81j3vug+UZE5gJzAcb0T+bzEB9JZrKQGG/VfCrqPVSHuA42IQ6yUo4fdZdUeULO/ZaeJKQmWtusaTQcqwvdRpGbKv7G7LIaDw0hZrNPSYBePYIf9Xs8TcTFRfswUjAiQJz/vik+CU9cMk3xPWhMSKUmpT/VqQOpSRlATUo/jCTAmk+iGkVlZSWFhYUATMptYkWDsHvLBmp2W3+LhvJ6BvYMfjIgse6Yv2xdkwm5HsCuzzZTWG61ahwsbgi9zXjj3yZAutSGXLeiZC+FhSUA7CxvCrpek8dDfFwc69esZk9q2/tU2xhHYebVZBTkk162gU17DpHFUbLMERJpeXCT0UNITrC2U91gqKgP/j0UgZzU49+t0hoPISpepCZCepI3ziZDWW3o73Z2igBR/D4YY2y9Yc0FaSJYfxHW8VlWwPLp3uXfa2sbI0aMMG62fPlyu0PoNG7eN2Pcv3/AWhOl3OHETqql3vvAoy/f88NdGItSyiZOTF4bvfdDA5bnBbyulHKxmE9eIpIlIs07Br0FVAOTA1adAmwybZxpVEq5Q0wnLxHJA4qBV33LjDFHgV8A3xeRk73rTQO+AtxmR5xKqa5n29lGEXkQq5F9sPd5kfels4wxvg4oNUAZ0OI8sTHmARGpBV4TkUagCfiqMebNLgleKWU7Ozup/ncY6xwgRHcnY81kpLMZKdVNxfRho1JKhaLJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVImryUUo6kyUsp5UiavJRSjqTJSynlSJq8lFKOpMlLKeVITpy3sUscO3aMkpISGhpCTNwYwzIyMti8OfgkqU7Xnn1LTEwkJyeHXr16dVJUyg6avII4duwYBw8eZMCAAaSkpCASenLSWFRRUUF6erjzcjtLpPtmjKGmpobi4mIATWAuooeNQZSUlDBgwABSU1Mdl7hUSyJCamoqAwYMoKSkxO5wVBRp8gqioaGBlJQUu8NQUZSSkuLIJgAVmiavELTG5S76ebqPJi+llCNp8lJKOZImLxUzFi5cyOjRoxERFi9ebHc4KsZp8uoGSkpKKCgoIDMzExGhoKCAp59+OqJtLFiwgPHjx7dYVlRUxN13383Ro0dbLK+rqyMvL4/HHnssoveYN28eb7zxRkRlVPelyasbyMnJoaioiJkzZwJW0rnhhhsi3sbgwYNbLCsqKuKee+45IXklJCQwePBgsrKyOha4Uq3QTqoqLLfeeiu33nprWOvGx8ezYsWKTo5IdXda8+rGfIeSQ4cO5c0332Tq1KkMHDiQ6dOns3fvXv963//+9xk8eDAiwq5duwC48847+dnPfgbAxRdfTEFBAbNmzeLw4cMUFBTQs2dPJk+e7N9GbW0tt99+O+PGjeNLX/oS+fn5zJs374Ram1Lh0uTVjfkOJY8cOcIHH3zAsmXL2LJlC7t37+YnP/mJf70nnniCe++9t0XZ+++/37/sjTfeoKioiBdeeIGsrCyKiopOaB87evQof/rTn3jllVdYv349q1evpqysjOuuu67zd1S5kh42RmDo7a+HfO2XV5zON8622oT+tmo3d7z8cch1dz1wif/xpY+/yyfFx4Kud/VZg/jVlfkAfLy3nNMHZrQn7DZVVFQwf/58AHr27Mn06dN56aWXovoe2dnZvP/++wwaNAiA5ORk5syZw0UXXcTBgwfJzc2N6vsp97MteYlIDvAI4PsX/TEw3xizN3Qpf9ldQLDjjduMMe9ELchuIjs7m8zMTP/zzMxMDh48GNX3SEhIYNu2bdx0000UFxeTkJBAZWUlADt27NDkpSJmS/ISkSTgbWAbcBpggD8Cy0XkDGNMZVvbMMYUdG6UJ2peY2rNN84e7K+FteW1WyaGtV5n1boAUlNTWzyPi4vD4/FE9T3efPNNLrnkEh5++GF++MMfIiIUFhYyZcoU6urqovpeqnuwq83rW0A+sMAY02iMaQIWACcD37MpJtWJnn32WXr27Mn8+fP1OkMVFXYlr6uA3caYHb4FxpgDwCbva8oBEhMTAWvMLIAlS5ZQVlYWdN26ujri4lp+3Q4cONC5ASpXsyt55QM7gyzfCZwezgZE5DcislZEtonIUhGZGdUIVZvy8vIA2Lt3LxUVFVxxxRVUVFQEXfeSSy6hvLzc37O/oqKCRx99tMtiVe5jV4N9NrAuyPJjQKqIpBhjalopXwJ8BNwFeIC5wKsicosx5vfBCojIXO969O3bl8KzYvXhAAAOFklEQVTCwpAbz8jICPkjdIKmpqYW8R86dIjLL7+cPXv2AJCfn8/cuXP5xz/+wccff0xVVRX5+fn861//4oEHHuDll1/2r/fwww/z97//nbfeeguACy+8kJtvvplvfetbnH766Vx33XVce+21pKSkMH/+fIwx5Ofns2PHDv82XnnlFWbNmsXOnTu55557eOSRR+jXrx+TJ09m1apVzJkzh3nz5pGUlMTChQsBuOuuu3j33XdPSHCB+xaJ2traVj/3WFBZWRnzMcYK8VX5u/RNReqBJcaYywKW/xX4BpDaRvIKts3XgYlAjjGmtrV1R44cabZu3Rry9c2bN3PqqadG8vYxRYeBDs4Jn2thYWGLzr1uIyLrjDHj216zbXYdNpYCwb6B6UB1pInLa5W3/GkdCUwp5Qx2Ja+NwNAgy/Ow+nuFJCIpItIzyEtN3vv4joWmlHICu5LXS8AQERnqWyAiucCpwIvNVxSRXBFpHuds4KEg2xwH1GGdsVRKuZxdyWsxVg3r1yKS4E1OD2CdbXzKt5KInAvsA54IKH+1iJzZbL3ZwOXAb8Lp4KqUcj5bzjYaY+pFZDrW5UGbsHrYfwJMDUg+lUA5sL/ZsjeBB4EnRSQR6A0cAeYZYxZ1RfxKKfvZdm2jMeYg1pnF1tbZAGQGLDsI/MJ7U0p1UzokjlLKkTR5KaUcSZOXUsqRNHkppRxJk5dSypE0eSmlHEmTl2rh0Ucf5ZVXXunS92xoaOC5555jxowZjBs3jtNOO43x48fz5z//mWADB6xbt45JkyYxZswYRo4cyW233UZtbavX4isX0uSlWrAjea1bt45rrrmG66+/nnXr1vHpp59y++238+1vf9s/vZrP9u3bmTJlCldeeSWffPIJq1atYsmSJVx//fVdGrOynyYvFRMmTJjA1Vdf7X8+a9YszjvvPB577LEWta9f/epXZGZm8oMf/ACA3r1787Of/Yznn3+eNWvWdHncyj6avFzuX//6FwUFBYgId911FwsWLGDcuHEMHDiQO++807/enj17KCgoYN++ff4yBQUFvPNO50/GdPbZZ7N8+fITlp900klUVVXR0NAAQGNjI2+88QaTJk1qMQ7+1KlTAXjxxRdP2IZyL5230eVmzpzJzJkzERGeffZZXnzxRX7961+zdOlSvvKVrzBp0iRmzJjBoEGDKCoqYujQoUyePJnFixe3ue21a9dyww03tLne+PHj/cM/ByMi/vHwm9u2bRvnnHMOSUlJgDVFWlVVlX/4aZ+srCzS09PZuHFjm7Eo99DkFa67O2/qsYjcXd7uogUFBYwbNw6AGTNm0LNnTwoLC5kxY0a7tjd+/HiKioraHU9rVq9ezcaNG1vUyEpLSwGCjqTaq1cvDh8+3CmxqNikh43dyIgRI1o879OnT9Qnl42GyspKvvOd73Dfffdx/vnnh1XGjuHMlb205hWuDtR4YkWwyWWbmppCrG2P+vp6rrrqKmbMmMFPf/rTFq9lZ2cDBJ2Ao6KigqysrC6JUcUGTV6q3aLV5uVTX1/PlVdeyejRo3nooRMHyz355JNJS0tj165dLZYfPnyYiooK8vPzw45dOZ8mL9VCYmKi/xDsiy++oLi4mAkTJgRdN5ptXr4a1/Dhw3nkkUf8y7/73e9y9913079/fxISErjoootYsWIFxhj/GUdfu9iVV14ZlViUM2jyUi3k5eWxd+9eABYtWsT+/ftDJq9oqa+vZ9asWezYsYPZs2fzl7/8xf/aypUrqaur8z+/4447mDRpEr///e+55ZZbKC8v59577+XrX/86Z511VqfGqWKMMabb3UaMGGFas2nTplZfj3XHjh3zP165cqUZO3asAUxubq656aabzJEjR8zYsWNNYmKi6dOnj7ngggv867/33ntm1KhRZsyYMebss882W7du7fR4X331VYM1FHjQ286dO1vs25o1a8z5559vRo8ebYYPH25+/OMfm5qamjbfxwmf6/Lly+0OoVMBa02Ufsda83K5iRMnBj20C3W4N2HCBDZv3tzZYbUwc+bMiM4Wjh8/nhUrVnRiRMoJtKuEUsqRNHkppRxJk5dSypE0eSmlHEmTl1LKkTR5KaUcSZNXCJGculexTz9P99HkFURiYiI1NTV2h6GiqKamJuiYYcq5NHkFkZOTQ3FxMdXV1fof2+GMMVRXV1NcXExOTo7d4ago0h72QfTq1QuAffv2+YcgdpLa2lqSk5PtDqNTtGffEhMTyc3N9X+uyh00eYXQq1cvx37ZCwsLOeOMM+wOo1O4ed9UZPSwUSnlSLYlLxHJEZG/ishW7+0FERkYZtlEEfmFiGwRkU9E5H0ROa+zY1ZKxQ5bkpeIJAFvA0nAacBooApYLiI9w9jE48BsYKIxZgzwR+BtESnopJCVUjHGrprXt4B8YIExptEY0wQsAE4GvtdaQREZCcwFHjDGHAIwxjwN7ADu79SolVIxw67kdRWw2xizw7fAGHMA2OR9rTVXAAIEzlK6DJgRZs1NKeVwdiWvfGBnkOU7gdPDKOsBdgcpm4B1CKqUcjm7ukpkA+uCLD8GpIpIijEmVBf3bKDae6gZWBYg6PxXIjIX63AToE5EPokwZifJBkrtDqKTuHnfwP37NzJaG4q1fl7SWWWNMYuARQAistYYM74D7xXT3Lx/bt436B77F61t2XXYWAqcOGe7tay6lVqXr2yqiMQHKQugc74r1Q3Ylbw2AkODLM8DPg6jbBwwKEjZRqBrZ49QStnCruT1EjBERIb6FohILnAq8GLzFUUkV0Sax/ky1pRYkwO2OQVYaow5cS74Ey2KPGRHcfP+uXnfQPcvbGLHqAneTqprsWpJ12CdPXwGOA84wxhT6V3vXGAlsMgY871m5RdiJatzjTGlIjIHeAI4xxgTnSmclVIxzZaalzGmHpgONGH17doM9AKm+hKXVyVQDuwP2MQtwD+B97xnDW8EZmjiUqr7sKXmpVQwItIf+BPwFWNMR848xxw375tdXDOqhNsv9O7g/u0SkaIgt2mdHXe4ROQK4ANgWDvLzxeRTSKyUUTWi8jl0Y2w/TqybyJS6N2vwM/uuuhHGjkRKRCRP4jIOhHZ4I31dyLSN4yyHfvdGWMcf8O6wHsD1qFkAhAP/Bn4DOgZRvmFwDagr/f5DUANUGD3vkVp/3bZvQ9hxLgKGA4str6WEZW9HasLzTDv8+lAA3CR3fsVhX0rBIbavQ+txLcF6yRbmvf5AO+ybUBKG2U79Luzfeej9Ae8EesM5MnNlvXDalP77zbKjsQ6YTAnYPmnwOt271tH98+77i679yGMGBO89xH9wIHeWCOS3Buw/HXgU7v3qyP75i3jhOR1SsCy73i/r1e1Uq7Dvzu3HDa6/ULvjuyfIxhjGttZ9EIgleCf32gRGdWhwKKgA/vmBPnGmO0By/Z57/u0Uq7Dvzu3JC+3X+jdkf0DQER+IyJrRWSbiCwVkZlRjdA++d77wL/PzoDXnezHIrLa2za0UkSutzsgH2P1HAg0AqvmtbKVoh3+3bkleWUDwTqn+i/0bqNsxBd6d7GO7B9ACfARcC7W4I+vAq+KyM1RjdIe2d77wL9PLH1+HXEU2I7Vr/E04HfAUyLyW1ujCsF72d4c4BljzLZWVu3w7y7WLsyOtk670DtGhBWjMeasgEVPiMjFwC9F5GljTG30Q7OdEz6/NhljAs+aviAiU4AficjvjDGBNRe73YV1md6P2lk+7M/NLTUvt1/o3ZH9C2WVt/xpHQksBviGjwn8+8TS5xdtq7B+u2faHUhz3sPZr2Gd5a1sY/UO/+7ckrzcfqF3u/dPRFJCNH76quuBXx6n2ei9HxqwPC/gdccRkSQRyQjyUsx9diLyTeBWrKtkSsIo0uHfnVuSl90Xene2juzfbOChINscB9RhnbF0DBHJ8l4b6/MWUE3wz2+TMWZLV8XWUUH2bQLwjyCrjvPef9T5UbVNRK7FmoNimvcsOCJyqXcAUN860f/d2d1PJEp9TZKwMvnfsdrx4rAuxWjRiROrwboJeCqg/EJgK5DtfT6H2Ouk2q79A76N1Qh6ZrNls7HO9NzbVfsQwb4uJkRfKKz/yrXAmwHLbwcO4e0HB0wjhjqptnffvD/sRuCSgGVVwP+ze3+88Vzj/a3cBlzb7Pa/wN0mxPfSu7xDvztXNNgbY+pFZDrwCFZNwgCfENmF3j/HutC7AevMVcxc6N3B/XsTeBB4UkQSsTp1HgHmGWt02ZggIg9i9Ywf7H3u+9ufZY6fjq8ByjjejwgAY8wDIlILvCYijVg/lK8aY97skuDb0IF9Ww/8BLhDRH4JpAH1wH1Yn2kseBxIJng893jvO+V3pxdmK6UcyS1tXkqpbkaTl1LKkTR5KaUcSZOXUsqRNHkppRxJk5dSypE0eSmlHEmTl1LKkTR5KaUcSZOXUsqRNHkppRxJk5eKad65Ko+JiEdE3vEue0JEjojIThG5we4YlT30wmwV80Tkq1jjWt1ojHlaRIYA/wdMMG2P2KlcSpOXcgQReQlrnK4C4I/AL40xS+2NStlJk5dyBBHphzWWWRPwmjEmZqb/UvbQNi/lCMYaXvgerCmzAicqVd2Q1ryUI3jHPy8EUrBGJD3NGFPaaiHlalrzUk7xQ6wpvy7HGnb4MXvDUXbTmpeKeSIyDHgB6+xijYh8F2vyhsuMMa/ZG52yi9a8VEwTkfuB/wD9sGaXAbjJe/9XEXnBlsCU7bTmpZRyJK15KaUcSZOXUsqRNHkppRxJk5dSypE0eSmlHEmTl1LKkTR5KaUcSZOXUsqRNHkppRxJk5dSypH+P8gKx9i0ZHKsAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "linear_convection(61)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEbCAYAAACLNQJjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNXd+PHPN8lkIwTIBsgioICAxiioFVRwAYsLj1vr41JbraK1tVK1xepPH7V1ad3r+ljb2to+j221LrWgaCVi1QeFGpFdBGSThCVknWwz5/fHvTMkw0ySSSa5c2++79drXpc5c86d72Vmvjn33nPPFWMMSinlNilOB6CUUl2hyUsp5UqavJRSrqTJSynlSpq8lFKupMlLKeVKmryUUq6U5sSbikgJ8H3gaDsGH/A28DNjzK4O2m4G9kV56SZjzNsJDlUplaTEiUGqIrIWWAVcZoypE5FhwD+xeoJHGmP87bTdbIwZ1TuRKqWSlZO7jfONMXUAxpjtwP3AWOAMB2NSSrmEI7uNQLExpimibIe9HNTbwSil3MeRnleUxAUwDjDAko7ai8gvRWSZiKwXkUUiMifhQSqlkppTPa82RCQVuAL4jTFmfQfVK4BPgNuAIDAXeFVErjPGPN7Oe8y165KZmTl55MiRCYk9GQWDQVJSvHki2cvbBt7fvvXr1+82xhQmYl2OHLA/IAiRO4CzgenGmNoutP8HcCJQZIxp6Kj++PHjzbp16+KO0y1KS0uZMWOG02H0CC9vG3h/+0RkuTFmSiLW5XiKF5HLgW8Cs7uSuGxLgf7ApIQFppRKao4mLxH5FnAjcIoxpqIT9bNEJCfKSwF7mZrI+JRSycux5CUilwLzgdOMMTvtsrPsY1OhOoNFpHWMFwIPRlndZKARWN2DISulkogjyUtELgF+DTwHnCYil9rJ7GzgILvONKzhE09ENL9IRI5pta4LgXOAX3Zjt1Mp5TJOnW18DMjEGpga6U57WQtUAV+1em2h3eZJEfEBA4FK4BpjzDM9F65SKtk4kryMMXmdqPMpkBdRVg78zH4opfowx882KqVUV2jyUkq5kiYvpZQrafJSSrmSJi+llCtp8lJKuZImL6WUK2nyUkq5kiYvpZQrafJSSrmSJi+llCtp8lJKuZImL6WUK2nyUkq5kiYvpZQrafJSSrmSJi+llCtp8lJKuZImL6WUK2nyUkq5kiYvpZQrafJSSrmSJi+llCtp8lJKuZImL6WUK2nyUkq5kiYvpZQrafJSSrmSJi+llCtp8lJKuZImL6WUKzmWvESkRER+LSLLReRTEVktIr8SkcJOtPWJyM9EZK2IrBSRD0TkhN6IWymVHJzseb0A5AEnGWOOBGYCs4D3RSSrg7aPARcCJxpjDgd+C7wlIiU9GbBSKnk4vds43xhTB2CM2Q7cD4wFzojVQETGA3OB+4wxu+y2zwIbgbt7PGKlVFJIc/C9i40xTRFlO+zloHbanQsIsDii/B3gGhHJMcbUJihGVwoaw6bddQSNCZcNG5hFpi8VgIrqBmoaW6K2zUhLYfigbACMMWzcXRfzfQr7Z5Cb6QOgqr6Z3XWNMeseUpgT/vfWvfU0BYJR6+Vm+ijsnwFAQ3OA7fv8bV7f1xi9HfV7oX5PzPdPKEnZv0xNh7QM65GeAyK9E4NyLnlFSVwA4wADLGmnaTEQBLZElG/C2p6JwEeJiNGtfruyiX+9WdqmbOH1JzJhaC4Av3xzHS8u3xa17VEjB/LytdMAaAkaTn3w3Zjv88sLivnmlBEAvFK2nf96bVXUemkpwoZ79nemr/rDMtburIla97LjD+au/zgcgFU7qjn/qQ/avC5A/phdnDi21aHRPV/AU9OgpW2i63XpOTBolPUYPR2O/hb4OjoCorrKyZ5XGyKSClwB/MYYs76dqgVAvTEmEFFebS/zY6x/LtbuJoWFhZSWlnYv4CS2eV8zIORnCj67k1C2fBnlOdYT/94mhmRH7yH4mmrC/zctQROzHsCXG9ZRWvsFANt3tMSsm5JCm//vjEBDzLo1u3ZQWrrbWn91oE29mmZDXTO8/q8yAtt94fJx6x7noBY/Tb4BtKT1ixlv4hjEGCBISrCFlGAzKcFGUptqoXyl9Vj7Oo3/vIctI8/nq6GzCKZmdGrNtbW1nv5uJpKYVrsWThKRO4Czgent7faJyCLgeGNM/4jyq4BngDOMMQvbe6/x48ebdevWdT/oJPW1uxaws97w9g0ncWhR/44buMT9b67licVfcOPMcVx36lirsGYnPHIEBJrhB8ug4FBngjMG/JVQuRl2rYX/ewp2rrBeyxsDV70DWe0dDbGUlpYyY8aMHg3VSSKy3BgzJRHrcvqAPQAicjnwTWB2J45X7Qay7Z5aa6FfaS8d+EheTfZhoYy0yP8idzty+EBmDE9j0rDc/YX/9xQEmmDCWc4lLrCOdWXnwbCjoeRiuHoJXPQC5I+FvRvhjVuci82jHE9eIvIt4EbgFGNMRSearMCKe0RE+WigBViT2Ajd5+riDJ69bEr4wLdXzJo0hO8cnsEphw22ChqqYNlvrX9Pm+dcYNGIwPjZVgJLy4RP/wfWtbtDoOLkaPISkUuB+cBpxpiddtlZ9vGpUJ3BItI6zpexDurPiFjdycAiY0z0I8F9yPi8VE6bODh8dtGzlj8HjdVw8AkwPCF7IolXcCicerv1779fb50VVQnh5Aj7S4BfA88Bp4nIpXYyOxs4yK4zDWv4xBOhdsaYdVjHtn4qIgV2vSuAQ4Bbe3MbVO+q8jezuSrAxl210NIIHz5pvXBCkvW6Ih13DYz4GtSWwxs3Ox2NZzjZ83oMyMQamPp8q8fcVnVqgSrgq4i21wF/xRqNvxK4CphljCnr6aCTXXMgyJ/XNfHYPz93OpSEe+/zXdzxYQMPLFoHa/8BtTuhaBIceprTobUvJRXOeRLSsmDFn2Hbcqcj8gQnx3nldaLOp1iXEEWWNwP/z36oVuobAyzc1My/vtq4/4ycR2TZu8ENzUH4yv47NeFsdwwMzT8Ejr4MPvpvWLcAhk92OiLXc/yAvUqshhZr+JsXj3eFtsnfFIAK+7zM4IkORhSncbOs5Ya3nI3DIzR5eUxDcyh5ee+jDW1TQ0sAKtZahYUTHIwoTgefYO06fvUp1JQ7HY3ree8b3sf57eSV5eGelzTWQNUW67rCvDEORxUHXyaMPsn694a3nY3FAzR5eUxDszVC1cu7jUOavrQKCsZDatJc4dY5Y2day88XORuHB2jy8pj9u43eTV7DmjZbBUUu2mUMCZ0Z/WIxBKLP7KE6R5OXx6SIUJAlDMnNdDqUhCvMyeD24zP5/iR7QpKiw5wNqCvyRluXDDVWwbY+PflJt2ny8phjR+fxwPRsfnXRUU6HknDpaSmMGZDKwNoNVkGRi840tjbWPuv4uZ517A5NXsp9Qmca3bjbCDDW3nXU5NUtmryUq7y8cjfU7sT4+sGAkU6H0zUHTwNfNpR/BtU7Oq6votLk5THPf7iZ7/+zjgcXeXO+sr07NwMQyB9nzXLoRmkZ1kyrAF+842wsLubST1/FUtPYQl0zMeeId7sJqdb01U354x2OpJtGHmctd650Ng4X0+TlMeFxXh6biDBkrFjJyz/Q5ckrdGXArrXOxuFimrw8JjTOKyvdm8nrULYCUDvA5RedF46zlru8uXvfGzR5eUx4kGqaBz9aYxiD1fOq7u/glM+JMPBga4bVmh3WjLAqbh78hvdtXh5hT20FA6ilymRT4ytwOpruSUmFArv3uKu9m2WpWDR5eYzfPublyd3GitUAbE8fTWa6y65pjKbQvkJAj3t1iQe+Aaq1OUceRHr9LiYOze24stvYc3hNLD4WRnZ8G7GkV2ifdNDk1SWavDxm5sTB+CrSGTvYO/drDNttH9wudOE1jdFoz6tbdLdRuUflZgCaB46m2Qvj2MLJS884doUmL49ZvLaCZTtbqG5odjqUxLOT1+m/38IfPvzS2VgSYdBoSPFB1VZo7PN37IubJi+PuXvBGh4va2RnVYPToSRWoAX2bSWIsN0UhM+qulpq2v4zjrv1jGO8NHl5zP5xXh4721i9DUyAqtQ8Gkn3RvKCVgftddcxXpq8PCacvNI99tHau4yVaYMBPJS89KB9V3nsG648O4e9nbz2+azk5fdM8tKeV1dp8vKYBq/ePchOXlXpRcD+JO162vPqMk1eHtIcCNISNKQI+FI99tHayasmYwjgoZ5X3iEgqVD5JTTVOx2Nq+ggVQ8J9bq8drgLCCevfgOH8Ivzj+Dg/H7OxpMoaemQf4h1tnHP505H4yqavDwkJyONVXeezjvvvud0KIlnJ6/MQQdx4TEunf45lsLxVvLatQ4ocjoa1/Di3+g+S0Tol5FG/3RxOpTE8u8DfyX4smn2DXA6msTT415d4mjyEpGhIvKGiBgn41BJbp89mn7QKGqbrXn6X1q+zdGQEqrAPuO4W3cb4+HYbqOInAs8DMR9HYuIlGL1r5siXnrIGPOH7kfnThsqarj15ZX0DzQyY4bT0SSQvcvIoFFUNRlue3UVYwr7cf7k4Y6GlTB5o61l5WYY7GgkruLkMa+bgZnArUBXpsU8wxizOaERudzeumaWbtrL2IEeOxrQKnml2yMkGr0yVAJg0ChrWemB6zV7kZPJa5oxpkXEY8dnHBQ+2+ixIV5tktc+6/vimRH2ANn5kJ4DjVWkNdc6HY1rOPYn2hjT4tR7e5U/nLw89gehdfKyE7NnxnkBiIR7X5kNO52NxUXcvH9xg4h8JCJrRWSJiFzudEBOC/VGfG7+VKNplbxC29bQHMAYD53nsZNXlr/c2ThcxK3jvPYBG4CfAg3AucAfRWSSMeamaA1EZC4wF6CwsJDS0tJeCrX3lG2zzn2kBFu8s30mwEmVX5ICLFmxmQZ/M6kiBAz8c3EpaSne6GUeUpvKCCClaot3Prse5srkZYw5J6LoRRE5GfiRiPzKGLMlSptngGcAxo8fb2Z46nSc5csPNsPKVWRn+vDM9u3bAu8GoP9QTjr1dEpLS+mX0UhL0HDs1BPIzfQ5HWFiZH8O215lQLCSCV757HqYK5NXDEuBa4FjgAOSV18wMj+bs4qHMji4x+lQEmfvJmsZOiMHfPpfs/DciR495hU31yUvEUkHsowxkXfqDB3B9dq5tk47eXwRJ48v8tZuR6vjXSGeS1ygx7y6IOkP7YpIvp2wQqYCf4lSdbK9/KTno1K9Jkry8qQBIwAhs6HCmvJadSipk5eIjAa2A69GvHSqiJzZqt4M4GrgeWNMn73Gory6gU276/C3eOgsXJTkNe+FTzjlwVJW76h2JKQe4cuE3IMQgtaU16pDjiUvEblfRMqAOfbzMvvRupflB/YCO1qV/Rv4CXCLiHwqIhuAJ4GfA1f0TvTJ6eG31nPyA6Us/cpDf7n3brSWrZLX9n1+Nu6q894dksIj7Tc7GYVrOHbMyxjz407U2QkcFFFWDTxkP1Qrfq+N8woG91+sXDAuXBya4tpTo+zBSl5fvq/Jq5O88jVXtL48yCMHtKu3QXMd9CuC7LxwsaeTF2jy6iRNXh4SmtfdM9c2hm5KEbpJhS0rnLw8dHE2aPKKkyYvDwlf2+iRUefhyflCk/XZMu39Yk9d3wiavOKkyctDGr02q0Q4ecXqeWny6stcN0hVxbZ/t9ErPa/QbmPbntfXxuSTkiJMGJrrQFA9qF8hgZQMUv2V1tTXWQOdjiipafLykAe/eSRV/mZqNn/mdCjdZ0zM5DX7iKHMPmKoA0H1MBH8WUPIqfvSmvpak1e7dLfRQw4fNoBphxaQmeaBnlfNV9BYDVl50K/A6Wh6TUOmPQ+07jp2SHteKjm1PlgfcS3j7tpGNu2uY1B2OocW5TgQXM/xZ2ny6izteXnInX9fxb0L1hAIeuDyoBjDJADeXl3ON57+kGeWfNHLQfW8hkzrjuCavDqmPS+PCAQNv3t/MwBfOz3b2WASIcYwCWg9SNVj47wAf5Ymr87SnpdHNLZYwwYyfSnemDKmnZ6XZ8d50Sp57d7gbCAukNDkJSJPJnJ9qvP8TdYPOTQGytWMgYo11r/b7Xl5L3k1ZA6F1HSo2gKNNU6Hk9Ti2m0Ukcs6qHJGN2JR3dDQYu1CZXohedXtgoZ9kDEA+g854OXQNnrq3o02k5JqXYRevtLqfQ6f4nRISSveY17PtfOaB44Su1eoF+KJnlfrkfVRdoFDycuLu40AFE2wklfFGk1e7Yg3ea3hwN5VP+Aw4GLgiUQEpeIX2m3M8KWyf0Zsl2rneBd4+PKgkKIJ1jK066yiijd5/dAYE+2e5KtFZCHwArC4+2GpeKWmCGOLchiVnw24/FhJO2caAQ7Oz+bNeSeRk+nRk+WFoeS12tk4klxcn74x5p/tvOYXkejfNtXjJgzN5a0bpgO4/wYcMS4LCsn0pTJ+SP9eDKiXac+rU+I9YH9StGJgEHAO1g1gleo6Y/b3OGLsNnrewIPBlw21O6F+b5uJGNV+8fa7S4l+YF6AbcCl3Q1I9XHlK6F+D/QfCgOGR60SCBp+/NdPaQoEefzio3s5wF6QkmIl7h2fWLvQB091OqKkFG/y+gK4MqIsAFQAXxhjPHoENfm9uHwbt7z8Gd+YPJyZg5yOphs+X2QtDz0t6plGgBSBV8q2EzTwyIVB0lI9ONa6aKKVvCpWa/KKId7k9YQx5t0eiUR1i785QFNL0P3jVT5/21qOnRmzioiQ6UulvilAQ0uQHE8mr9Bxr7XOxpHE4vrUjTGP9FQgqnsavTDOy78Pti6FlDQYM6PdqqHtDA0R8Rw9aN8hD/7J6ptCP+JMN9/3bONiMAEY8TXIHNBuVS9fIgS0HS5hXN+f7hEu/qar1hpCF2anubjn9flb1rKdXcaQUJIOXZDuObkHWZdH+fdCbYXT0SQlTV4e4W+yrvPLcuvdN4JB2NDx8a4QL0+LA1gnK0K7jrt01zEaTV4eEep5Zbj1mNfOFVBbDrnDrDNtHSgePpDjx+STnubhr3CRPUhXj3tF5dHrK/qes4qHckhhDlMOHkT5uk1OhxO/DfYuYztDJFq797wjejigJBBK4nqZUFSavDxi6iEFTD3EulFF+TqHg+mK8PGuWc7GkUz0jGO7PNznVq5RtQ22fQwpPhgzvVNNmgNBqvzN3j3bCK16XmuhpcnZWJKQJi+PeGdtOS9/so3dtY1OhxIfY+Dv14MJwmFnQkbnLrj+6d8+48g7F/Hapzt6OEAH9Suwhkw01cDKl5yOJuk4mrxEZKiIvCEiOpClmx57ZwM/+vOnfLmnzulQ4vPJ89ZZxsyBMPsXnW4WGirh6Z4XwNTrrOX7j1pnZFWYY8lLRM4FPgQO6WL7eSKyWkRWiMi/ReScxEboLuHJCN00zmvfVnjzVuvfZzwQdcrnWDw/IWHIEd+A/gdZwyVC130qwNme183ATOD9eBuKyM3A/wPONsYUA/OBv4rI7MSG6B6NLS4b5xUMwmvXWXfFPuwsOOKCuJqHp4Ju8nhvJC0djr/W+vf7jzobS5JxMnlNM8Z8Hm8jERkI3AY8aYz5AsAY8xawCHggsSG6R6gHkvQ34AgGYdXL8PQ063KgrDw46+FODY9oLTxI1asj7Fub/B1rtP2WD2DrR05HkzQcGyphjGnpYtOvA9kcON30O8ADInKYMabdS/GDLY1sWrU06mtDBmSFd0n21DZS3dActV5aagojBu2/ueum3XXEugdJXr8MBmT5AKhpaGF3bew5Gw/O60dKivVD3rHPH/Pyl34ZPor6ZwBWryu3ej25Av0q19KvdjOUr4r5HgnT5po7Yz+3l8EABBqhpdHqXVV+ad1I9cv390/znDsc/uNxyCmK+61Dyavav//z2VZZT21j9K9Vti+NkfnW5xUIGj6viD1V9kEDs8jNtD6vPbWN7IpxEiRFhHGD959g2FBRS0uM41J52ekU5WYCUNfYwtbK+qj1ttYEaWoJhgffWttkKJh4KQWfPEHtW/ey67ibAchMS2PowMzwNm3ZG/t4Z2H/THIyrJ/7vvomKuujn70UEUbl9ws/37K3nkCMbcrN9JGfY30H65sClFf7Y77/8EHZ+BI8+4cbx3kV28vIkZibWr3ebvIa4N/K6L92PJ4o3350xuhO1utvPzrjoE7WywDezLCf/B6OAVjWycZOyB0OJ94AR10KaRkd148i9AdmbFFOuOyeBWtY8NnOqPWnHZrPn678GgC1DS18/ZH3Yq77iYuP5szioQD8Zdk2fvFG9K9TbmYaK+44Pfz88uc+Yuve6D/gq6eP4aezrXFbn2zZx6W/if7HE2D61IZwog1tUyGT+FeGj5wt75Cz5Z0D2qTS+e/gQPvRGSM7WS87jvdPFDcmrwJ7Gfmns9peRs03IjIXmAtw+NBMvmBE1JXnZQq+VKvnU9MUpD56x4u0FMjP2v+XpKIu9lxa/dOFbJ+1Tn+Lobox9snVwdkS3oXa6w8S69K9rDTIzbDevyVo2OM3ZKTBwIwUgsEAKSlO7D4KIBgBI6kY8RFM8RFIzaQhczD+rMH4s4ZSOehITJ0P/vVh3O9QW1tLaWkpGf4g4wel0L96E6Wl1j1hAjWNDM+JvvuZ2lAVntu/vtnErAewcd0qSvdaI30rtjXHrJuZFmxzv4CBKU2YGHWry7dRWloOwIbKQMx1BoJBPv7o/9hof7f2b9MgngpezFnBdwi1TE+FQZlWPWMMFfWxv1cDMoTMNKtlXbOhtil63RSBwuz93+vd/iCBGN/Bfj4hJ91aZ1PAUNkQ+/0LsoTUFAE+i1knbsYYRx9Y94I0cdR/Bmv/LD+ifKZd/r2O1jFu3DjjZYsXL3Y6hB7j5W0zxvvbBywzCcodbhykutteRu59hZ7v6cVYlFIOcWPyWmEvR0WUj454XSnlYUmfvEQkX0TSWxW9AdQDMyKqngysNh2caVRKeUNSJy8RGQ1sB14NlRlj9gE/A74vImPseqcBpwM3ORGnUqr3OXa2UUTuxzrIPtJ+Xma/dKwxJjQIxQ/sBdpcfWuMuU9EGoDXRaQF6/Zr3zDGLOyV4JVSjnNykOqPO1FnJzGGOxnrTkZ6NyOl+qik3m1USqlYNHkppVxJk5dSypU0eSmlXEmTl1LKlTR5KaVcSZOXUsqVNHkppVxJk5dSypU0eSmlXEmTl1LKlTR5KaVcSZOXUsqVNHkppVxJk5dSypU0eSmlXMmN923sFdXV1VRUVNDcHOPGjUlswIABrFmzxukwekRXts3n81FUVERubm4PRaWcoMkriurqasrLyxk2bBhZWVmIxL5BaTKqqamhf//O3pfbXeLdNmMMfr+f7du3A2gC8xDdbYyioqKCYcOGkZ2d7brEpdoSEbKzsxk2bBgVFRVOh6MSSJNXFM3NzWRlZTkdhkqgrKwsVx4CULFp8opBe1zeop+n92jyUkq5kiYvpZQrafJSSePpp59m4sSJiAjPPfec0+GoJKfJqw+oqKigpKSEvLw8RISSkhKeffbZuNYxf/58pkyZ0qasrKyMO+64g3379rUpb2xsZPTo0Tz66KNxvcc111zDggUL4mqj+i5NXn1AUVERZWVlzJkzB7CSzpVXXhn3OkaOHNmmrKysjDvvvPOA5JWWlsbIkSPJz8/vXuBKtUMHqapOufHGG7nxxhs7VTc1NZV33323hyNSfZ32vPqw0K7kqFGjWLhwIaeccgrDhw9n5syZbNu2LVzv+9//PiNHjkRE2Lx5MwC33nort99+OwBnnHEGJSUlXHDBBezZs4eSkhJycnKYMWNGeB0NDQ3cfPPNTJ48maOPPpri4mKuueaaA3ptSnWWJq8+LLQrWVlZyYcffsg777zD2rVr2bJlCz/5yU/C9Z544gnuuuuuNm3vvvvucNmCBQsoKyvjxRdfJD8/n7KysgOOj+3bt4/f/e53vPLKK/z73//mo48+Yu/evVx22WU9v6HKk3S3MQ6jbv5HzNfuOfcILj7OOib0P0u3cMvLn8Wsu/m+M8P/Puux91i5vTpqvYuOHcG95xUD8Nm2Ko4YPqArYXeopqaGefPmAZCTk8PMmTP529/+ltD3KCgo4IMPPmDEiBEAZGZmcsUVVzB79mzKy8sZPHhwQt9PeZ9jyUtEioCHgdCf6M+AecaYbbFbhdtuBqLtb9xkjHk7YUH2EQUFBeTl5YWf5+XlUV5entD3SEtLY/369Vx77bVs376dtLQ0amtrAdi4caMmLxU3R5KXiKQDbwHrgUmAAX4LLBaRo4wxtR2twxhT0rNRHqh1j6k9Fx83MtwL68jr153YqXo91esCyM7ObvM8JSWFYDCY0PdYuHAhZ555Jg899BDXX389IkJpaSknn3wyjY2NCX0v1Tc4dczr20AxMN8Y02KMCQDzgTHA9xyKSfWg559/npycHObNm6fXGaqEcCp5nQ9sMcZsDBUYY3YCq+3XlAv4fD7AmjML4M0332Tv3r1R6zY2NpKS0vbrtnPnzp4NUHmaU8mrGNgUpXwTcERnViAivxSRZSKyXkQWicichEaoOjR69GgAtm3bRk1NDeeeey41NTVR65555plUVVWFR/bX1NTwyCOP9FqsynucOmBfACyPUl4NZItIljHG3077CuAT4DYgCMwFXhWR64wxj0drICJz7XoUFhZSWloac+UDBgyI+SN0g0Ag0Cb+Xbt2cc4557B161YAiouLmTt3Ln/5y1/47LPPqKuro7i4mNdee4377ruPl19+OVzvoYce4s9//jNvvPEGAF//+tf5wQ9+wLe//W2OOOIILrvsMi699FKysrKYN28exhiKi4vZuHFjeB2vvPIKF1xwAZs2beLOO+/k4YcfZsiQIcyYMYOlS5dyxRVXcM0115Cens7TTz8NwG233cZ77713QIKL3LZ4NDQ0tPu5J4Pa2tqkjzFZSKjL36tvKtIEvGmMOTui/E/AxUB2B8kr2jr/AZwIFBljGtqrO378eLNu3bqYr69Zs4YJEybE8/ZJRaeBjs4Nn2tpaWmbwb1eIyLLjTFTOq7ZMad2G3cD0b6B/YH6eBOXbandflJ3AlNKuYNTyWsFMCpK+Wis8V4xiUiWiOREeSlgL1O7F5pSyg2cSl73We3EAAAMxklEQVR/Aw4WkVGhAhEZDEwAXmpdUUQGi0jrOC8EHoyyzslAI9YZS6WUxzmVvJ7D6mH9QkTS7OR0H9bZxqdClURkGrADeCKi/UUickyrehcC5wC/7MwAV6WU+zlyttEY0yQiM7EuD1qNNcJ+JXBKRPKpBaqAr1qVLQTuB54UER8wEKgErjHGPNMb8SulnOfYtY3GmHKsM4vt1fkUyIsoKwd+Zj+UUn2UTomjlHIlTV5KKVfS5KWUciVNXkopV9LkpZRyJU1eSilX0uSl2njkkUd45ZVXevU9m5ub+d///V9mzZrF5MmTmTRpElOmTOH3v/890SYOWL58OdOnT+fwww9n/Pjx3HTTTTQ0tHstvvIgTV6qDSeS1/Lly7nkkku4/PLLWb58OatWreLmm2/mO9/5Tvj2aiEbNmzg5JNP5rzzzmPlypUsXbqUN998k8svv7xXY1bO0+SlksLUqVO56KKLws8vuOACTjjhBB599NE2va97772XvLw8fvjDHwIwcOBAbr/9dl544QU+/vjjXo9bOUeTl8e99tprlJSUICLcdtttzJ8/n8mTJzN8+HBuvfXWcL2tW7dSUlLCjh07wm1KSkp4++2evxnTcccdx+LFiw8oP+igg6irq6O5uRmAlpYWFixYwPTp09vMg3/KKacA8NJLLx2wDuVdet9Gj5szZw5z5sxBRHj++ed56aWX+MUvfsGiRYs4/fTTmT59OrNmzWLEiBGUlZUxatQoZsyYwXPPPdfhupctW8aVV17ZYb0pU6aEp3+ORkTC8+G3tn79eo4//njS09MB6xZpdXV14emnQ/Lz8+nfvz8rVqzoMBblHZq8OuuOnrv1WFzuqOpy05KSEiZPngzArFmzyMnJobS0lFmzZnVpfVOmTKGsrKzL8bTno48+YsWKFW16ZLt37waIOpNqbm4ue/bs6ZFYVHLS3cY+ZNy4cW2eDxo0KOE3l02E2tpavvvd7/Lzn/+ck046qVNtnJjOXDlLe16d1Y0eT7KIdnPZQCAQo7YzmpqaOP/885k1axY//elP27xWUFAAEPUGHDU1NeTn5/dKjCo5aPJSXZaoY14hTU1NnHfeeUycOJEHHzxwstwxY8bQr18/Nm/e3KZ8z5491NTUUFxc3OnYlftp8lJt+Hy+8C7Yl19+yfbt25k6dWrUuok85hXqcY0dO5aHH344XH711Vdzxx13MHToUNLS0pg9ezbvvvsuxpjwGcfQcbHzzjsvIbEod9DkpdoYPXo027ZtA+CZZ57hq6++ipm8EqWpqYkLLriAjRs3cuGFF/LHP/4x/NqSJUtobGwMP7/llluYPn06jz/+ONdddx1VVVXcdddd/Od//ifHHntsj8apkowxps89xo0bZ9qzevXqdl9PdtXV1eF/L1myxBx55JEGMIMHDzbXXnutqaysNEceeaTx+Xxm0KBB5tRTTw3Xf//9981hhx1mDj/8cHPccceZdevW9Xi8r776qsGaCjzqY9OmTW227eOPPzYnnXSSmThxohk7dqy54YYbjN/v7/B93PC5Ll682OkQehSwzCTod6w9L4878cQTo+7axdrdmzp1KmvWrOnpsNqYM2dOXGcLp0yZwrvvvtuDESk30KESSilX0uSllHIlTV5KKVfS5KWUciVNXkopV9LkpZRyJU1eMcRz6l4lP/08vUeTVxQ+nw+/3+90GCqB/H5/1DnDlHtp8oqiqKiI7du3U19fr3+xXc4YQ319Pdu3b6eoqMjpcFQC6Qj7KHJzcwHYsWNHeApiN2loaCAzM9PpMHpEV7bN5/MxePDg8OeqvEGTVwy5ubmu/bKXlpZy1FFHOR1Gj/Dytqn46G6jUsqVHEteIlIkIn8SkXX240URGd7Jtj4R+ZmIrBWRlSLygYic0NMxK6WShyPJS0TSgbeAdGASMBGoAxaLSE4nVvEYcCFwojHmcOC3wFsiUtJDISulkoxTPa9vA8XAfGNMizEmAMwHxgDfa6+hiIwH5gL3GWN2ARhjngU2Anf3aNRKqaThVPI6H9hijNkYKjDG7ARW26+151xAgMi7lL4DzOpkz00p5XJOJa9iYFOU8k3AEZ1oGwS2RGmbhrULqpTyOKeGShQAy6OUVwPZIpJljIk1xL0AqLd3NSPbAkS9/5WIzMXa3QRoFJGVccbsJgXAbqeD6CFe3jbw/vaNT9SKkm2cl/RUW2PMM8AzACKyzBgzpRvvldS8vH1e3jboG9uXqHU5tdu4Gzjwnu1WWX07va5Q22wRSY3SFkDv+a5UH+BU8loBjIpSPhr4rBNtU4ARUdq2AL179willCOcSl5/Aw4WkVGhAhEZDEwAXmpdUUQGi0jrOF/GuiXWjIh1ngwsMsYceC/4Az0Tf8iu4uXt8/K2gW5fp4kTsybYg1SXYfWSLsE6e/gb4ATgKGNMrV1vGrAEeMYY871W7Z/GSlbTjDG7ReQK4AngeGNMYm7hrJRKao70vIwxTcBMIIA1tmsNkAucEkpctlqgCvgqYhXXAX8F3rfPGl4FzNLEpVTf4UjPS6loRGQo8DvgdGNMd848Jx0vb5tTPDOrhNcv9O7m9m0WkbIoj9N6Ou7OEpFzgQ+BQ7rYfp6IrBaRFSLybxE5J7ERdl13tk1ESu3tivzsLkt8pPETkRIR+bWILBeRT+1YfyUihZ1o273fnTHG9Q+sC7w/xdqVTANSgd8DnwM5nWj/NLAeKLSfXwn4gRKnty1B27fZ6W3oRIxLgbHAc9bXMq62N2MNoTnEfj4TaAZmO71dCdi2UmCU09vQTnxrsU6y9bOfD7PL1gNZHbTt1u/O8Y1P0H/gVVhnIMe0KhuCdUztxx20HY91wuCKiPJVwD+c3rbubp9dd7PT29CJGNPsZVw/cGAg1owkd0WU/wNY5fR2dWfb7DZuSF6HRpR91/6+nt9Ou27/7ryy2+j1C727s32uYIxp6WLTrwPZRP/8JorIYd0KLAG6sW1uUGyM2RBRtsNeDmqnXbd/d15JXl6/0Ls72weAiPxSRJaJyHoRWSQicxIaoXOK7WXk/8+miNfd7AYR+cg+NrRERC53OqAQY40ciDQOq+e1pJ2m3f7deSV5FQDRBqeGL/TuoG3cF3r3su5sH0AF8AkwDWvyx1eBV0XkBwmN0hkF9jLy/yeZPr/u2AdswBrXOAn4FfCUiDzgaFQx2JftXQH8xhizvp2q3f7dJduF2YnWYxd6J4lOxWiMOTai6AkROQO4R0SeNcY0JD40x7nh8+uQMSbyrOmLInIy8CMR+ZUxJrLn4rTbsC7T+1EX23f6c/NKz8vrF3p3Z/tiWWq3n9SdwJJAaPqYyP+fZPr8Em0p1m/3GKcDac3enf0m1lne2g6qd/t355Xk5fULvbu8fSKSFePgZ6i7HvnlcZsV9nJURPnoiNddR0TSRWRAlJeS7rMTkW8BN2JdJVPRiSbd/t15JXk5faF3T+vO9l0IPBhlnZOBRqwzlq4hIvn2tbEhbwD1RP/8Vhtj1vZWbN0VZdumAn+JUnWyvfyk56PqmIhcinUPitPss+CIyFn2BKChOon/3Tk9TiRBY03SsTL5n7GO46VgXYrRZhAn1gHrAPBURPungXVAgf38CpJvkGqXtg/4DtZB0GNalV2Idabnrt7ahji29TlijIXC+qvcACyMKL8Z2IU9Dg44jSQapNrVbbN/2C3AmRFldcAfnN4eO55L7N/KTcClrR7/DdxhYnwv7fJu/e48ccDeGNMkIjOBh7F6EgZYSXwXev8X1oXezVhnrpLmQu9ubt9C4H7gSRHxYQ3qrASuMdbssklBRO7HGhk/0n4e+r8/1uw/He8H9rJ/HBEAxpj7RKQBeF1EWrB+KN8wxizsleA70I1t+zfwE+AWEbkH6Ac0AT/H+kyTwWNAJtHjudNe9sjvTi/MVkq5kleOeSml+hhNXkopV9LkpZRyJU1eSilX0uSllHIlTV5KKVfS5KWUciVNXkopV9LkpZRyJU1eSilX0uSllHIlTV4qqdn3qqwWkaCIvG2XPSEilSKySUSudDpG5Qy9MFslPRH5Bta8VlcZY54VkYOBvwNTTcczdiqP0uSlXEFE/oY1T1cJ8FvgHmPMImejUk7S5KVcQUSGYM1lFgBeN8Ykze2/lDP0mJdyBWNNL3wn1i2zIm9Uqvog7XkpV7DnPy8FsrBmJJ1kjNndbiPladrzUm5xPdYtv87Bmnb4UWfDUU7TnpdKeiJyCPAi1tlFv4hcjXXzhrONMa87G51yiva8VFITkbuBfwFDsO4uA3CtvfyTiLzoSGDKcdrzUkq5kva8lFKupMlLKeVKmryUUq6kyUsp5UqavJRSrqTJSynlSpq8lFKupMlLKeVKmryUUq6kyUsp5Ur/H4T7S+lFwR2SAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "linear_convection(71)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So far so good—as we refine the spatial grid, the wave is more square, indicating that the discretization error is getting smaller. But what happens when we refine the grid even further? Let's try 85 grid points." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEbCAYAAACLNQJjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXl8XFX5/99Pksnepk26AF1oWVrWEGkpP8uWAi2bVjYFFVFqrYBbBbS4gBQFUZRNcamo+BUXFBAQKSDStAhYaKGULrSUtrRN9zX7Ouf3x70zmZncmcyaO3f6vF+vvG7mzLl3npO595PnPOc554gxBkVRFK+R57YBiqIoyaDipSiKJ1HxUhTFk6h4KYriSVS8FEXxJCpeiqJ4EhUvRVE8SYEbHyoiNcCXgJNtG3zAi8D3jTG7+jh3I7Df4a2bjDEvptlURVGyFHEjSVVE3gVWAlcbY5pFZATwHyxP8CRjTGuMczcaY8b0j6WKomQrbnYb5xhjmgGMMfXA3cDRwIUu2qQoikdwpdsIVBtjOiLKttrHwf1tjKIo3sMVz8tBuADGAQZY1Nf5IvJjEVkiImtF5AURmZ52IxVFyWrc8rzCEJF8YAbwW2PM2j6q7wTeAm4B/MAs4CkR+Yox5ucxPmOWXZfi4uIJo0ePTovt2Yjf7ycvLzcHkr3etoKuFkpat9JVUEpBVwsgNA44Mvi+19vXF2vXrt1tjBmajmu5ErDvZYTIbcBHgbOMMU1JnP8v4AxgmDGmra/648ePN2vWrEnYTq9QV1dHbW2t22ZkBM+3bc18+MuVcPR58N7zVtn39oMIkAPt6wMRWWqMmZiOa7ku8SJyDfAJ4IJkhMtmMTAAOD5thilKJvB3W8e8ApC88DIlIVwVLxH5DHAjcLYxZmcc9UtEpNzhrcC3n59O+xQl7fi7rGNeHoh9uxoVr2RwTbxE5CpgDnCuMWa7XfYROzYVqDNcREJtvAL4qcPlJgDtwKoMmqwoqWNCPK88W7zU80oKV8RLRD4N/AZ4GDhXRK6yxeyjwGF2ndOw0icejDj9kyJySsi1rgAuBn6cQrdTUfqHgFBJvnpeKeLWaOPPgGKsxNRI5trHJuAAsC3kvfn2Ob8QER8wCNgHXGuMmZc5cxUlTQRjXvnqeaWIK+JljKmMo87bQGVE2Q7g+/aPongP4xCwN3737PEwro82KspBRSBgL3nqeaWIipei9CdhqRIa80oFFS9F6U/CYl4F4WVKQqh4KUp/4pgq0eWePR5GxUtR+pOwVIlAwF49r2RQ8VKU/iQ0wz7oeeloYzKoeClKf2I0YJ8uVLwUpT8J7TZqqkRKqHgpSn+iqRJpQ8VLUfqTYMwr34p7gXpeSaLipSj9iQnJ81LPKyVUvBSlP3GMeeloYzKoeClKf+KUYa+eV1KoeClKf+KUKqEZ9kmh4qUo/UlwVQlNlUgVFS9FySR+P2z6H3S22q9DA/Y6PSgVVLwUJZOsfgp+dx68bG+9YJxWUtWAfTKoeClKJtmx0jo2bLWOuoZ92lDxUpRMcqDeOnbZeyH7nZbEUfFKBhUvRckkBzZbx6526xiaYa+eV0qoeClKJmmI8LzCYl46PSgVVLwUJVMYE9JtjPC8NOaVMipeipIpmndDty1awZiXPbIYFvPS0cZkUPFSlEwRiHdBlG6jTg9KBRUvRckUB7b0/B6r26jTg5JCxUtRMkUgWA8OqRIasE8VFS9FyRSxPC9NlUgZFS9FyRRh4hWIeTkF7FW8kkHFS1EyhaPn5TQ9SEcbk0HFS1EyRWTMy5iQbqN6Xqmi4qUomaCrAxq3W8veBDys7o6QVIk8jXmliGviJSI1IvIbEVkqIm+LyCoReUBEhsZxrk9Evi8i74rIChF5VURO7w+7FSUuGrcBBsoPgcIyq6yrLWIxQh1tTAU3Pa+/ApXAmcaYk4CpwDTgFREp6ePcnwFXAGcYY04Afgf8W0RqMmmwosRNIN5VMRIKiqzfu9rDM+zV80oJt7uNc4wxzQDGmHrgbuBo4MJoJ4jIeGAWcJcxZpd97kPAeuCOjFusKPEQJl7F1u+hnldohr1OD0qKAhc/u9oY0xFRZq/YxuAY510CCLAgovwl4FoRKTfGNKXJRs+xq7GdbU1+3t/V8ycoLypg+EDrAero8rN5X0vU8w+tKKa00Lot9jS1s7+107GeLy+P0VWlwdcbdjfjN8axbmVpIYPLCgFobOtkZ2N71M8fU1VGfp5Adxc7PlhNe1e4V9K1fzPsfi/q+RkntI15+ZDvs37f+pa13HN7I5x3BzQExGsEbF9u/d7VHr4BR55m2KeCa+LlIFwA4wADLIpxajXgBzZFlG/Aas9xwOvpsNFrLP1gLx//1Wv4DfDfhcHyi048lAc/fTIA2w+0cc5PF0a5Ajzy+VM5/eghAMx7eT2/Xrjesd6oyhJe/ubZwdcf+/l/aWhzfgjnnH8M19UeCcDCtbv48p/fivr5b986jYpSH/zpcoavj/z/BKMBlkU9PTsoKOpJf6gYFeF5BVIlNGCfKm56XmGISD4wA/itMWZtjKpDgBZjen3jDfaxKsr1Z2F1Nxk6dCh1dXWpGZyF1G3uxG+gON8wqKgnItDRsCvY3j2tfg4plajXWL3ibbrqrYdq3/aOqHXLaQ/7G1YV+inNc667Y/N66uqsScrrd3XF/PxXXvkvpT7hw1uWUQRsYRjdWPb4jfUzsFAo80W9RD8QsN9Pnr8TMd20lI6iccBRjNr8JLzxW9qKh1MCrNi0j9Et7QwElr7+Kse3NFMM/O/1JQzb+QFHAB9s3MAG+2/Z1NSUk/dmJsga8QJuAbqAryd5fvQnAjDGzAPmAYwfP97U1tYm+THZy4ZXNsDKVZw+wsdD158Xtd5lF8R3vUT+RPHWrQW+Gk/F1/OhA0be9AqUDwPgR8+9yy/r3ucbU8bzpSlHxW9cP1CEHet4/jB47eeUtG0H4ITTzoemRdD4HhOqj4d1PmiH/zf5NFheDxvg8FEjONz+A9bV1ZGL92YmyArxEpFrgE8AtXHEq3YDpSKSH+F9DbCPezJhoxeYftJhTDh8MGuWv+m2KaljQjLRba6rPZJjZRvTTh/rklFxUHszrHgCGu3wbdhoY2SqhCappoLbo42IyGeAG4GzjTE74zhlOZbdoyLKx2J5bqvTa6F3qCovonrkIIaWuv61pk4wpaCnLQOLfQwsEop9+VFOygKKBsD5d1q/+8qgtCok5tUevgGHTg9KCVfvchG5CpgDnGuM2W6XfcSOTwXqDBeRUDv/gRXUr4243BTgBWNMY2atVvoFB8/LMxx3MUy7A6Y/ACKQb420WlOEnPZtVM8rGVzrNorIp4HfYMW6zhUJhqzOALbZdU7DGnmcB1wHYIxZIyLzgG+JyDPGmN0iMgM4Eriqf1uRXTz99lZee38Po0x3L2X3HKHrXtksXLuLny5pY13+emaecYRLhsWBCEz+cs9rR89Ll8RJFTdjXj8DirESUyOZax+bgAPYYhbCV4DvYWXjdwKNwDRjTLYPomeUJRv38pfXN3HVsYVum5I6oROYbXY2tPHO7m7GbfOYcx0W8wrxKHV6UEq4medVGUedt7GmEEWWdwLftX8Um/ZOK3biy4GQl1O3sciOdUUmrmY9YZ6XLkaYLnLhNldsAg+1Lz9m1kj2Y0zIon094lWYb92u7V0eC3CHel5hGfaB6UEqXsmg4pVDBB5qz3tewdE3seJHNkU+r4qXg+elqRIp4/XbXAkh8FAXeP1bdQjWAxTZDWvv9NjDHvC8OgNzSkXX80oDXr/NlRCC3cYo03Q8Q5Q0iaKCQMzLo55Xhy1eAVFWzyslsiLDXkkPh1eVsa+5kzKf05x3DxHF86oqK+RDw/KZdESfYz3ZRcDz6mi2joFYVyB9UT2vpFDxyiHuvOREAO9P7I3ieY0ZUsbXTi6mtvZYF4xKgYDn1WmLl6jnlQ6026hkH0HPK0duz2C3MeB52aKl04NSIkfuDgWgs9uPibIgoKcIPMwRnle337C3zc/mvdEXU8xKgt1GjXmlExWvHOLMHy/giG8/y55Wj/8nd8iuB2sV1hvqWrnogZddMCoFonUbdbQxJVS8coj2Lj/G5MBoY9RUCa+ONkYJ2Ov0oJRQ8cohOgJJqh5ciCGMKAH7QjvPq8Nr3eOoqRK2iKnnlRQqXjlET56Xy4YkSuN2a5PWAFEC9vl5Qr5Ys4c6u70kXhGeV2S3UT2vpPDaba5Eodtvgg+0p6Y2NmyDe4+Hxz/fUxYlYA89wuypydmRMS8N2KcFFa8cIdBlLCrII2RttOxn/yYrQL9vY09ZlJgXhIqXh+JeAc8rdEUJ0IB9iqh45QgBT6TIaxMbu9qsY6j3EWMV1cCKGd4Sr+Lw18GAvXpeqaAZ9jlCsS+fH19Wbe2h1PS+2+bET5e9AW2o9xHD85p5YhHV1SdRVeahBRcDnleAYMxLpwelgopXjlDsy+cTp1h7ktTVeUm8Ap5XyIa1MTyv46ryg5vieoZenldkzMtDXmQW4bE+hpJzBDwvv5PnlSO3Z6TnpTGvtKCeV46wq7Gd+Su2MXxgMUV9V88eAp5X6AMcY7Rx0ZZOXpu/mk9NGs3hVWX9YGAayMuHPB/4O63XOjE7LeTIvzZl094Wbn1qJb/0UpcRnAP2UaYHAfxvWxe/XrieTZ6b3xjSdQwuiaOeVyqoeOUI3h1tjNVtdEqVsEcbOz0WJwrtOmqeV1rw2J2uRCOQOlDktblBCQbsPZnnBRGel4pXOlDxyhFCk1Q9RcxUid5tCYhXR7fHHvhQz0tXlUgLHrvTlWi0e1a8Ap5XiCcVT5Kq57qNDjEv9bxSwmN3uhKNwI46hZ4Tr0DMK6Tb6O+9Z2OAAs92Gx1iXup5pYTH7nQlGl1+g0jPmleewTFVIrrnVVEojBxcQrHXls4I9bw0VSItaJ5XjvDJSaO58pRR+A28vGih2+bEj6PnFX208WNHFXLvzNrM25VuHD0vnR6UCipeOYSIeGs5HOhjYrbHvKtYxBxt9FgXOEvIobtD8SQBzwvT8xDH8Lw8S5jnpUmq6cBV8RKRQ0XkORHx0LKY2clDL6/ngvtf5u9LNrttSmIEPC/oeYiD4tW7Y7BwSycnzX2BHzyzqh+MSyMa80o7rnUbReQS4F6gM4lz64BhQOTW0PcYY/4vdeu8x9b9baze1sCB1k6Gum1MIgQ9L6yHON8XM2BvDBxo7aSpvavXe1mNjjamHTdjXjcDU4HvAEclcf6FxpiNabXIw4RND/LSsxDqeQWC9rm2kipEiXnZj5/fY0KcJbgpXqcZY7o8tWRxFtOTpJoP7X1UziZCPa+ABxJzelBgJVUvKTTOGfahMwj8/txZAqifcO2vZYzRfzdppGduo8cegDDPKzLm5TA9yH7ucyLDHrTrmAIeu9PDuEFEXheRd0VkkYhc47ZBbtLh9VUloEe04piY3dHtNfFyiHmF/q5B+4Txap7XfmAd8C2gDbgEeEREjjfG3OR0gojMAmYBDB06lLq6un4ytX/YtsPyYNasXsmRJW2ead/k1kYCq9G/+srLdBRVMmLLGo4G6rdt572IdnS1twHCjt17PdNGgJGb64OB3c3123jftv0MA/nAokV1+POLaWpq8lS73MST4mWMuTii6DERmQJ8XUQeMMZscjhnHjAPYPz48aa2tjbzhvYjm4o2cuL2Ri6cPIb61UvxTPte7fGgJv+/SVAxEv63GtbBiJGjGRHRjr3PvcQ3zhvNyMEl1NaM6GdjU+D198BeJ3LU6MMZFWjXq4XQ0cGZp02G4oHU1dV557tzGU+KVxQWA9cDpwC9xCvXufrDY4K/1692z46EiRnz6t1trCzO49LaZAanXcZptBE05pUCHguQgIgUikiFw1uBbz+H0rJzHH93z7ru0JMykPPTg0J8hsCghE4RSpisvztEpEpEQjfpmwz8zaHqBPv4Vuatyj7e2XKAtzfvp63TQ//BuyJyOkzf04M6ug3z39nGcyu2Z9i4NOOUKhH6u3peCZPV4iUiY4F64KmIt84RkYtC6tUCXwT+aIx5r/8szB6+9Oc3+diDr7Cjoa3vytlCV4StvZJUe0c1Wrvguj+9yXeffCfDxqWZqJ6XjjYmi5vTg+7GyrAfbb9eZr81yRgTmPbTCuwFtoac+ibwTeDbInInUIY1TegHwN39YHpW0pNh76Fec6TnlUCqhPfyvEI6D6H5a+p5JY1r4mWM+UYcdbYDh0WUNQD32D+KTSBJ1VMrqUZ6Xr0mZufSSqrRPC+dIpQsHrrTlVgEPBFPJan28rwiA/axk1T9fg8tRhIt5hUM2KvnlSgeutOVaBhjvLlvY6+YV2TAvndbRCToXXoqy77PVAkPtSVL8NCdrkSjy2/wG8jPEwryPfSVJuF5QY9Aeyru5bQYIWjAPgU8dKcr0fDuno3RYl7Rdw+CnkGJdi/t3Ri2GKEG7NNBLmXYH7QU+/J54etn0tXtoRgQJDXaCPDiDWfiy8+jtNBDI6uaKpF2VLxygPw8YdzwAW6bkTh95nk5i9Og0kLH8qwm2qoSuoNQ0nisn6HkFNEy7HNxelB+tNFG3UEoWXLo7jh4qd/fyg2PLuPef69125TESCLDHmDuP1dyxa9fY832xgwal2byC3rao4sRpgUVrxxgT1M7T7xVz3/e3eG2KYnRS7z6TlIFWLm1gcUb9rKvJXL/lSwnEPcKTQHRmFfSqHjlAB2h69d7iVRTJTyXZW93HTXDPi2oeOUA7TmXKhHb8wqmSnhpBQ3o8bxEA/bpIK13u4j8Ip3XU+IjkF3vqXmN4OB5RQbsDwbPS7uNyZJQqoSIXN1HlQtTsEVJEk/Oa4Q4AvbO7Qm0s8Nz4hWIeen0oHSQaJ7XwzHe81iGZO7Qnisxrzj2bYSe7d0863mJBuzTQaLitZre3lUZcAzwKeDBdBilJEZFiY+TRw9i7JAyt01JjKieV+zpQSeNHERjWxejKksyaFwGCHpemiqRDhIVr68aYz5wKF8lIvOBvwILUjdLSYQpxwxjyjHD3DYjcbrtVIeCYkvI4pwedOWk0Vw5aXQ/GJhmTroSjIHDPtRTpp5X0iQkXsaY/8R4r1VEjkndJOWgIeB5+Uqt3+NYw97TTPic9ROKjjYmTaIB+zOdioHBwMVYG8Aq/UxwtDE/DxFx2ZoECMS8CsuhdW/ceV4NbZ3saeqgvKiAoQOKHOt4BvW8kibRbmMdzoF5AbYAV6VqkJI4D760jgdeWsfsc49m9rnj3DYnfgKeV2GpdQxm2NsiFsXzevT1zdzx7Gpmnj6W737kuAwbmWF0tDFpEhWv94GZEWXdwE7gfWPU93UDz482FtoDDXEG7D072uhE0PPSDPtESVS8HjTGLMyIJUrSeD7D3md7XnGmShTmB8QrB/5XBqcH5UBb+pmE7nZjzH2ZMkRJnqB4+bwmXpGeV5zTg3weTVJ1QgP2SeOxu11xwpN7NkJIzCtCvPqcHmTPbcwF8dKAfdKoeOUA3l3D3va8IruNcU4Pygnx0oB90njsblec8OSGsxDieZVbxzhTJQoLcinmpZ5Xsuga9jnAjNPGMu244ZwwosJtUxIjGPOKTJWIPdp44ogK/jTzVCrLPLiWfSQ6PShpVLxygA8fWQVUuW1G4iQZ8xpUWshpRw3JsHH9hHpeSeOxfoaSM/j9IXMb7QnWcS5GmFPoaGPSqOeVAzy+dAtN7V1MP+kwBnulK9VtdxkLinvnOgUe5CgbcOxr7uCXC9+nrLCAr517dIYNzTDqeSWNel45wC/q1vG9p1eyp7m978rZQqDLWFDUO8s8cIyy9VlLZzfzFq3nr29syrCR/YCONiaNq+IlIoeKyHMiogsZpoAnpwd1hXpeEUHrvqYH5VKqhG7AkTSuiZeIXAK8BhyZ5PmzRWSViCwXkTdF5OL0WugdPDM9aNXT8NNjYON/IzyvKN3Gvtaw99oGHE5otzFp3LzbbwamAq8keqKI3Ax8F/ioMaYamAP8XUQuSK+J3iDwEGd1nlfLXnhmNjRugzXzwz0viXiA4909KBc8Lw3YJ42bd/tpxpj3Ej1JRAYBtwC/MMa8D2CM+TfwAvCT9JroDTzRbfzP7dCyx/q9oT52zKsPz8uXL4hAl9/Q7fd4xEE9r6RxbbTRGJNsJ/98oJTey02/BPxERI4xxrwb6wL+rnY2rFzs+N4hFSWU+Kwbak9TOw1tnY71CvLzGDW4NPh6w+5mou1BUllWREWJD4DGti52N0Vfs/HwyjLy8qwFBbfub42aRV5W5GPYgCKMMYzp3ggChXtWQ55Q1rQRdqyM+hkZxRjwd0JXh5W/NexY2LoMlj7cU+fAlj5iXrE9LxGhqCCP9i4/HV1+Sgrz2by3heYO51uqrLCAUZXWd9XZ7ef9XU1RzR8xqIQBxdZ3tauxPeogSEGecNSwAcHX63Y20hVFSKvKioKLJja2dVK/vzX4XmVTJ8OAA7u3snflYooaN8COoQBsP9BGa6dzm0p8BRxSURxs05Z9LVHbNGxAMWVF1qO+t7mDA63OO43n5wmjK3v2Qdi0tznqP4eKksJgknBzexc7G6Pf0yMHl+LLT7+f5MVUiWr7uCGifEPI+zHFq6J1M2P/Pq3PD0ok9XNsnPUG2D/xcFgcdQR4PrCY6K+twykAS+L8kExTOMDeNcfA8ZfCyifgQKjn5dBtDIy8RRltBCj25fODi0+ky+8H8pn7z5W8uHqnY90p44fy+2smAbCvpYPz73s56nV/c/VEph43HIA/L97EvS+udaw3pLyIJd89N/j60w8tZkeDs9B99eyjuGHaeADe2LiXGQ/3fDnX5G/iez6oWPsYFWsfs+6jpdZ7h0S1Mhwf8d9/lfZPPMS7S0BZAp+fTrwoXoHU6saI8gb76Kg3IjILmAVwwqHFvM8ox4tXFgu+fMvzaezw0+LseFGQB1UlPQ/XzmZ/1L3fBhQKpT7rmq1dhob26F2d4aUC9lLOe1v9dEYJ65QUwMAi6/O7/IaWThN87fd3k+digqeRAvx5BRR27KekbQd0NNJWNIQlgy7jNP4BjdtYsfQ1TgT2NDSzbdW7nADs3rmDFXV1nNHZQT7w8iuv0l1QGnbtpqYm6urqOGcEFO5+j6X/W2d9ZnM7I8udl8CWln3U1dUB0NhhotYDWLd6Bb6dqwHYs7Uzat0yX2fwmgCVBZ34otTdu20TdXXbAHhvb3fYNVeZk3m7+3XKsTynqmJDQb713e1v89MepTdZlA+Diu3v2xh2tUS/pwYVCUUF1mc2dxiaOp3r5gkMLe25p3e1+InWKy/zCeWF1jXbuw3726J//tASCfYm4J2o9RLGGOPqD9ZekCaB+vOw+mdVEeVT7fLr+rrGuHHjTC6zYMECt03o4cBWY1b905i9G63XPznGmO8NNObVn1vHv3zKmHeftX7/0yesOt8fZr1ub+51uaxqWwbI9fYBS0yatMOLntdu+zgA2BNSHuiN7UHJHgYeCgM/0vO6YgQ0boU971uvw1IlInfMzuIBCMV1snhsPSrL7eOYiPKxEe8r2UjFSOu4x+rukV/UE9vqtQGHF/+3Kv1F1ouXiFSJSOiEveeAFqA2ouoUYJXpY6RRcZmBI6zj3vXWMTRVwnRbo5WB6GGMgL2iZPXdISJjgXrgqUCZMWY/8H3gSyJyhF3vXOA84CY37FQSoMIeKDmw2TpGTswOeF+SFxy4UBQnXPPLReRurCD7aPv1MvutScaYQCJKK7AX2Bp6rjHmLhFpA54RkS6s7dc+boyZ3y/GK8lTMSL8dUFReKpEHwmqihLAzSTVb8RRZztR0p2MtZOR7mbkNQIxrwBhnleXBuuVuMnqbqOSgwyMFK+ino02jHpeSvyoeCn9S9kQa4QxQFiGvXpeSvyoeCn9i0h41zEsz8sf19QgRQEVL8UNQoP2kROz1fNS4kTFS+l/KkLmlUZm2GvMS4kTFS+l/xkY4XmFZtgHPS/Nrldio+Kl9D9hMa/IJNXA1CD1vJTYqHgp/U9YzCtyelBIhr2ixEDvEKX/CYt5RSxG2MfOQYoSQMVL6X8GRnpeGrBXEkfFS+l/isqheJD1u6ZKKEmi4qW4wwmXwtBjYdDo8B101PNS4kTHoxV3+Mi9Pb+HrqAa9Lz0/6oSG71DFPcJ3bdRPS8lTlS8FPcJBOyNjjYq8aPipbiPLkaoJIHGvKLQ0NDAzp076eyMsnFjFlNRUcHq1avdNiMxzvubdTxQZP1eUAwObUimbT6fj2HDhjFw4MB0WKpkCSpeDjQ0NLBjxw5GjBhBSUkJ4rG11BsbGxkwIN59ubOEre2AgcojrYW/C8thyNG9qiXaNmMMra2t1NfXA6iA5RDabXRg586djBgxgtLSUs8Jl/eJvvNyMogIpaWljBgxgp07d6b12oq7qHg50NnZSUlJidtmHFwE/kkYE/46TZSUlHgyBKBER8UrCupx9TcB8fKHv07X1fX7zDlUvJQsI73dRiV3UfFSsoPIbmOaPS8l91DxUrIE4Vf/9xjHTTgNGXEyD//lMbcNUrIcFa+DgJ07d1JTU0NlZSUiQk1NDQ899FBC15gzZw4TJ04MK1u2bBm33XYb+/fvDytvb29n7Nix3H///Ql9xrVXX86zjz9iv1LPS4mNitdBwLBhw1i2bBnTp08HLNGZOXNmwtcYPXp0WNmyZcuYO3duL/EqKChg9OjRVFVVxf8BwW5jsCAh+5SDD01SVeLixhtv5MYbb4yrbn5+PgsXLkzug4KpEsmdrhw8qOd1EBPoSo4ZM4b58+dz9tlnM3LkSKZOncqWLVuC9b70pS8xevRoRISNGzcC8J3vfIdbb70VgAsvvJCamhouv/xy9uzZQ01NDeXl5dTW1gav0dbWxs0338yECRM4+eSTqa6u5tprrw3x2gJq5UdR4kHF6yAm0JXct28fr732Gi+99BLvvvsumzZt4pvf/Gaw3oMPPsjtt98edu4dd9wRLHv22WdZtmwZjz32GFVVVSxbtqxXfGz//v38/ve/58knn+TNN9/k9ddfZ+/evVx99dVWhV5Jqplps5I7aLcxAcbc/K+o7915yYl86lQrJvTnxZv49j/eiVp3410XBX//yM9eZkV9g2O9T06SWVCVAAAQ3ElEQVQaxQ8vrQbgnS0HOHFkRTJm90ljYyOzZ88GoLy8nKlTp/LEE0+k9TOGDBnCq6++yqhR1uYbxcXFzJgxgwsuuIAdO3YwPChWmiqhxIdr4iUiw4B7gcC/6HeA2caYLdHPCp67Edjv8NZNxpgX02bkQcKQIUOorKwMvq6srGTHjh1p/YyCggLWrl3L9ddfT319PQUFBTQ1NQGwfv16hh852K6p4qXEhyviJSKFwL+BtcDxWHfs74AFIvIhY0xTX9cwxtRk1srehHpMsfjUqaODXlhfPPOVM+KqlymvC6C0tDTsdV5eHn5/emNP8+fP56KLLuKee+7ha1/7GiJCXV0dU6ZMob293SFJVVFi41bM67NANTDHGNNljOkG5gBHANe5ZJOSQf74xz9SXl7O7NmzY88zVPFS4sQt8boM2GSMWR8oMMZsB1bZ7ykewOfzAdaaWQDPP/88e/fudazb3t5OXsSmGtu3bw95FRA0FS8lPtwSr2pgg0P5BuDEeC4gIj8WkSUislZEXhCR6Wm1UOmTsWPHArBlyxYaGxu55JJLaGxsdKx70UUXceDAgWBmf2NjI/fdd19PBZ3bqCSIWwH7IcBSh/IGoFRESowxrTHO3wm8BdyClRg0C3hKRL5ijPm50wkiMsuux9ChQ6mrq4t68YqKiqgPoRfo7u4Os3/Xrl1cfPHFbN68GYDq6mpmzZrF3/72N9555x2am5uprq7m6aef5q677uIf//hHsN4999zDo48+ynPPPQfA+eefz5e//GU++9nPcuKJJ3L11Vdz1VVXUVJSwuzZszHGUF1dzfr164PXePLJJ7n88svZsGEDc+fO5d577+WQQw6htraWxYsXM2PGDL7y+U9Skg8P/O5RAG65825eXroiXOAc2pYIbW1tMb/3bKCpqSnrbcwWxLgQYxCRDuB5Y8xHI8r/BHwKKO1DvJyu+S/gDGCYMaYtVt3x48ebNWvWRH1/9erVHHvssYl8fFbhyWWgd6+DjkYoLIOOZigfDgMP61UtlbZ54Xutq6sLS+7NNURkqTFmYt81+8atbuNuwOkOHAC0JCpcNovt849PxTDFJYIhL+02KvHhlngtB8Y4lI/FyveKioiUiEi5w1v2nlnonlmeRDPslcRwS7yeAA4XkTGBAhEZDhwLPB5aUUSGi0ionVcAP3W45gSgHWvEUvEckaONql5KbNwSr4exPKwfiUiBLU53YY02/jJQSUROA7YCD0ac/0kROSWk3hXAxcCP40lwVbKQoHbpxGwlPlwZbTTGdIjIVKzpQauw/t2uAM6OEJ8m4ACwLaRsPnA38AsR8QGDgH3AtcaYef1hv5IJMrt7kJJ7uDa30RizA2tkMVadt4HKiLIdwPftHyVXEO02KomhS+IoWYLObVQSQ8VLyS40VUKJExUvJTuQiJVUVbuUPlDxUrIUVS8lNipeSpYQKVYqXkpsVLyU7CAyNUK1S+kDFS8ljPvuu48nn3yyXz+zs7OTvzz2JNM+eT0Tzv8Ux0+5nImnn8Mf/vAHnBYOWLp0KWeddRYnnHAC48eP56abbqKtLeZcfCUHUfFSwnBDvJYuXcqnv/BVrrliOkuf+zMrFzzGzTfO5nOf+1xwe7UA69atY8qUKVx66aWsWLGCxYsX8/zzz3PNNdf0q82K+6h4KVnB5FMn8smLzw++vvzSj3H66adz//33h3lfP/zhD6msrOSrX/0qAIMGDeLWW2/lr3/9K2+88Ua/2624h4pXjvP0009TU1ODiHDLLbcwZ84cJkyYwMiRI/nOd74TrLd582ZqamrYunVr8JyamhpefDHzmzGdeuqpLHjm7xGlwmGHHUZzczOdnZ0AdHV18eyzz3LWWWeFrYN/9tlnA/D444+jHDzovo05zvTp05k+fToiwh//+Ecef/xxfvSjH/HCCy9w3nnncdZZZzFt2jRGjRrFsmXLGDNmDLW1tTz88MN9XnvJkiXMnDmzz3oTJ04MLv/shIjg8xVaa4KEsHbtWj784Q9TWFgIWFukNTc3B5efDlBVVcWAAQNYvnx5n7YouYOKV7zclrmtxxLitgNJn1pTU8OECRMAmDZtGuXl5dTV1TFt2rSkrjdx4kSWLVuWtD1hRIw2vr7kTZYvX86CBQuCZbt37wZwXEl14MCB7NmzJz22KJ5Au40HEePGjQt7PXjw4LRvLpsOmppb+Px1X+EHP/gBZ555ZlznuLGcueIu6nnFSwoeT7bgtLlsd3d3lNr9jeV5dXR0ctkXbmLaOWfzrW99K6zGkCFDABw34GhsbKSqqirzZipZg4qXkjTpinkBIEJHRyeXzryR444+gp/+6I5eVY444gjKysrYuHFjWPmePXtobGykuro6EfMVj6PipYTh8/mCXbAPPviA+vp6Jk+e7Fg3nTGvjo4OLvvCTRw9djT3zr0pGAP74he/yG233cahhx5KQUEBF1xwAQsXLsQYExxxDMTFLr300rTYongDjXkpYYwdO5YtW7YAMG/evL49pjTQ0dHB5Z+ZyfpN9UyoPpZHHv8Xj/z5UR555BEWLVpEe3vPMOS3v/1t9uzZw89/bm3PeeDAAW6//XauvPJKJk2alHFblSzCGHPQ/YwbN87EYtWqVTHfz3YaGhqCvy9atMicdNJJBjDDhw83119/vdm3b5856aSTjM/nM4MHDzbnnHNOsP4rr7xijjnmGHPCCSeYU0891axZsybj9j711FMGawlVx58NGzaEte2NN94wZ555pjnuuOPM0UcfbW644QbT2tra5+d44XtdsGCB2yZkFGCJSdNzrN3GHOeMM85w7NpF6+5NnjyZ1atXZ9qsMKZPn45p2gUHNvcUDhlnbUDrwMSJE1m4cGE/WadkK9ptVLID3XBDSRAVLyVLUTFTYqPipWQJket5qXgpsVHxUrIDFSslQVS8lCxBl4FWEkPFS8lOVLuUPlDxioLRib79S69uY3rVS7/P3EPFywGfz0dra6vbZhxkZFa8Wltb8fl8ab2m4i4qXg4MGzaM+vp6Wlpa9D+2W6RJu4wxtLS0UF9fz7Bhw9JzUSUr0Ax7BwYOHAjA1q1bg0sQe4m2tjaKi4vdNiMxutqgaWfP630+yMvvVS2Ztvl8PoYPHx78XpXcQMUrCgMHDvTszV5XV8eHPvQht81IjI3/hcc/0fP6G+9D2ZBe1TzZNiUjaLdRyQ4kwssSvTWV2Lh2h4jIMBH5k4issX8eE5GRcZ7rE5Hvi8i7IrJCRF4VkdMzbbOSQfIiOgEOXUZFCcUV8RKRQuDfQCFwPHAc0AwsEJHyOC7xM+AK4AxjzAnA74B/i0hNhkxWMk1exK0Y6YkpSgRueV6fBaqBOcaYLmNMNzAHOAK4LtaJIjIemAXcZYzZBWCMeQhYD/ReO1jxBup5KQnilnhdBmwyxqwPFBhjtgOr7PdicQnWQPqCiPKXgGlxem5KttEr5qXipcTGLfGqBjY4lG8ATozjXD+wyeHcAqwuqOI11PNSEsStVIkhwFKH8gagVERKjDHRUtyHAC12VzPyXADH/a9EZBZWdxOgXURWJGizlxgC7HbbiJSYG/XW9H7bYpPr7RufrgtlW55XKnnVMc81xswD5gGIyBJjzMQUPiuryeX25XLb4OBoX7qu5Va3cTfQe892q6wlhtcVOLdUpFdQJHA93fNdUQ4C3BKv5cAYh/KxwDtxnJsHjHI4twvo390jFEVxBbfE6wngcBEZEygQkeHAscDjoRVFZLhIWLr1P7C2xKqNuOYU4AVjTO+94HszL3GTPUUuty+X2wbavrgRN1ZNsJNUl2B5SZ/GGj38LXA68CFjTJNd7zRgETDPGHNdyPm/whKr04wxu0VkBvAg8GFjTHq2cFYUJatxxfMyxnQAU4FurNyu1cBA4OyAcNk0AQeAbRGX+Arwd+AVe9TwC8A0FS5FOXhwxfNSFCdE5FDg98B5xpicWgg6l9vmFjkzdT/XJ3qn2L6NIrLM4efcTNsdLyJyCfAacGSS588WkVUislxE3hSRi9NrYfKk0jYRqbPbFfndXZ1+SxNHRGpE5DcislRE3rZtfUBEhsZxbmrPnTHG8z9YE7zfxupKFgD5wB+A94DyOM7/FbAWGGq/ngm0AjVuty1N7dvodhvisHExcDTwsHVbJnTuzVgpNEfar6cCncAFbrcrDW2rA8a43YYY9r2LNchWZr8eYZetBUr6ODel5871xqfpD/gFrBHII0LKDsGKqX2jj3PHYw0YzIgoXwn8y+22pdo+u+5Gt9sQh40F9jGhBxwYhLUiye0R5f8CVrrdrlTaZp/jBfE6KqLs8/b9elmM81J+7nKl25jrE71TaZ8nMMZ0JXnq+UApzt/fcSJyTEqGpYEU2uYFqo0x6yLKttrHwTHOS/m5yxXxyvWJ3qm0DwAR+bGILBGRtSLygohMT6uF7lFtHyP/Phsi3vcyN4jI63ZsaJGIXOO2QQGMlTkQyTgsz2tRjFNTfu5yRbyGAE7JqcGJ3n2cm/BE734mlfYB7ATeAk7DWvzxKeApEflyWq10h8BC95F/n2z6/lJhP7AOK6/xeOAB4Jci8hNXrYqCPW1vBvBbY8zaGFVTfu6ybWJ2usnYRO8sIS4bjTGTIooeFJELgTtF5CFjTFv6TXMdL3x/fWKMiRw1fUxEpgBfF5EHjDGRnovb3II1Te/rSZ4f9/eWK55Xrk/0TqV90Vhsn398KoZlAYHlYyL/Ptn0/aWbxVjP7iluGxKK3Z39BNYob1Mf1VN+7nJFvHJ9onfS7RORkijBz4C77vVV/5bbxzER5WMj3vccIlIoIhUOb2XddycinwFuxJols7Ov+qThucsV8XJ7onemSaV9VwA/dbjmBKAda8TSM4hIlT03NsBzQAvO398qY8y7/WVbqji0bTLwN4eqE+zjW5m3qm9E5CqsPSjOtUfBEZGP2AuABuqk/7lzO08kTbkmhVhK/ihWHC8PaypGWBInVsC6G/hlxPm/AtYAQ+zXM8i+JNWk2gd8DisIekpI2RVYIz2391cbEmjrw0TJhcL6r9wGzI8ovxnYhZ0HB5xLFiWpJts2+8HuAi6KKGsG/s/t9tj2fNp+Vm4Crgr5+TVwm4lyX9rlKT13ORGwN8Z0iMhU4F4sT8IAK0hsovf3sCZ6d2KNXGXNRO8U2zcfuBv4hYj4sJI69wHXGmt12axARO7Gyowfbb8O/O0nmZ7h+FZgLz15RAAYY+4SkTbgGRHpwnpQPm6Mmd8vxvdBCm17E/gm8G0RuRMoAzqAH2B9p9nAz4BinO2Zax8z8tzpxGxFUTxJrsS8FEU5yFDxUhTFk6h4KYriSVS8FEXxJCpeiqJ4EhUvRVE8iYqXoiieRMVLURRPouKlKIonUfFSFMWTqHgpiuJJVLyUrMbeq7JBRPwi8qJd9qCI7BORDSIy020bFXfQidlK1iMiH8da1+oLxpiHRORw4J/AZNP3ip1KjqLipXgCEXkCa52uGuB3wJ3GmBfctUpxExUvxROIyCFYa5l1A88YY7Jm+y/FHTTmpXgCYy0vPBdry6zIjUqVgxD1vBRPYK9/XgeUYK1IerwxZnfMk5ScRj0vxSt8DWvLr4uxlh2+311zFLdRz0vJekTkSOAxrNHFVhH5ItbmDR81xjzjrnWKW6jnpWQ1InIH8F/gEKzdZQCut49/EpHHXDFMcR31vBRF8STqeSmK4klUvBRF8SQqXoqieBIVL0VRPImKl6IonkTFS1EUT6LipSiKJ1HxUhTFk6h4KYriSVS8FEXxJP8fSt8c/WgCo8oAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "linear_convection(85)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oops. This doesn't look anything like our original hat function. Something has gone awry. It's the same code that we ran each time, so it's not a bug!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### What happened?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To answer that question, we have to think a little bit about what we're actually implementing in the code when we solve the linear convection equation with the forward-time/backward-space method. \n", "\n", "In each iteration of the time loop, we use the existing data about the solution at time $n$ to compute the solution in the subsequent time step, $n+1$. In the first few cases, the increase in the number of grid points returned more accurate results. There was less discretization error and the translating wave looked more like a square wave than it did in our first example. \n", "\n", "Each iteration of the time loop advances the solution by a time-step of length $\\Delta t$, which had the value 0.025 in the examples above. During this iteration, we evaluate the solution $u$ at each of the $x_i$ points on the grid. But in the last plot, something has clearly gone wrong. \n", "\n", "What has happened is that over the time period $\\Delta t$, the wave is travelling a distance which is greater than `dx`, and we say that the solution becomes *unstable* in this situation (this statement can be proven formally, see below). The length `dx` of grid spacing is inversely proportional to the number of total points `nx`: we asked for more grid points, so `dx` got smaller. Once `dx` got smaller than the $c\\Delta t$—the distance travelled by the numerical solution in one time step—it's no longer possible for the numerical scheme to solve the equation correctly!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![CFLcondition](figures/CFLcondition.png)\n", "#### Graphical interpretation of the CFL condition." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the illustration above. The green triangle represents the _domain of dependence_ of the numerical scheme. Indeed, for each time step, the variable $u_i^{n+1}$ only depends on the values $u_i^{n}$ and $u_{i-1}^{n}$. \n", "\n", "When the distance $c\\Delta t$ is smaller than $\\Delta x$, the characteristic line traced from the grid coordinate $i, n+1$ lands _between_ the points $i-1,n$ and $i,n$ on the grid. We then say that the _mathematical domain of dependence_ of the solution of the original PDE is contained in the _domain of dependence_ of the numerical scheme. \n", "\n", "On the contrary, if $\\Delta x$ is smaller than $c\\Delta t$, then the information about the solution needed for $u_i^{n+1}$ is not available in the _domain of dependence_ of the numerical scheme, because the characteristic line traced from the grid coordinate $i, n+1$ lands _behind_ the point $i-1,n$ on the grid. \n", "\n", "The following condition thus ensures that the domain of dependence of the differential equation is contained in the _numerical_ domain of dependence: \n", "\n", "$$\n", "\\begin{equation}\n", "\\sigma = \\frac{c \\Delta t}{\\Delta x} \\leq 1\n", "\\end{equation}\n", "$$\n", "\n", "As can be proven formally, stability of the numerical solution requires that step size `dt` is calculated with respect to the size of `dx` to satisfy the condition above. \n", "\n", "The value of $c\\Delta t/\\Delta x$ is called the **Courant-Friedrichs-Lewy number** (CFL number), often denoted by $\\sigma$. The value $\\sigma_{\\text{max}}$ that will ensure stability depends on the discretization used; for the forward-time/backward-space scheme, the condition for stability is $\\sigma<1$.\n", "\n", "In a new version of our code—written _defensively_—, we'll use the CFL number to calculate the appropriate time-step `dt` depending on the size of `dx`. \n", " " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def linear_convection_cfl(nx, L=2.0, c=1.0, sigma=0.5, nt=20):\n", " \"\"\"\n", " Solves the 1D linear convection equation\n", " with constant speed c in the domain [0, L]\n", " and plots the solution (along with the initial conditions).\n", " Here, the time-step size is calculated based on a CFL constraint.\n", "\n", " Parameters\n", " ----------\n", " nx : integer\n", " Number of grid points to discretize the domain.\n", " L : float, optional\n", " Length of the domain; default: 2.0.\n", " c : float, optional\n", " Convection speed; default: 1.0.\n", " sigma : float, optional\n", " CFL constraint; default: 0.5.\n", " dt : float, optional\n", " Time-step size; default: 0.025.\n", " nt : integer, optional\n", " Number of time steps to compute; default: 20.\n", " \"\"\"\n", " # Discretize spatial grid.\n", " dx = L / (nx - 1)\n", " x = numpy.linspace(0.0, L, num=nx)\n", " # Compute the time-step size based on the CFL constraint.\n", " dt = sigma * dx / c\n", " # Set initial conditions.\n", " u0 = numpy.ones(nx)\n", " mask = numpy.where(numpy.logical_and(x >= 0.5, x <= 1.0))\n", " u0[mask] = 2.0\n", " # Integrate the solution in time.\n", " u = u0.copy()\n", " for n in range(1, nt):\n", " u[1:] = u[1:] - c * dt / dx * (u[1:] - u[:-1])\n", " # Plot the solution along with the initial conditions.\n", " pyplot.figure(figsize=(4.0, 4.0))\n", " pyplot.xlabel('x')\n", " pyplot.ylabel('u')\n", " pyplot.grid()\n", " pyplot.plot(x, u0, label='Initial',\n", " color='C0', linestyle='--', linewidth=2)\n", " pyplot.plot(x, u, label='nt = {}'.format(nt),\n", " color='C1', linestyle='-', linewidth=2)\n", " pyplot.legend()\n", " pyplot.xlim(0.0, L)\n", " pyplot.ylim(0.0, 2.5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, it doesn't matter how many points we use for the spatial grid: the solution will always be stable!" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEbCAYAAACLNQJjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VNX5+PHPk2SysyUQZIugsogSoqDWPaigiFK31rb6tZVaam1tqVrRWv2qra2tdaultWir/Wm/3dSqteJWCVC1KCiyBRABgbCEPfs65/fHnRmSYSaZZO6dO3d43q9XXpO5c5fnZGaenHvuueeIMQallPKaNLcDUEqpntDkpZTyJE1eSilP0uSllPIkTV5KKU/S5KWU8iRNXkopT8pw46AiUgp8GzgxEIMPeAv4sTFmVxfbbgL2R3jpFmPMWzaHqpRKUuJGJ1URWQOsAq4xxtSJyBDg31g1wfHGmIZOtt1kjBmemEiVUsnKzdPG2caYOgBjTCXwADASuNDFmJRSHuHKaSNQYoxpDlu2LfDYL9HBKKW8x5WaV4TEBTAKMMDCrrYXkV+IyBIRWScib4jIdNuDVEolNbdqXh2ISDowA/i9MWZdF6tXAR8BdwJ+YCbwkojcaIz5dSfHmBlYl+zs7AnFxcW2xJ6M/H4/aWmpeSE5lcsGqV++devW7TbGDLBjX6402B8ShMjdwMXA2caY2h5s/y/gTKDIGNPY1fqjR482a9eu7XacXlFeXk5ZWZnbYTgilcsGqV8+EVlqjJlox75cT/Eici3wRWBqTxJXwGKgF3CcbYEppZKaq8lLRP4HuBk4xxhTFcP6OSKSH+GltsBjup3xKaWSl2vJS0SuBmYD5xljdgSWXRRomwquM1BE2sd4JfBghN1NAJqA1Q6GrJRKIq4kLxG5CngCeBo4T0SuDiSzi4HBgXVOx+o+MSds8y+LyEnt9nUlcAnwizhOO5VSHuPW1cbHgGysjqnh7gk81gIHgO3tXpsX2OY3IuID+gL7gOuNMXOdC1cplWxcSV7GmIIY1vkYKAhbthP4ceBHKXUYc/1qo1JK9YQmL6WUJ2nyUkp5kiYvpZQnafJSSnmSJi+llCdp8lJKeZImL6WUJ2nyUkp5kiYvpZQnafJSSnmSJi+llCdp8lJKeZImL6WUJ2nyUkp5kiYvpZQnafJSSnmSJi+llCdp8lJKeZImL6WUJ2nyUkp5kiYvpZQnafJSSnmSJi+llCdp8lJKeZImL6WUJ2nyUkp5kiYvpZQnafJSSnmSJi+llCdp8lJKeZJryUtESkXkCRFZKiIfi8hqEfmViAyIYVufiPxYRNaIyEoReVdEzkhE3Eqp5OBmzesvQAFwljFmPDAZmAK8IyI5XWz7GHAlcKYx5njgD8CbIlLqZMBKqeTh9mnjbGNMHYAxphJ4ABgJXBhtAxEZDcwE7jfG7Aps+ySwAbjP8YiVUkkhw8VjlxhjmsOWbQs89utku0sBAeaHLX8buF5E8o0xtTbF6Dm7aprYXuvn010H/wT5WRkM7J0NQHOrny376qNuP6hPNrmZ1sdiT20T+xtaIq7nS0ujuDA39Hzj7jr8xkRctyA3k355mQDUNLZQtb+a9KYapLUejB8xbRDYdkjfXNLTBICd1Y00tbZ12Ffr/i2w+5NDDyJp1k9aBmT3hsxekOb2/2blJNeSV4TEBTAKMMDCTjYtAfzA5rDlG7HKMxZ4344YvWbpZ3v5wuPv4TfAfxaElk8bN4g5V50IwI4DjZz74IIoe4Bnv34KZ4zsD8DcRRv43YINEdcbVpDDolvPCT3//K//Q3Vja8R1Hzq1ict6rYbP3iVn60cc3RY9ebY3MMKyYoBlsWwt0GcoFJ8KR54GYz8PuQUxHVd5g5s1rw5EJB2YAfzeGLOuk1X7A/XGmLaw5dWBx8Io+5+JdbrJgAEDKC8vjy/gJFS+pQW/gex0Q9+sg7WO5updofLuafBzRK5E3UfFyo9prUwHYN+O5qjr5tPU4W9YmOknN63julmmie+ZZ7jso3+HlmUAraRTQy6NZNFGGoY0/FjbFmYLwd3sbzK0+A/uz2+sn96ZQp6vYzxi/ICfNH8r6W31ZLQ1woEtsGILrPgbTW/cQ8Wx32d/v/FRy54MamtrU/Kz6QQxUar6iSYidwMXA2d3dtonIm8ApxpjeoUt/wYwF7jQGDOvs2ONHj3arF27Nv6gk8xT72zknn+u5rziDJ684Xx3g6laA3//GuyqgPRMmDgDRpwFw06B3EKQ6Ak0mp+/tobfln/KD84fzbcnHdP5ym2tsHstfPYuLP8rbP0AEDhjFpxzJ6Sl96hYTisvL6esrMztMBwjIkuNMRPt2FdS1LxE5Frgi0BZDO1Vu4FcEUkPq30Fk9keJ2L0gunjBzPhyH6sXf6hu4E018Gzl0P1Vug/Ci7/PQwqiXu33yo7mmNlO1POGNH1yukZMPA462fiDFj4S1hwP/znYcjuA2d8P+54lLtcb9EUkf8BbgbOMcZUxbDJcqy4h4UtHwG0AhX2RugdhflZlAzty4Bcl9/WRQ9aiWvQeJhZbkviAuid7aN3lpDt62atKS0dymbDlc9azxc8ANXbOt9GJT1XP+UicjUwGzjPGLMjsOyiQPtUcJ2BItI+zn9gNeqXhe1uEvCGMabG2ahVp/Z8Cu8+Zv1+4YOQmeduPO2NmQZjLoKWOnjzLrejUXFys4f9VcATwNPAeSJydSCZXQwMDqxzOlb3iTnB7Ywxa7Hatm4Xkf6B9WYARwN3JLIMyeblj7dx+wsrWL0n/FpGAr12O7Q1w/ivwLCTbN31gnW7eHBJI08uinwFNCbn3wcZ2bDi71Z7mPIsN2tejwHZWB1Tn2n3M7PdOrXAAWB72LY3An/H6o2/EvgGMMUYE9NF9FS1ZNNe/vz+ZrbV+rte2Qnr34JPXrf6WJ13t+27r6puZMXuNiq2x1G57jccTv+e9fu8W0P9y5T3uNnPq8tON8aYj7FuIQpf3gL8KPCjApoC/Qp8bv1L+vAZ6/GMWdArUi+t+GQF2rrCO6522+mzYMlTsGMFbP8YButdZV7keoO9sk/wS+1L7343hLi1NMInb1q/l3zRkUNkplsf16bWOGuWmblw7MXW7xX/jDMq5RZNXikk+KV2pea1Yb7VED5oPPQtduQQWT6bkhfAsRdZj2teiX9fyhWavFJI8Eud4ca7WhFIAsEajQOyAgVrarHhgsTwM63+XrvWRL5XUiU9TV4pJHTamJbg08a2Vlj7qvX7GCeTV7DNy4aaV7oPRk21ftdTR0/S5JVCjizMY+yg3ofc9+e4ze9Cw14oPAYGjHbsMIV5mZxQlM4pR9l0g7WeOnpaUtwepOzx00vHAST+xt72p4w9uGcxVsP75/G9E7MpKzvWnh0efS5k5EDlUjhQCX2G2LNflRBa81LxMeZgzcXBU0ZHZObCMedav6/5l7uxqG7T5JVCWtr8JHyUkKrVUF0JvQbB4BMcPVSb37C30c+WvbGNBxaT4AWG9W/at0+VEJq8UshZv5jPUT98lT0NCexhv3WJ9Vh8quMjl9Y0tnBTeQPTfrXIvp0Wn2o9bl2ive09RpNXCmlq9WNMgq82Vi61HofaMkRTp2y92hjUtxjyBlgXHPZtsm+/ynGavFJIc7CTaiLH2asMjB02ZILjh8oM9PNqtvP0WORg7MFErDxBk1cKOdjPK0EHbK6z2rwkHY6wZ8yuzqSnCelind21tNl4ihdKXi4P4qi6RZNXimjzm9AXOmG3Nm7/GEwbDBxrXblLgGBijvvm7PaGWJOTULnEvn0qx2nyShHBU8asjDTEwb5WHQRPsxJwyhh0MHnZ2O41OJC8tn8MbZGnelPJR5NXigjWRLISeWNjKHk531gfFBwxw9bklVsABUdDa6N1Gqw8QXvYp4hsXzq/uLzEmo639tPEHNSFmtd147IoKRlPYWASW9sMmQB7P7XKNCi5p0dTFq15pYhsXzpfPGkYX5wYPi+JQ2p3wf7N4Mtz9H7GcGML0zljZP/uT8LRlWAC3qpXHL1Ck5fqmWCta/AJSTsHYrdodwnP0eSVInbVNPH/3tvE66t2JOaAoVPGExNzvICFW1v42bwKPttTZ++OjxgHaT5rfK8mnYDKCzR5pYjNe+u566VV/LY8we1dCehZ395/t7fyuwUb2Gzn/Y0Avmw44njAwLaP7N23coQmrxSR8KuNO1ZYjwlu3A7e+hScbMRWwbLsWGn/vpXtNHmliGDXgaxE3BtUtxvqqiAzH/o4M159NI708woqOs561O4SnqDJK0W076TquOCXe8AYx0eSCBdMXs1tDkysWxQY5FCTlydo8koRTQlNXhXW48Cxzh8rTKiTqhOnjUWB8lStAb9LE/eqmGnyShHBGXUyE1nzKkp88spw8rQxrxDyB1pTuB3YbP/+la00eaWIVr9B5OCYV47a6V7y6pMpDO2XQ7ZTQ2cEy7RTTx2Tnd4elCK+fHIxXzppGH4DixYucO5Axhw8bXQheX3+mEwevq7MuQMUjbUm0K1aDWMudO44Km6avFKIiDg/HM6BLdBcY40+mj/A4YO5INiOp432SU9PG1X3hGpdNk0/lmxCVxwr3I1DdcnV5CUig0TkNRHRmQ/i9OSiDUx9dBF/X7LF2QO52FgPsGBrC+PveYOfvOJQzWjAGEBg9zpobXbmGMoWriUvEbkUeA84ugfblovIahFZFvZzjf2ResO2/Y1UbK/mQIPDg+m52FgPVpPbgYYWaptanTlAZh70Gw7+Vtiz3pljKFu42eZ1GzAZuAM4pgfbX2iM2WRrRB7W4fYgB/pvhrjYWA8O97APKhoL+zZatUwX+rKp2Lh52ni6MeYTF4+fUg52UnWwq0RbK+xea/1eNMa543QidG+jnWPYhws12mu7VzJzLXkZYxyq9x+eDt7b6OBbuvdTaGu25jrM6uXccToRvHXTkR72QXqbkCd4+WrjTSLyvoisEZGFInKt2wG5qTkRo0qEGuuPc+4YXTh4b6OTyUtv0PYCr/bz2g+sB24HGoFLgWdF5DhjzC2RNhCRmcBMgAEDBlBeXp6gUBNj+85GANZWrOLonEZHyjd846sMBz5rzGWjS3+/1qZGQNi5e69j76H4WzlTMkjbt4lFb82jLSPHkeNEUltbm3KfTad4MnkZYy4JW/SciEwCvi8ivzLGHHJjmjFmLjAXYPTo0aasrMz5QBNoc9Ymxu2o4cLThlNZsRRHyrfjCQCOnDiVI0sc2H8M9r72Nj84v5ih/XIoKx3i3IHWjIaqVZw5pgiGJm6CkfLycmfeuxTkyeQVxWLgBuAk4LC7q/aaU4eHfq90qp3ZxdEkggqy07isrCcXp7tp4FioWmWdOiYweanYea7NS0QyRaRPhJeCl59SYDaIJNRcD3s3QFoGFI50OxrnaaN90kv65CUihSLSfpK+04C/RVg1+O/xsByAfMXWA3y8ZT+NLQ51Idi9FjBQeAxk2DxnYjc0txnmrdjOaysdnmikSO9xTHZJnbxEZARQCbwU9tK5IjKt3XplwDeBZw7XvmPf/r8P+fycd9hZ3ejMAVzuWR/U0Arf+tOH/OjFFc4eSIfGSXqutXmJyANYPeyLA8+XBV462RgTvKmsAdgLbGu36YfArcAPReSnQB7QDPwEeCABoSelgz3sHTprdvmexqBQD3sn+3kB9BlmjdFfV2WN2Z/X39njqW5zLXkZY34Qwzo7gMFhy6qBhwI/KiDYSdWxkVSDycvl22UcHUm1vbQ0q91r6wdW2Uec5ezxVLcl9Wmjil2wJuJYJ9UkGQqnfSdVv9/hwUh0eJykpskrBRhjnJ23sX4v1GwHXy70HW7//rtBREK1S0d72YM22ic5TV4poNVv8BtITxMy0h14S4M1DxemOoskmKAdb/fSRvuk5v4nUcXN8Tkbk6SxPih4UaLJibkb2ytqN7qE0fEyk00q9bA/bGX70nnj+2fR2ubQFyxJGuuD3rrpLHzpaeRmOtwfOX+ANVZ/3S5r7P6+iZ0dXHVOk1cKSE8TRg10cIiaJGmsD+qbm8BOskXHwsZd1t9Ak1dS0dNG1TljkqaDqitC7V6r3I1DHUKTVwqo3N/ATX9dxsNvrrN/59WV0HQAcgqs2aSTwD3/XMWVv3uPtTtqnD+YJq+kpckrBeypbeKFjyr595qd9u+88kPrcXApiNOTQsZm1bZqFm/cy776BMzuM7jUetx2WN4ym9Q0eaWAZifHr69caj0OSZ5hYUJdJZzuZQ9WzSsj2xoCu36v88dTMdPklQKanOwqkZTJK9BVwqkRNNpL98Gg8dbv2z50/ngqZrZ+2kXkN3buT8Um2Lve9vsa/W0HT5eSKnklsOYFB8teqckrmXSrq0QMk7peGEcsqoccu69x9zporoU+xZBfZO++4xAsZ3PCk9fSxBxPxaS7/bye7uQ17YLsEsfmbAydMp5o737jFJzeLfE1r6VW15EkuXBxuOtu8qrg0NpVHjAG+Aowx46gVPf0yfFxYnFfRvTPs3fHSdjeBTB+aF9qGlsZVpCgWX36Dbe6imhP+6TS3eT1XWPMZxGWrxaRecBfgPnxh6W6Y9KYIiaNceC0busS6zHJkteXTi7mSycnMIGIWH+D9W9afxNNXkmhW40kxph/d/JaA1YNTKWClgarY6akHbzadjjTdq+k090G+0jDSQrQD7gEawJYlWChq43paYhd7THbl4Npg4HHQ1a+Pfu0SXVjC3tqm8nPymBAr6zEHHToROtRrzgmje6eNpYTuWFegK3A1fEGpLpvztvr+dXb65l13khmnTfKnp0maWM9wF/f38J9r1Zw3Rkj+NFFCbrfcnDg77B9GbS1QrqOaeC27r4DnwLXhS1rA6qAT40xCeg1qMI5crWxMjnbu8CFq40AeYVWw/2+TbCrAo4Yl7hjq4i6m7zmGGMWOBKJ6jHbe9j722BD4G0eerI9+7RRZnoweSX4f+Wwz1nJ69P5mrySQHcb7B9xKhDVc6Hk5bMpeW1ZDPW7rZpGkozh1V6wnAnrpBo0eqr1uOaVxB5XRaT3NqYA2+dsrPin9XjsxUnZITN0b2Oik9cx51k3aW9ZDDUOz9ituqTJKwXYOoa9MVARqFmMuTj+/Tkg4fc2hg6cD0efY/2+5l+JPbY6hCavFGDrhLPbP4YDm62BB4eeFP/+HJCZ4VKbF8CYi6zHYO1UuUav96aAGaePYMrYgRw/pE/8Owu254yZlhTTnEUybkgf/nTdKRTkJXAs+6DRU0HSYdMiaNgHOf0SH4MCNHmlhFOPLgQK7dlZ+/auJNU3N5PTj+nvzsFzC2D46bBxIax7HcZ/yZ04lJ42qnZ2r4ddayC7Dww/0+1oklewLVBPHV2lySsFPL90K398dxP76uIc0/2DJ6zHURdYI4gmqX11zfz01QoefesTdwIYM816/ORNq9+XcoUmrxTwm/L1/O/Lq9hT19TznexcDe8/Yd2Ifep37AvOAfUtbcxduIG/fLDZnQD6DIFxX4S2Jnj9DndiUO4mLxEZJCKviYgOZBiHuG8PMgbm3WrdiD1xBgwqsTE6+7nWVaK9yfdCZr51gWN91MFWlINcS14icinwHnB0D7efJSKrRWS5iHwoIpfYG6F3xH170OoXratnOQUwKflrEqHklYgJOKLpPQjO+oH1+7zZ0JqAadhUB27WvG4DJgPvdHdDEbkN+BFwsTGmBJgN/F1EptobojcEv8Q96ue1YwW8drv1+7l3WlfTkpxrPezDfe4GKDwG9nwCr2kCSzQ3k9fpxphut7iKSF/gTuA3xphPAYwxbwJvAL+0N0Rv6NFpozHw38fhiXOgZjsMOwVO/KpDEdrLly6IQKvf0OZ3scUhIxOmPWj1+1ryB/jDFNjzqXvxHGZc6+dljGnt4aYXALkcOtz028AvRWSMMWZNZzvwtzaxcdXiiK8d0SeHHJ+VBPbUNlHd2BJxvYz0NIb1yw0937i7jmhzkBTkZdEnx7p6V9PYyu7a6GM2HlmQR1qadT/htv0NUXuR52X5KOqVhTGG4W2bQCBzTwUI5NVugh0rrXj8bVZbVksDNFZDXRVs+cA6TdwfGNF74gyYch+kOTBprQNEhKyMNJpa/TS3+snJTGfL3nrqmiN/pPIyMxhWYL1XLW1+Pt1VG3XfQ/rm0Cvbeq921TRFvQiSkSYcU9QLjiqDGa/R8rcZ+LZ9hH/O52gceCL1gz9HU7+R+LN6k9+7gH75uZCWTl2zn6ra6DW0rJqNsHMAADsONNLQErlMOb4MjuiTHSrT1n31UfdZ1CubvCzrq763rpkDDZGPn54mFBccnAdh8966qP8c+uRkhjoJ1zW1UlUT/TM9tF8uvnT760le7KQabE3eGLZ8Y7vXO01efRq2MOLvU7o8UHe6fo6Icb1egZ9YDI5hHQFeDw4m+jvr4SSAJTFsnNsfLnoYxk6PMaLkke1L5yeXjKPV7wfSueefq3iroiriupNGD+Cpa62hffbVN3PBI4ui7veJayYyeexAAP5v8WYefmtdxPX652ex5EfnWU+Gncz5jffx3bbHuYR3yd3+X3K3/zfidnl0/lkZARAYB/KITtZrz9fFPtsrCPzEItaR+rsqk1O8mLyCXatrwpZXBx4j5hsRmQnMBDh+UDafMizizguyBV+6VfOpafZTH7niRUYaFOYc/G9SVeePOvdbr0wh12fts6HVUN0U/VRnYK6ERnLY2+CnJUqzTk4G9M6yjt/qN9S3mNBzv99PWuDWHiNpQBpt6Vm0ZuTRmpFHbf7R7O87lrq8EZiqdKgqjxpPsqmtraW8vJxzh0Dm7k9Y+t/1AJi6JobmRx4BQ+r3UV5eDkBNs4m6HsD6ipX4qioA2LOtJeq6eb6W0D4BsnyZ/NLcyBPmq4wz6yg1FRSxh3zqGZTZQF+fQYyfpjY/+xr9RIugMNuQkW7VgPc3+mmKck0iKx36Zgfeb2PYVR/9M9U3S8jKsI5Y12yobYm8bprAgNyDn+ld9X6inZXn+YT8TGufTW2G/Y3Rjz8gR0JnE7Ai6nrdZoxx9QdrLkjTjfXnYp2fFYYtnxxY/q2u9jFq1CiTyubPn+92CI5J5bIZk/rlA5YYm3KHFzup7g48hp99BZ/vSWAsSimXeDF5LQ88Dg9bPiLsdaVUCkv65CUihSLSfuyT14B6oCxs1UnAatPFlUalVGpI6uQlIiOASuCl4DJjzH7gx8C3ReSowHrnAecDt7gRp1Iq8Vy72igiD2A1shcHni8LvHSyMSbYEaUB2Atsa7+tMeZ+EWkEXhGRVqzp175gjJmXkOCVUq5zs5PqD2JYZwdRujsZayYjnc1IqcNUUp82KqVUNJq8lFKepMlLKeVJmryUUp6kyUsp5UmavJRSnqTJSynlSZq8lFKepMlLKeVJmryUUp6kyUsp5UmavJRSnqTJSynlSZq8lFKepMlLKeVJmryUUp7kxXkbE6K6upqqqipaWqJM3JjE+vTpQ0VFhdthOKInZfP5fBQVFdG7d2+HolJu0OQVQXV1NTt37mTIkCHk5OQgEn2S0mRUU1NDr16xzsvtLd0tmzGGhoYGKisrATSBpRA9bYygqqqKIUOGkJub67nEpToSEXJzcxkyZAhVVVVuh6NspMkrgpaWFnJyctwOQ9koJyfHk00AKjpNXlFojSu16PuZejR5KaU8SZOXUsqTNHmppPH4448zduxYRISnn37a7XBUktPkdRioqqqitLSUgoICRITS0lKefPLJbu1j9uzZTJw4scOyZcuWcffdd7N///4Oy5uamhgxYgSPPvpot45x/fXX8+qrr3ZrG3X40uR1GCgqKmLZsmVMnz4dsJLOdddd1+19FBcXd1i2bNky7rnnnkOSV0ZGBsXFxRQWFsYXuFKd0E6qKiY333wzN998c0zrpqens2DBAocjUoc7rXkdxoKnksOHD2fevHmcc845DB06lMmTJ7N169bQet/+9rcpLi5GRNi0aRMAd9xxB3fddRcAF154IaWlpVxxxRXs2bOH0tJS8vPzKSsrC+2jsbGR2267jQkTJnDiiSdSUlLC9ddff0itTalYafI6jAVPJfft28d7773H22+/zZo1a9i8eTO33npraL05c+Zw7733dtj2vvvuCy179dVXWbZsGc899xyFhYUsW7bskPax/fv389RTT/Hiiy/y4Ycf8v7777N3716uueYa5wuqUpKeNnbD8Nv+FfW1n146jq+cYrUJ/d/izfzwHyuirrvp/mmh3y96bBErK6sjrvflk4fxs8tKAFix9QDjhvbpSdhdqqmpYdasWQDk5+czefJkXnjhBVuP0b9/f959912GDRsGQHZ2NjNmzGDq1Kns3LmTgQMH2no8lfpcS14iUgQ8DAT/Ra8AZhljtkbfKrTtJiDS+cYtxpi3bAvyMNG/f38KCgpCzwsKCti5c6etx8jIyGDdunXccMMNVFZWkpGRQW1tLQAbNmzQ5KW6zZXkJSKZwJvAOuA4wAB/AOaLyAnGmNqu9mGMKXU2ykO1rzF15iunFIdqYV155cYzY1rPqVoXQG5ubofnaWlp+P1+W48xb948pk2bxkMPPcT3vvc9RITy8nImTZpEU1OTrcdShwe32ry+CpQAs40xrcaYNmA2cBTwLZdiUg565plnyM/PZ9asWXqfobKFW8nrcmCzMWZDcIExZgewOvCa8gCfzwdYY2YBvP766+zduzfiuk1NTaSldfy47dixw9kAVUpzK3mVABsjLN8IjItlByLyCxFZIiLrROQNEZlua4SqSyNGjABg69at1NTUcOmll1JTUxNx3WnTpnHgwIFQz/6amhoeeeSRhMWqUo9bDfb9gaURllcDuSKSY4xp6GT7KuAj4E7AD8wEXhKRG40xv460gYjMDKzHgAEDKC8vj7rzPn36RP0SekFbW1uH+Hft2sUll1zCli1bACgpKWHmzJn87W9/Y8WKFdTV1VFSUsLLL7/M/fffzz/+8Y/Qeg899BB//etfee211wC44IIL+M53vsNXv/pVxo0bxzXXXMPVV19NTk4Os2bNwhhDSUkJGzZsCO3jxRdf5IorrmDjxo3cc889PPzwwxxxxBGUlZWxePFiZsyYwfXXX09mZiY6AJKYAAANeklEQVSPP/44AHfeeSeLFi06JMGFl607GhsbO33fk0FtbW3Sx5gsJFjlT+hBRZqB140xF4ct/xPwFSC3i+QVaZ//As4EiowxjZ2tO3r0aLN27dqor1dUVHDsscd25/BJRYeBjswL72t5eXmHzr2pRkSWGmMmdr1m19w6bdwNRPoE9gLqu5u4AhYHtj8unsCUUt7gVvJaDgyPsHwEVn+vqEQkR0TyI7zUFnhMjy80pZQXuJW8XgCOFJHhwQUiMhA4Fni+/YoiMlBE2sd5JfBghH1OAJqwrlgqpVKcW8nraawa1s9FJCOQnO7Hutr42+BKInI6sA2YE7b9l0XkpHbrXQlcAvwilg6uSinvc+VqozGmWUQmY90etBqrh/1K4Jyw5FMLHAC2t1s2D3gA+I2I+IC+wD7gemPM3ETEr5Ryn2v3NhpjdmJdWexsnY+BgrBlO4EfB36UUocpHRJHKeVJmryUUp6kyUsp5UmavJRSnqTJSynlSZq8lFKepMlLdfDII4/w4osvJvSYLS0t/PnPf2bKlClMmDCB4447jokTJ/LHP/6RSAMHLF26lLPPPpvjjz+e0aNHc8stt9DY2Om9+CoFafJSHbiRvJYuXcpVV13Ftddey9KlS1m1ahW33XYbX/va10LTqwWtX7+eSZMmcdlll7Fy5UoWL17M66+/zrXXXpvQmJX7NHmppHDaaafx5S9/OfT8iiuu4IwzzuDRRx/tUPv62c9+RkFBAd/97ncB6Nu3L3fddRd/+ctf+OCDDxIet3KPJq8U9/LLL1NaWoqIcOeddzJ79mwmTJjA0KFDueOOO0LrbdmyhdLSUrZt2xbaprS0lLfecn4yplNOOYX58+cfsnzw4MHU1dXR0tICQGtrK6+++ipnn312h3HwzznnHACef/75Q/ahUpfO25jipk+fzvTp0xERnnnmGZ5//nl+/vOf88Ybb3D++edz9tlnM2XKFIYNG8ayZcsYPnw4ZWVlPP30013ue8mSJVx33XVdrjdx4sTQ8M+RiEhoPPz21q1bx6mnnkpmZiZgTZFWV1cXGn46qLCwkF69erF8+fIuY1GpQ5NXrO52buqxbrn7QI83LS0tZcKECQBMmTKF/Px8ysvLmTJlSo/2N3HiRJYtW9bjeDrz/vvvs3z58g41st27dwNEHEm1d+/e7Nmzx5FYVHLS08bDyKhRozo879evn+2Ty9qhtraWr3/96/zkJz/hrLPOimkbN4YzV+7Smles4qjxJItIk8u2tbVFWdsdzc3NXH755UyZMoXbb7+9w2v9+/cHiDgBR01NDYWFhQmJUSUHTV6qx+xq8wpqbm7msssuY+zYsTz44KGD5R511FHk5eWxadOmDsv37NlDTU0NJSUlMceuvE+Tl+rA5/OFTsE+++wzKisrOe200yKua2ebV7DGNXLkSB5++OHQ8m9+85vcfffdDBo0iIyMDKZOncqCBQswxoSuOAbbxS677DJbYlHeoMlLdTBixAi2bt0KwNy5c9m+fXvU5GWX5uZmrrjiCjZs2MCVV17Js88+G3pt4cKFNDU1hZ7/8Ic/5Oyzz+bXv/41N954IwcOHODee+/lS1/6EieffLKjcaokY4w57H5GjRplOrN69epOX0921dXVod8XLlxoxo8fbwAzcOBAc8MNN5h9+/aZ8ePHG5/PZ/r162fOPffc0PrvvPOOGTNmjDn++OPNKaecYtauXet4vC+99JLBGgo84s/GjRs7lO2DDz4wZ511lhk7dqwZOXKkuemmm0xDQ0OXx/HC+zp//ny3Q3AUsMTY9D3WmleKO/PMMyOe2kU73TvttNOoqKhwOqwOpk+f3q2rhRMnTmTBggUORqS8QLtKKKU8SZOXUsqTNHkppTxJk5dSypM0eSmlPEmTl1LKkzR5RdGdS/cq+en7mXo0eUXg8/loaGhwOwxlo4aGhohjhinv0uQVQVFREZWVldTX1+t/bI8zxlBfX09lZSVFRUVuh6NspD3sI+jduzcA27ZtCw1B7CWNjY1kZ2e7HYYjelI2n8/HwIEDQ++rSg2avKLo3bu3Zz/s5eXlnHDCCW6H4YhULpvqHj1tVEp5kmvJS0SKRORPIrI28POciAyNcVufiPxYRNaIyEoReVdEznA6ZqVU8nAleYlIJvAmkAkcB4wF6oD5IpIfwy4eA64EzjTGHA/8AXhTREodClkplWTcqnl9FSgBZhtjWo0xbcBs4CjgW51tKCKjgZnA/caYXQDGmCeBDcB9jkatlEoabiWvy4HNxpgNwQXGmB3A6sBrnbkUECB8ltK3gSkx1tyUUh7nVvIqATZGWL4RGBfDtn5gc4RtM7BOQZVSKc6trhL9gaURllcDuSKSY4yJ1sW9P1AfONUM3xYg4vxXIjIT63QToElEVnYzZi/pD+x2OwiHpHLZIPXLN9quHSVbPy9xaltjzFxgLoCILDHGTIzjWEktlcuXymWDw6N8du3LrdPG3cChc7Zby+o7qXUFt80VkfQI2wLonO9KHQbcSl7LgeERlo8AVsSwbRowLMK2rUBiZ49QSrnCreT1AnCkiAwPLhCRgcCxwPPtVxSRgSLSPs5/YE2JVRa2z0nAG8aYQ+eCP9Tc7ofsKalcvlQuG2j5YiZujJoQ6KS6BKuWdBXW1cPfA2cAJxhjagPrnQ4sBOYaY77VbvvHsZLV6caY3SIyA5gDnGqMsWcKZ6VUUnOl5mWMaQYmA21YfbsqgN7AOcHEFVALHAC2h+3iRuDvwDuBq4bfAKZo4lLq8OFKzUupSERkEPAUcL4xJp4rz0knlcvmlpQZVSLVb/SOs3ybRGRZhJ/znI47ViJyKfAecHQPt58lIqtFZLmIfCgil9gbYc/FUzYRKQ+UK/y9u8b+SLtPREpF5AkRWSoiHwdi/ZWIDIhh2/i+d8YYz/9g3eD9MdapZAaQDvwR+ATIj2H7x4F1wIDA8+uABqDU7bLZVL5NbpchhhgXAyOBp62PZbe2vQ2rC83RgeeTgRZgqtvlsqFs5cBwt8vQSXxrsC6y5QWeDwksWwfkdLFtXN871wtv0x/wG1hXII9qt+wIrDa1H3Sx7WisCwYzwpavAv7ldtniLV9g3U1ulyGGGDMCj936ggN9sUYkuTds+b+AVW6XK56yBbbxQvI6JmzZ1wOf18s72S7u712qnDam+o3e8ZTPE4wxrT3c9AIgl8jv31gRGRNXYDaIo2xeUGKMWR+2bFvgsV8n28X9vUuV5JXqN3rHUz4AROQXIrJERNaJyBsiMt3WCN1TEngM//tsDHvdy24SkfcDbUMLReRatwMKMlbPgXCjsGpeCzvZNO7vXaokr/5ApM6poRu9u9i22zd6J1g85QOoAj4CTsca/PEl4CUR+Y6tUbqjf+Ax/O+TTO9fPPYD67H6NR4H/Ar4rYj80tWoogjctjcD+L0xZl0nq8b9vUu2G7Pt5tiN3kkiphiNMSeHLZojIhcCPxWRJ40xjfaH5jovvH9dMsaEXzV9TkQmAd8XkV8ZY8JrLm67E+s2ve/3cPuY37dUqXml+o3e8ZQvmsWB7Y+LJ7AkEBw+Jvzvk0zvn90WY313T3I7kPYCp7NfxLrKW9vF6nF/71IleaX6jd49Lp+I5ERp/AxW18M/PF6zPPA4PGz5iLDXPUdEMkWkT4SXku69E5H/AW7GukumKoZN4v7epUrycvtGb6fFU74rgQcj7HMC0IR1xdIzRKQwcG9s0GtAPZHfv9XGmDWJii1eEcp2GvC3CKtOCDx+5HxUXRORq7HmoDgvcBUcEbkoMABocB37v3du9xOxqa9JJlYm/ytWO14a1q0YHTpxYjVYtwG/Ddv+cWAt0D/wfAbJ10m1R+UDvobVCHpSu2VXYl3puTdRZehGWZ8mSl8orP/KjcC8sOW3AbsI9IMDziOJOqn2tGyBL3YrMC1sWR3w/9wuTyCeqwLflVuAq9v9/A6420T5XAaWx/W9S4kGe2NMs4hMBh7GqkkYYCXdu9H7f7Fu9G7BunKVNDd6x1m+ecADwG9ExIfVqXMfcL2xRpdNCiLyAFbP+OLA8+Df/mRz8HJ8A7CXg/2IADDG3C8ijcArItKK9UX5gjFmXkKC70IcZfsQuBX4oYj8FMgDmoGfYL2nyeAxIJvI8dwTeHTke6c3ZiulPClV2ryUUocZTV5KKU/S5KWU8iRNXkopT9LkpZTyJE1eSilP0uSllPIkTV5KKU/S5KWU8iRNXkopT9LkpZTyJE1eKqkF5qqsFhG/iLwVWDZHRPaJyEYRuc7tGJU79MZslfRE5AtY41p9wxjzpIgcCfwTOM10PWKnSlGavJQniMgLWON0lQJ/AH5qjHnD3aiUmzR5KU8QkSOwxjJrA14xxiTN9F/KHdrmpTzBWMML34M1ZVb4RKXqMKQ1L+UJgfHPy4EcrBFJjzPG7O50I5XStOalvOJ7WFN+XYI17PCj7oaj3KY1L5X0RORo4Dmsq4sNIvJNrMkbLjbGvOJudMotWvNSSU1E7gP+AxyBNbsMwA2Bxz+JyHOuBKZcpzUvpZQnac1LKeVJmryUUp6kyUsp5UmavJRSnqTJSynlSZq8lFKepMlLKeVJmryUUp6kyUsp5UmavJRSnvT/AXOngZpjWceFAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "linear_convection_cfl(85)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAS8AAAEbCAYAAACLNQJjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl8VOW9+PHPN8mEJIQEkrDJImgFQY1pQa1YBRewauW6Xtu6tFKL1lZr1VutXr0uXWy91drWXn900db21rZal2tFrS0Rqy0ICKhsIiAQ9iUb2TPf3x9nZgjDJJnlTM6c4ft+veY1yZnzPPN9MnO+ec5znnOOqCrGGOM3OV4HYIwxybDkZYzxJUtexhhfsuRljPElS17GGF+y5GWM8SVLXsYYX8rz4k1FpAr4KvCJUAwB4DXgflXd2UvZDUBtjJduVdXXXA7VGJOhxItJqiKyCngfuEpV94nICOBvOD3B41W1uYeyG1R1TN9EaozJVF7uNt6mqvsAVLUGeBA4CjjXw5iMMT7hyW4jUKmqbVHLtoSeB/V1MMYY//Gk5xUjcQGMAxSY31t5EfmBiCwSkTUi8qqIzHQ9SGNMRvOq53UAEckFZgG/VNU1vay+A3gHuAsIArOB50XkBlX9aQ/vMTu0LgUFBZNGjx7tSuyZKBgMkpOTnQeSs7ltkP3tW7NmzS5VHexGXZ4M2B8UhMg9wPnAVFVtTKL8X4BTgSGq2tLb+uPHj9fVq1cnHKdfVFdXM23aNK/DSItsbhtkf/tEZLGqTnajLs9TvIhcDfw7cE4yiStkATAAOMa1wIwxGc3T5CUiVwK3AGeo6o441i8UkeIYL3WGnnPdjM8Yk7k8S14icgVwG3CWqm4LLftMaGwqvM5QEeka42XAD2NUNwloBVakMWRjTAbxJHmJyOXAz4EngLNE5IpQMjsfOCy0zik40ycejSr+ORE5oUtdlwEXAD9IYbfTGOMzXh1t/AlQgDMxNdq9oedGoA7Y2uW1uaEyPxORADAQ2Atcp6pz0heuMSbTeJK8VLUsjnWWAWVRy7YD94cexphDmOdHG40xJhmWvIwxvmTJyxjjS5a8jDG+ZMnLGONLlryMMb5kycsY40uWvIwxvmTJyxjjS5a8jDG+ZMnLGONLlryMMb5kycsY40uWvIwxvmTJyxjjS5a8jDG+ZMnLGONLlryMMb5kycsY40uWvIwxvmTJyxjjS5a8jDG+ZMnLGONLlryMMb5kycsY40uWvIwxvmTJyxjjS5a8jDG+ZMnLGONLlryMMb5kycsY40ueJS8RqRKRn4vIYhFZJiIrROTHIjI4jrIBEblfRFaJyHsi8paIfKov4jbGZAYve15PAWXAaap6PDAdmAG8KSKFvZT9CXAZcKqqHgv8CviriFSlM2BjTObwerfxNlXdB6CqNcCDwFHAud0VEJHxwGzgAVXdGSr7C2Ad8J20R2yMyQh5Hr53paq2RS3bEnoe1EO5CwEB5kUt/ztwnYgUq2qjSzH6zq7GVrY2Bvlwp/MnyMsRDi/vH3l9w659dKrGLDuoKJ+y/vkANLZ2sL2+pdv3ObysiLzcHFBl265dtNbvhmA7op0HrFcYyGPwgH4AtHcG2VrX3G2dFcX9KMp3vpK1ze3UN0d/PaCjdhPs+mD/AslxHvnFUFQGObnd1m+yi2fJK0biAhgHKDC/h6KVQBDYGLV8PU57JgIL3YjRbxZ/tJdLH3uLoAL/eB2AYSUF/OuOMyPrXPLYW+xqjPWnh5unj+PGM48C4M21u7j2ycXdvJOy7LJOSj94Btb+jWFt8f2vCACj42zLwNAj2miApd0UkhwYdhwcezF8/EonmZms5WXP6wAikgvMAn6pqmt6WLUCaFKN+hcP9aHn8m7qn42zu8ngwYOprq5OLeAM9PqmdoIKBbnKwH7OiEBJbtsBbR2Y10FekcQsv7NmA9XVNQCs3dXJsFjrqXKj/o7S51+KLGoln1qKaSOAIij7y/XLhZJ85/dOhT0tsXt9AKX5Qn6o47SvHZo6Dlw3qM6jJF/oH3CWiQYRDZLb2UygowG2LoOty2j6x2Msrfoubf166sRnnsbGxqz8bqZDxiQv4C6gA/hGkuVjb5EhqjoHmAMwfvx4nTZtWpJvk7k2/XMDvP8+Jx8W4PGvnh1znXibPQ34aqwX/v5tmP8S5ATg9DvguEvpVzqSodLjnz9iQHxvT1GMZSu31jN3/kIuOWsKg8tjrNHeDGtfg3nfo2jH+0xZ+3344l+gf8z/ZxmpurqabPxupoPXA/YAiMjVwL8D58QxXrULKAr11LoKbxe73Y7PLxRnjCsvXZ/qxn/B/AdBcuHSx+HUm2HgKIgzcaVqwvASPjE0j9GxEhdAoBAmnA9feAEGHw07V8Jrd/dJbKbveZ68RORK4BbgDFXdEUeR5Thxj4paPhan57bS3Qj946qTx7D2u+fy+aPz0/MGC3/uPJ9yo5MkMlX/Cvjs/zo/v/s0NO3xNh6TFp4mLxG5ArgNOEtVt4WWfSY0PhVeZ6iIdI3zWZxOxrSo6k4HXlXVhvRGnfkkHT2hxh2w4nlnUHzyLPfrj8Pij/bw/No23vhgZ+8rlx8JR54JHS2w9H/TH5zpc17OsL8c+DnwBHCWiFwRSmbnA4eF1jkFZ/rEo+FyqroaZ+zqWyJSEVpvFnAkcGdftuGQ8s6TEGyHcZ+GgfEeM3TXog17eXZtO298sCu+AidcEyr4SwgG0xeY8YSXA/Y/AQpwJqZGuzf03AjUAVujXr8B+C+c2fjtQAMwQ1W7O4h+SPjlP9bzp0WbOLmi/aBuaUqCQVj0uPPzCV9ys+aE5OU6/2vbOuJMROPOhtJRsGcdrJsHHzuz9zLGN7yc59XrJBxVXYZzClH08nbgP0MPE7K9voVV2xqoLAm4W/Gu1VC3CQYMhyPOcLfuBOTnOrvDHfH2onJy4fjPOgcZ1v7NkleW8XzA3rinvdPZqHNzXB7z2hSa8zv6k5Dj3Vcm3PNq7+h+rthBRp/sPG8+JOctZzVLXlkkkrzcHq8Pb/gjT3S54sQEwsmrM4Hxq5GTAXEmr3a0picw4wlLXlmko9Ppkbg+z2vT287zKK+Tl5OV24MJ9LwKSp05X51tTgIzWcOSVxZpS0fPq3mvM+aV2w+GVbpYceKK8vMoyYfifgmefD3qBOd5k+06ZhNLXlkk3PNydcxrc+jk7OHHQ16aJr/GafrEofz4jP5876IEk+jIUPKyca+skknnNpoUnXpUBaWFAYbLdvcqDW/wHu8ypiQ8Vhfe/TVZwXpeWeTSyaO4/4JjGVvq4jWtwrta4d6LH1WMc8a+GrZA3WavozEuseRluqcKW95xfs6A5PXOxr18Y14T1/x6UWIFc3JgxCTn53B7jO9Z8soiq7bVs2TjXpoTmQfVk6Y90FIL/Uqg5DB36kxBUGFvq7KrMYkpD4MnOM+717oblPGMjXllkbufe5+FG/Zw+4kF7lQY3tDLj+yzy970JJDoDPuuyo90ni15ZQ3reWWR9qDLUyUiyetjLlWYmkAyM+zDwm3Y/aGLERkvWfLKIuGZ565NUs245BWepJpMzyucvKznlS0seWWRyDyvbO95JXJ6UNiA4RAogn07obnW5ciMFyx5ZZH9PS+Xsld4Fys8XuSxcPIKJ+mE5ORAWagde2zXMRtY8soi7W72vILB/Rt5WWYkr5LCAOcfEWDWKWOTqyAyaG/JKxvY0cYs0hG5JI4LldXXOJdQLh4KBSUuVJi64n55XDwun2mnHZFcBTbulVUseWWR315zEi3tQbas6u5msQnIsPEuV1jyyiq225hFjhhczMTDStwZ8+o6xytDBIPKsp0dvLYiyXM3LXllFUteJrbIYH3m9LwUeHhxK19+chGqycz16jLmlUx5k1FstzGL3PqnZQhwdrkLG2YG7jbm5giCk3c6g0peokcmisqgsAya90DDNigZnpY4Td+wnlcWeWbJZv60eDOuzJTYu8F5HpTkkb00CU/A7UjkaqpdlYXaU/uROwEZz1jyyhKdQUUVcgRyUj0PUXX/pWMGRt+Y3FvhzlZbMhNVwbkVGkDtJncCMp6x5JUlIhNU3ZgnsW8XdDRDwUDoNyD1+lwU6XklM1EV9ifjuo3uBGQ8Y8krS4STV8CNfcbwhp1hvS7Yf4nrpE4RAigN3e271pKX31nyyhLh2fUBN87KDu9ShTf0DJIX3m2M967Z0QbabmO2sKONWaIjcl6jG8krc3te//nJAk45ZQrl/fslV0F4zKvOkpffWfLKEjk5wuTDB1FaGAD2pVZZeMMemHk9r0EFOQwZkMLFFrv2vFQz4iKLJjmWvLJERXE/nv7KFACqq6tTqyyy25h5Pa+UFZRCv1JorYOm3dC/wuuITJJszMscLNLzyrzk9dSqNq5+fCEf7mxMvpJI78sG7f3M0+QlIsNF5GURsXM1UtQZVFraO+lMdvJmVxk8YL+2tpN5q3dS29SWfCU27pUVPNttFJELgYeB9iTKVgNDgOhv8EOq+pvUo/Of92rq+LdH3+S4EaXcclwKFbXUObtUgSLndJoME5mkmsodkuyIY1bwcszrdmA6cCeQzAl056rqBlcj8rHwHXUSPt8vWtfxrgwczA7PBEl6nhdYzytLeJm8TlHVDsnADcSPwj2RQKoz7DN4vAsgV1K4/VmYjXllBc/GvFS1w6v3zkbhjTngZs8rA4Vzc0q7jZFZ9tbz8jM/H228WUQWisgqEZkvIld7HZCXwuf6pd7zytwJqrB/zMuVnped3+hrfp3nVQusBb4FtAAXAr8VkWNU9dZYBURkNjAbYPDgwanPhcowS7Y7Hdm6vXtoHNyRdPsmfrCEIcCKLY3syMC/0YjCDjoG57F57Uqq96xJrhJVTs3JJ7eljjdee4nOvCJ3g0xBY2Nj1n0308WXyUtVL4ha9LSInA58Q0R+rKoH/UtV1TnAHIDx48frtGnT0h9oH9q3fCu8s4RhQwZTXNxA0u374D4AJp58NhNHf9K9AF1TnXzbunpvNOxey6nHHQ5Dj0m9PpdUV7vUvkOAn3cboy3Aac8JXgfihcqRpfzw0uO58uTDU6sow8e8XDPQxr38znc9LxHJBwpVtS7qpc7Qc24fh5QRRpUVMarM2f2pTnZ7bG+BfTsgJw8GDHMvOBc1dyjb6looLsijuF8KX1+bLuF7Gd/zEpHyUMIKmwL8Mcaqk0LP76Q/qiwVvnpqyQjIycz/AX/+oI1Pfu9vPLUwxcF2my7hexmdvERkLFADPB/10pkicl6X9aYB1wJPquoHfRdh5nivpo5f/WM9C9btTr6SyJHGzDstKCwvcjHCFE+DCk+XsJ6Xb3mWvETkQRFZCswM/b409Ojay2oG9gBbuixbAnwTuENElonIWuBnwLeBWX0Tfeb517rd3PfiCl5+f1vylfhgvCt8McKOVGbYg50ilAU8G/NS1f+IY51twGFRy+qBh0IPExK+m05+KvO8ajO/55XrxulBYGNeWSCjdxtN/No7XDi3McNPDYL9Pa/2VK+eMWA4SC40bncOVBjfseSVJcIbc0qXgfbBbmPkBhzJXsM+UlGec2ACoL4mxaiMFyx5ZYnwblR+Kjfg8FHPK+mbznZlRxx9zXfzvExs+2/AIZDMdt3ZAfVbAIGSka7G5qbjh+Ry+omTGV3uwik9A0fDR2/auJdPWfLKEiJCfm6O0/NK+PKOQMMW0E5nLCgvv/f1PTKkKIdpE4e6U5ndPdvXLHlliTvOncAd504AoLr6o8Qr8MF4l+sG2hFHP7MxL+PwwXgXQE1DkIf/uoYXlm3pfeXeWM/L1yx5GUd40DrDe15b9gV55G8f8NLyralXFjk52wbs/ch2G7PE/S+u4I0PdkZ2HRO2+0PnuewI94JKA1euYR9WOsqZ61W3yZnrFUjhZramz1nPK0vU7G1mzfZGmto6e185lt1rnefyZO6F0ndy3ZqkCs6BiUGHAwp716den+lTlryyRLgnktRloFVhd+h89gxPXnluTVINC7c3nLyNb7iavETkZ27WZ+IXmWGfzOlBTXuc+zXmD4DiIS5H5q5Iz8uN3Uaw5OVjCY15ichVvaxybgqxmBSEJ6kGcnJI+LZMkV3GIzPyXo1dRU7MdmO3EZw2gyUvH0p0wP6JHl5z6dtkkrF/t1FSSF6ZvcsI0C9XGDygHwMLA+5UGOl5fehOfabPJJq8VnJw76o/cDTweeBRN4IyiQtfnC8vmTEvHyWvUQNyePvOM9yr0HYbfSvR5HWjqsaavr1CROYCTwHzUg/LJOozlcM5fmQpQ0v60ZBoYR8lL9cNOAzyCmHfTmiuhcKBXkdk4pRQ8lLVv/XwWrOIHJ16SCYZ15y6f35Wwn2I8C5TePznUJKT47R7+3uw50MYMan3MiYjJDpgf1qsxcAg4AKcG8AaPwkGnY0WfJG8aluCnPid1yjrn8/LN8X6OiYhnLx2W/Lyk0R3G6uJPTAvwGbgilQDMslZsnEvqsqxI0oTK1hfAx0t0H8IFCRY1gMiwo6GVjrdOtoINu7lU4kmrw+Ba6KWdQI7gA9VNcnp3SZV1/92CdvqW3jr9gQHs3023uX6PC+w5OVTiSavR1X19bREYlLSEUzyGvbb3nWeB49zOaL02H9uo4s9r4rxznP4b2F8IaHj6qr6o3QFYlLTFjpdJuG7B21+23keeYLLEaVHOHmFk7Urhh0Huf1g1xpo3utevSat7NzGLNERTGKel2qX5HViGqJy3/7dRkXVpd5XXj4cVuX8vHmxO3WatLPklSU6QrtRgUR2G+s2Q8NWKBjomzEvEYmcnO3KTTjCwj3PzQvdq9OklV3PKwuoKm1dzm2MW3hDHTnZme/kEzeeeRSun4EZTl6bLHn5hSWvLBCeNpAjkJOTwGa9yV+7jGE3nnmU+5WOCv0NahY7c998lMwPVZa8skCOCK/dPDXxuU/hntcofwzWp1XJYc4t3+o3w85VMHSi1xGZXti/lyyQkyN8bEgx44cNiL9QewtsXQ4IjJicttjS4c21u3j5vW00J3vV2O6MsnEvP7Hkdaja8A8ItsOQCVBQ4nU0Cfnm08u57reL2dXY6m7Fo05yntd2ewqvySCWvLJAbVMbX3/qHe77vxXxF1r0K+f52IvTE1Qa5Ycme7k6yx5gwkznhhyrX4KGbe7WbVznafISkeEi8rKI2IUMU9DQ0sHzS7fwyvtxbnC1m2DNXMgJwCd6uzhu5knLVAmA0hEw/hwIdsCS37hbt3GdZ8lLRC4E/gkkdSkDEblJRFaIyHIRWSIiF7gboX+EN+K453gt+TVoECb+W8Zfsz6W8E1G2ty6CUdXJ4RO3V38BHQmfE1a04e87HndDkwH3ky0oIjcDvwncL6qVgK3AX8SkXPcDdEfwrtPcc2ub9zpbJiwf0P1mXCSdr3nBTB2qjNht74Glv3e/fqNa7xMXqeo6geJFhKRgcBdwM9U9UMAVf0r8Crw3+6G6A9x3/asaQ/85t+cq4aOPAFGf7IPonNfOEm7PuYFzvyuKTc6P794E6ye6/57GFd4Ns9LVZPtk38aKOLgy03/HfhvETlaVVf1VEGwo5X17y+I+dqw0kIKA7kA7G5spb6lPeZ6ebk5jBpUFPl9/a59dHcPkrL+/SgN3TCioaWDXY3dX7Px8LL+kYmmW2qbae2IPR2gf78AQwb0A6Bzay3jZSMf0wGwbRD9Gzc4V0gIdkBHq5OsNi+Cd//k9CgqxsFnf5/xdwrqTrjn1d4RZF9rB5v2NnW77tiK/vTLcz7PmtpmGrr5PIsCeYwudz7Pzqorqd24gvJlj6F/uJKGMTNoPPwsOvoPI5g/gPKS/gwoyAegtrmdvU1tMesUEcaU94/8vmlvU+QuT9FKCgKUF/ejf+MGmjYtZ3t9c7dtGjGoKHIC/vaGVppaY7epIC+P4QOdu4B3BpWNe/Z1W2fFgAIG9HPSQW1Tm6ttAmhq6+yxTcnw4yTVytBz9C2O13d5vcfkVdq8ibF/mtHrG5WHHvEYG+d6A0KPeBwW53qVwCv9gFrgMTgBYFE3Kw+eAFc+C8WD46w98/zs8kkEVSktDLD4o718ds6/ul33b7dM5cjBxQB8f+4qXli2JeZ6J44t44/XngxAS0eQSQtO5Vt5m7g27y+UrHMesQwMPeIxKo51wp9dvN+noXGul0v8dbrdJnB6G/G+f7z8mLwqQs/R95moDz3HzDciMhuYDXDs8AI+7ObPXlYgkf/sDW1BmmL/UyMvB8oL9++m7dgX7PbebwPyhaKAU2dzh1Lf2v1YzdAiifSI9jQHae9mz6gwD0r6Oe/fEVR2N2vkfYLBIDk5OajkEswJ0B4ooaVgKDsHT6GudAIsWQ2s7jaGTNbY2Mjyt9+K/L6+rpORxd33IJe8vZBNRaG/U31bt+sGWuuprq4GoLVTGVmcw5NcwSt6DmfqP5mgaxmk9RTRQlm/IIWhLaepXWloi/15isCQov3fkV3NQbrb0y0KwID8HILBTjo1hz0t3X9HKgqF3FDvvK41SEs3+zD5uTCowHl/VWVHU/d1lvYTCvIkLW0C5yogTpve7zaGhKmqpw+ce0FqAuvPwdk/K49aPj20/Cu91TFu3DjNZvPmzfM6hLTJ5rapZn/7gEXqUu7w4yTVXaHn6L2v8O+7+zAWY4xH/Ji8loeex0QtHxv1ujEmi2V88hKRchHJ77LoZaAJmBa16unACu3lSKMxJjtkdPISkbFADfB8eJmq1gL3A18VkSNC650FnA3c6kWcxpi+59nRRhF5EGeQfXTo96Whl05U1fAkk2ZgD3DA8W1VfUBEWoAXRaQD5/Zrl6qqzSg05hDh5STV/4hjnW10M91JnTsZ2d2MjDlEZfRuozHGdMeSlzHGlyx5GWN8yZKXMcaXLHkZY3zJkpcxxpcseRljfMmSlzHGlyx5GWN8yZKXMcaXLHkZY3zJkpcxxpcseRljfMmSlzHGlyx5GWN8yZKXMcaX/Hjfxj5RX1/Pjh07aG/v5saNGay0tJSVK1d6HUZaJNO2QCDAkCFDKCkpSVNUxguWvGKor69n+/btjBgxgsLCQkS6v6lpJmpoaGDAgHjvy+0vibZNVWlubqampgbAElgWsd3GGHbs2MGIESMoKiryXeIyBxIRioqKGDFiBDt27PA6HOMiS14xtLe3U1hY6HUYxkWFhYW+HAIw3bPk1Q3rcWUX+zyzjyUvY4wvWfIyxviSJS+TMR577DEmTpyIiPDEE094HY7JcJa8DgE7duygqqqKsrIyRISqqip+8YtfJFTHbbfdxuTJkw9YtnTpUu655x5qa2sPWN7a2srYsWN55JFHEnqP6667jpdeeimhMubQZcnrEDBkyBCWLl3KzJkzASfpXHPNNQnXMXr06AOWLV26lHvvvfeg5JWXl8fo0aMpLy9PLXBjemCTVE1cbrnlFm655Za41s3NzeX1119Pc0TmUGc9r0NYeFdyzJgxzJ07lzPOOIORI0cyffp0Nm/eHFnvq1/9KqNHj0ZE2LBhAwB33nknd999NwDnnnsuVVVVXHLJJezevZuqqiqKi4uZNm1apI6WlhZuv/12Jk2axCc+8QkqKyu57rrrDuq1GRMvS16HsPCu5N69e/nnP//J3//+d1atWsXGjRv55je/GVnv0Ucf5b777jug7He+853IspdeeomlS5fy9NNPU15eztKlSw8aH6utreXxxx/nueeeY8mSJSxcuJA9e/Zw1VVXpb+hJivZbmMCxtz+l25f++6Fx/H5k5wxof9dsJE7nn2323U3PHBe5OfP/OQN3qupj7ne504cxfcuqgTg3c11HDeyNJmwe9XQ0MBNN90EQHFxMdOnT+fPf/6zq+9RUVHBW2+9xahRowAoKChg1qxZnHPOOWzfvp2hQ4e6+n4m+3mWvERkCPAwEP4X/S5wk6pu7r5UpOwGINb+xq2q+pprQR4iKioqKCsri/xeVlbG9u3bXX2PvLw81qxZw/XXX09NTQ15eXk0NjYCsG7dOkteJmGeJC8RyQf+CqwBjgEU+BUwT0Q+rqqNvdWhqlXpjfJgXXtMPfn8SaMjvbDevHjDqXGtl65eF0BRUdEBv+fk5BAMBl19j7lz53Leeefx0EMP8fWvfx0Robq6mtNPP53W1lZX38scGrwa8/oCUAncpqodqtoJ3AYcAXzFo5hMGj355JMUFxdz00032XmGxhVeJa+LgY2qui68QFW3AStCrxkfCAQCgHPNLIBXXnmFPXv2xFy3tbWVnJwDv27btm1Lb4Amq3mVvCqB9TGWrweOi6cCEfmBiCwSkTUi8qqIzHQ1QtOrsWPHArB582YaGhq48MILaWhoiLnueeedR11dXWRmf0NDAz/60Y/6LFaTfbwasK8AFsdYXg8UiUihqjb3UH4H8A5wFxAEZgPPi8gNqvrTWAVEZHZoPQYPHkx1dXW3lZeWlna7EfpBZ2fnAfHv3LmTCy64gE2bNgFQWVnJ7Nmz+eMf/8i7777Lvn37qKys5IUXXuCBBx7g2Wefjaz30EMP8Yc//IGXX34ZgE9/+tN87Wtf4wtf+ALHHXccV111FVdccQWFhYXcdNNNqCqVlZWsW7cuUsdzzz3HJZdcwvr167n33nt5+OGHGTZsGNOmTWPBggXMmjWL6667jvz8fB577DEA7rrrLt54442DElx02xLR0tLS4+eeCRobGzM+xkwh4S5/n76pSBvwiqqeH7X8d8DngaJeklesOv8CnAoMUdWWntYdP368rl69utvXV65cyYQJExJ5+4xil4GOzQ+fa3V19QGTe7ONiCxW1cm9r9k7r3YbdwGxvoEDgKZEE1fIglD5Y1IJzBjjD14lr+XAmBjLx+LM9+qWiBSKSHGMlzpDz7mphWaM8QOvktefgcNFZEx4gYgMBSYAz3RdUUSGikjXOC8DfhijzklAK84RS2NMlvMqeT2B08P6vojkhZLTAzhHG/8nvJKInAJsAR6NKv85ETmhy3qXARcAP4hngqsxxv88Odqoqm0iMh3n9KAVODPs3wPOiEo+jUAdsLXLsrnAg8DPRCQADAT2Atep6py+iN8Y4z3Pzm1U1e04RxZ7WmcZUBa1bDtwf+hhjDlE2SVxjDG+ZMnLGONLlrzQGWXTAAAMEUlEQVSMMb5kycsY40uWvIwxvmTJyxjjS5a8zAF+9KMf8dxzz/Xpe7a3t/P73/+eGTNmMGnSJI455hgmT57Mr3/9a2JdOGDx4sVMnTqVY489lvHjx3PrrbfS0tLjufgmC1nyMgfwInktXryYyy+/nKuvvprFixfz/vvvc/vtt/PFL34xcnu1sLVr13L66adz0UUX8d5777FgwQJeeeUVrr766j6N2XjPkpfJCFOmTOFzn/tc5PdLLrmET33qUzzyyCMH9L6+973vUVZWxo033gjAwIEDufvuu3nqqad4++23+zxu4x1LXlnuhRdeoKqqChHhrrvu4rbbbmPSpEmMHDmSO++8M7Lepk2bqKqqYsuWLZEyVVVVvPZa+m/GdNJJJzFv3ryDlh922GHs27eP9vZ2ADo6OnjppZeYOnXqAdfBP+OMMwB45plnDqrDZC+7b2OWmzlzJjNnzkREePLJJ3nmmWf4/ve/z6uvvsrZZ5/N1KlTmTFjBqNGjWLp0qWMGTOGadOm8cQTT/Ra96JFi7jmmmt6XW/y5MmRyz/HIiKR6+F3tWbNGk4++WTy8/MB5xZp+/bti1x+Oqy8vJwBAwawfPnyXmMx2cOSV7zuSd+txxJyT13SRauqqpg0aRIAM2bMoLi4mOrqambMmJFUfZMnT2bp0qVJx9OThQsXsnz58gN6ZLt27QKIeSXVkpISdu/enZZYTGay3cZDyLhx4w74fdCgQa7fXNYNjY2NfOlLX+Lb3/42p512WlxlvLicufGW9bzilUKPJ1PEurlsZ2dnN2t7o62tjYsvvpgZM2bwrW9964DXKioqAGLegKOhoYHy8vI+idFkBkteJmlujXmFtbW1cdFFFzFx4kR++MODL5Z7xBFH0L9/fzZs2HDA8t27d9PQ0EBlZWXcsRv/s+RlDhAIBCK7YB999BE1NTVMmTIl5rpujnmFe1xHHXUUDz/8cGT5tddeyz333MPw4cPJy8vjnHPO4fXXX0dVI0ccw+NiF110kSuxGH+w5GUOMHbsWDZv3gzAnDlz2Lp1a7fJyy1tbW1ccsklrFu3jssuu4zf/va3kdfmz59Pa2tr5Pc77riDqVOn8tOf/pQbbriBuro67rvvPj772c9y4oknpjVOk2FU9ZB7jBs3TnuyYsWKHl/PdPX19ZGf58+fr8cff7wCOnToUL3++ut17969evzxx2sgENBBgwbpmWeeGVn/zTff1KOPPlqPPfZYPemkk3T16tVpj/f5559XnEuBx3ysX7/+gLa9/fbbetppp+nEiRP1qKOO0ptvvlmbm5t7fR8/fK7z5s3zOoS0AhapS9ux9byy3Kmnnhpz16673b0pU6awcuXKdId1gJkzZyZ0tHDy5Mm8/vrraYzI+IFNlTDG+JIlL2OML1nyMsb4kiUvY4wvWfIyxviSJS9jjC9Z8upGIofuTeazzzP7WPKKIRAI0Nzc7HUYxkXNzc0xrxlm/MuSVwxDhgyhpqaGpqYm+4/tc6pKU1MTNTU1DBkyxOtwjItshn0MJSUlAGzZsiVyCWI/aWlpoaCgwOsw0iKZtgUCAYYOHRr5XE12sOTVjZKSEt9+2aurq/n4xz/udRhpkc1tM4mx3UZjjC95lrxEZIiI/E5EVoceT4vIyDjLBkTkfhFZJSLvichbIvKpdMdsjMkcniQvEckH/grkA8cAE4F9wDwRKY6jip8AlwGnquqxwK+Av4pIVZpCNsZkGK96Xl8AKoHbVLVDVTuB24AjgK/0VFBExgOzgQdUdSeAqv4CWAd8J61RG2MyhlfJ62Jgo6quCy9Q1W3AitBrPbkQECD6LqV/B2bE2XMzxvicV8mrElgfY/l64Lg4ygaBjTHK5uHsghpjspxXUyUqgMUxltcDRSJSqKrdTXGvAJpCu5rRZQFi3v9KRGbj7G4CtIrIewnG7CcVwC6vg0iTbG4bZH/7xrtVUabN85J0lVXVOcAcABFZpKqTU3ivjJbN7cvmtsGh0T636vJqt3EXcPA9251lTT30usJli0QkN0ZZALvnuzGHAK+S13JgTIzlY4F34yibA4yKUbYD6Nu7RxhjPOFV8vozcLiIjAkvEJGhwATgma4rishQEeka57M4t8SaFlXn6cCrqnrwveAPNifxkH0lm9uXzW0Da1/cxIurJoQmqS7C6SVdjnP08JfAp4CPq2pjaL1TgPnAHFX9Spfyj+Ekq1NUdZeIzAIeBU5WVXdu4WyMyWie9LxUtQ2YDnTizO1aCZQAZ4QTV0gjUAdsjariBuBPwJuho4ZfBmZY4jLm0OFJz8uYWERkOPA4cLaqpnLkOeNkc9u8kjVXlcj2E71TbN8GEVka43FWuuOOl4hcCPwTODLJ8jeJyAoRWS4iS0TkAncjTF4qbROR6lC7oj+7q9yPNHEiUiUiPxeRxSKyLBTrj0VkcBxlU9vuVNX3D5wTvJfh7ErmAbnAr4EPgOI4yj8GrAEGh36/BmgGqrxum0vt2+B1G+KIcQFwFPCE87VMqOztOFNojgz9Ph1oB87xul0utK0aGON1G3qIbxXOQbb+od9HhJatAQp7KZvSdud54136A34Z5wjkEV2WDcMZU/uPXsqOxzlgMCtq+fvAX7xuW6rtC627wes2xBFjXug5oQ0cGIhzRZL7opb/BXjf63al0rZQGT8kr49FLftS6Pt6cQ/lUt7usmW3MdtP9E6lfb6gqh1JFv00UETsz2+iiBydUmAuSKFtflCpqmujlm0JPQ/qoVzK2122JK9sP9E7lfYBICI/EJFFIrJGRF4VkZmuRuidytBz9N9nfdTrfnaziCwMjQ3NF5GrvQ4oTJ2ZA9HG4fS85vdQNOXtLluSVwUQa3Jq5ETvXsomfKJ3H0ulfQA7gHeAU3Au/vg88LyIfM3VKL1REXqO/vtk0ueXilpgLc68xmOAHwP/IyL/7WlU3QidtjcL+KWqrulh1ZS3u0w7MdttaTvRO0PEFaOqnhi16FERORf4roj8QlVb3A/Nc374/HqlqtFHTZ8WkdOBb4jIj1U1uufitbtwTtP7RpLl4/7csqXnle0neqfSvu4sCJU/JpXAMkD48jHRf59M+vzctgBn2z3B60C6Cu3O/jvOUd7GXlZPebvLluSV7Sd6J90+ESnsZvAz3F2P/vL4zfLQ85io5WOjXvcdEckXkdIYL2XcZyciVwK34JwlsyOOIilvd9mSvLw+0TvdUmnfZcAPY9Q5CWjFOWLpGyJSHjo3NuxloInYn98KVV3VV7GlKkbbpgB/jLHqpNDzO+mPqncicgXOPSjOCh0FR0Q+E7oAaHgd97c7r+eJuDTXJB8nk/8BZxwvB+dUjAMmceIMWHcC/xNV/jFgNVAR+n0WmTdJNan2AV/EGQQ9ocuyy3CO9NzXV21IoK1P0M1cKJz/yi3A3KjltwM7Cc2DA84igyapJtu20IbdAZwXtWwf8Buv2xOK5/LQtnIrcEWXx/8D7tFuvpeh5Sltd1kxYK+qbSIyHXgYpyehwHskdqL3f+Gc6N2Oc+QqY070TrF9c4EHgZ+JSABnUude4Dp1ri6bEUTkQZyZ8aNDv4f/9ifq/sPxzcAe9s8jAkBVHxCRFuBFEenA2VAuVdW5fRJ8L1Jo2xLgm8AdIvJdoD/QBnwb5zPNBD8BCogdz72h57Rsd3ZitjHGl7JlzMsYc4ix5GWM8SVLXsYYX7LkZYzxJUtexhhfsuRljPElS17GGF+y5GWM8SVLXsYYX7LkZYzxJUtexhhfsuRlMlroXpX1IhIUkddCyx4Vkb0isl5ErvE6RuMNOzHbZDwRuRTnulZfVtVfiMjhwP8BU7T3K3aaLGXJy/iCiPwZ5zpdVcCvgO+q6qveRmW8ZMnL+IKIDMO5llkn8KKqZsztv4w3bMzL+II6lxe+F+eWWdE3KjWHIOt5GV8IXf+8GijEuSLpMaq6q8dCJqtZz8v4xddxbvl1Ac5lhx/xNhzjNet5mYwnIkcCT+McXWwWkWtxbt5wvqq+6G10xivW8zIZTUS+A/wDGIZzdxmA60PPvxORpz0JzHjOel7GGF+ynpcxxpcseRljfMmSlzHGlyx5GWN8yZKXMcaXLHkZY3zJkpcxxpcseRljfMmSlzHGlyx5GWN86f8DLdAETvhXFc4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "linear_convection_cfl(121)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice that as the number of points `nx` increases, the wave convects a shorter and shorter distance. The number of time iterations we have advanced the solution to is held constant at `nt = 20`, but depending on the value of `nx` and the corresponding values of `dx` and `dt`, a shorter time window is being examined overall. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "###### The cell below loads the style of the notebook." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "css_file = '../../styles/numericalmoocstyle.css'\n", "HTML(open(css_file, 'r').read())" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (MOOC)", "language": "python", "name": "py36-mooc" }, "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.5" } }, "nbformat": 4, "nbformat_minor": 1 }