{ "cells": [ { "cell_type": "raw", "metadata": {}, "source": [ "Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Lorena A. Barba, Gilbert F. Forsyth 2015. Thanks to NSF for support via CAREER award #1149784." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[@LorenaABarba](https://twitter.com/LorenaABarba)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "12 steps to Navier-Stokes\n", "=====\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We continue our journey to solve the Navier-Stokes equation with Step 4. But don't continue unless you have completed the previous steps! In fact, this next step will be a combination of the two previous ones. The wonders of *code reuse*!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Step 4: Burgers' Equation\n", "----\n", "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can read about Burgers' Equation on its [wikipedia page](http://en.wikipedia.org/wiki/Burgers'_equation).\n", "\n", "Burgers' equation in one spatial dimension looks like this:\n", "\n", "$$\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} = \\nu \\frac{\\partial ^2u}{\\partial x^2}$$\n", "\n", "As you can see, it is a combination of non-linear convection and diffusion. It is surprising how much you learn from this neat little equation! \n", "\n", "We can discretize it using the methods we've already detailed in Steps [1](http://nbviewer.ipython.org/urls/github.com/barbagroup/CFDPython/blob/master/lessons/01_Step_1.ipynb) to [3](http://nbviewer.ipython.org/urls/github.com/barbagroup/CFDPython/blob/master/lessons/04_Step_3.ipynb). Using forward difference for time, backward difference for space and our 2nd-order method for the second derivatives yields:\n", "\n", "$$\\frac{u_i^{n+1}-u_i^n}{\\Delta t} + u_i^n \\frac{u_i^n - u_{i-1}^n}{\\Delta x} = \\nu \\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{\\Delta x^2}$$\n", "\n", "As before, once we have an initial condition, the only unknown is $u_i^{n+1}$. We will step in time as follows:\n", "\n", "$$u_i^{n+1} = u_i^n - u_i^n \\frac{\\Delta t}{\\Delta x} (u_i^n - u_{i-1}^n) + \\nu \\frac{\\Delta t}{\\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Initial and Boundary Conditions\n", "\n", "To examine some interesting properties of Burgers' equation, it is helpful to use different initial and boundary conditions than we've been using for previous steps. \n", "\n", "Our initial condition for this problem is going to be:\n", "\n", "\\begin{eqnarray}\n", "u &=& -\\frac{2 \\nu}{\\phi} \\frac{\\partial \\phi}{\\partial x} + 4 \\\\\\\n", "\\phi &=& \\exp \\bigg(\\frac{-x^2}{4 \\nu} \\bigg) + \\exp \\bigg(\\frac{-(x-2 \\pi)^2}{4 \\nu} \\bigg)\n", "\\end{eqnarray}\n", "\n", "This has an analytical solution, given by:\n", "\n", "\\begin{eqnarray}\n", "u &=& -\\frac{2 \\nu}{\\phi} \\frac{\\partial \\phi}{\\partial x} + 4 \\\\\\\n", "\\phi &=& \\exp \\bigg(\\frac{-(x-4t)^2}{4 \\nu (t+1)} \\bigg) + \\exp \\bigg(\\frac{-(x-4t -2 \\pi)^2}{4 \\nu(t+1)} \\bigg)\n", "\\end{eqnarray}\n", "\n", "Our boundary condition will be:\n", "\n", "$$u(0) = u(2\\pi)$$\n", "\n", "This is called a *periodic* boundary condition. Pay attention! This will cause you a bit of headache if you don't tread carefully." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Saving Time with SymPy\n", "\n", "\n", "The initial condition we're using for Burgers' Equation can be a bit of a pain to evaluate by hand. The derivative $\\frac{\\partial \\phi}{\\partial x}$ isn't too terribly difficult, but it would be easy to drop a sign or forget a factor of $x$ somewhere, so we're going to use SymPy to help us out. \n", "\n", "[SymPy](http://sympy.org/en/) is the symbolic math library for Python. It has a lot of the same symbolic math functionality as Mathematica with the added benefit that we can easily translate its results back into our Python calculations (it is also free and open source). \n", "\n", "Start by loading the SymPy library, together with our favorite library, NumPy." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy\n", "import sympy" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We're also going to tell SymPy that we want all of its output to be rendered using $\\LaTeX$. This will make our Notebook beautiful!" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from sympy import init_printing\n", "init_printing(use_latex=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start by setting up symbolic variables for the three variables in our initial condition and then type out the full equation for $\\phi$. We should get a nicely rendered version of our $\\phi$ equation." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUcAAAAfBAMAAACGzRJ2AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEIl2mSJE3e9UMqtm\nzbsXyEShAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEVklEQVRYCcVWT2gcVRj/TTY7m52dbMaQ0oKW\nxA1oLWjTf5cScbDJoauH0SgVixCk4KmSgwHRQ6PgfVGhBFSWeKgXYYIHtfGwJh5EA+4lJwmMh/ZY\ntmLxIHX9vpk3b3Z252VmFho/mP3e++b7vt9v3r43vwEOw4za7WFhNp5xhi3NWVfFWzkrZLqnr8jx\nwx2MY9oaEqGt3RuyMn/ZjfwloqI8sJJ6U9yqBk8+Euv98TxNeX8t7TSNnTqwieL8BaxrDs8mt61K\nt9uiILvJc1a5cQTkCtRK24s16p+oYE9T4nGnP3spDPxs4QqNR8M5+7L9Bv3y/rLhVpc/x5sfUUIF\nvz1wePaH5o2++oofJLdX9MzuHjvKxqPBM3OXJFPBEry2OVDgioh+08IOjWMkx/xn4v119/UZDTVg\nFdXbG/jRAs86+MDA4xwkp3XQMR8BO6xDW5wcgOoNqGDLDUy9NNObSWPNFoHHStbIO62A5NLsr0F0\n+rOT/uAGjj4ANjzmo11vY/enJs2Ka3gPeouD5IjdfXO3zo7WutT9GweYGtbF9e5KX6U+JwJzJQsr\nGN36dut7NN4eCaLTzTEe0P7aPuUB15iPWfuXtin1ueaTHKf7RJLcLf0+cNN3Va46yNSwXkKZ2YK2\nsLBwyaidOKm1KcH/u88CUxRtTMD/w2l/ucXVAiYc4lPBZaeNDs/o78YxqiGS5Mxf1sYxwY4pH2xK\nWDCHfpOPhDGrbLcCkjN3zCCv0jxK2433Fx2GY/iK+VTpOTyjzTM6OLy6TJLc6fJcFV+zo5wUU8N6\nCZWaK4LaqcXiN05A8oczF4KodrGOy/7+ulpbNmdv4cu/GvrzT2DjnMOzqW0HV8FBdmfmYcwusqPV\nTjE1bMgn1qA/GDvdscwck/XUXBUsne4Eky8scc9IyMkdslMrVLDJa6QvpzbMncCKk2IqWFacyL64\nS2ah+xCtE6HJUTbYpErZInEgRF2orf8G3USGxUps1huMN56hWxUL8aXsTT9wTKLORrso1HZScJbn\nrPaJIlE2LrqUYZNsnXDisqwoHAyzqLNRI6nt9GJMP8CylYJk1Nico1ybrpKD5OMteykGJOrle/BI\nbaW289s79VUYtVOQjBo/+SIl23QRSV4M3ybPLopRBseivmZ4rBGhtjPJVFGJWock+2CjxiXgha33\ntxo+SS8o1FfxctQiZcSibrjjDqttm3L9FxmRTJXnqK8g2Qfb0/g459p08UoyCNlrDa0ejDL9jlmF\n83XQSkptH24lB2Bl40/5U8Gmi0l65Mm+25m3glGWXxL1ylP0ueZCajuLddY9+fT+/rv7+78nwEaN\nr5yn2zZ9AD53Se7Jf7JwG8hxg0ioW/lPtwpWNCaSbOHp/jOY5vxdCvJDbbezl4s9qYIVjZtBw3AR\nPgQK2SHCzLja5lEcQVIFG28cKs6z9BkbQh+GFyTzwRYu7h4GN4khSB42rMTPNBAkM+X+b0lH8iP/\nB/BpkElltnvsAAAAAElFTkSuQmCC\n", "text/latex": [ "$$e^{- \\frac{\\left(- 4 t + x - 6.28318530717959\\right)^{2}}{4 \\nu \\left(t + 1\\right)}} + e^{- \\frac{\\left(- 4 t + x\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}$$" ], "text/plain": [ " 2 2 \n", " -(-4⋅t + x - 6.28318530717959) -(-4⋅t + x) \n", " ──────────────────────────────── ─────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) \n", "ℯ + ℯ " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x, nu, t = sympy.symbols('x nu t')\n", "phi = sympy.exp(-(x-4*t)**2/(4*nu*(t+1))) + sympy.exp(-(x-4*t-2*numpy.pi)**2/(4*nu*(t+1)))\n", "phi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's maybe a little small, but that looks right. Now to evaluate our partial derivative $\\frac{\\partial \\phi}{\\partial x}$ is a trivial task. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzUAAAA/BAMAAAArsfskAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMol2mSJE71Sr\nZruYlGYbAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAM10lEQVR4Ae1cDYxcVRX+3vztzM7O7LOkLCHV\nnSwkFlNhCkEaIfRFUkmM2imwiyKwU4VuYyAdfpqtUuQpCkig2QBCNDG7RUmKKI5CxFjNDgIt4mIH\nTAshrh2CieBP2V1Qfgpdz7k/72dmfzo7nZldmZN0zr3nnvvOefd79757v7cpUK2Eel+ttkvLv0Ej\nkMTdDYrUClPtCATQbVbbp+XfqBF4pFGBWnFmHIF7tPVlWchIFaEpY6zUbS3djBGIZlXUxE5EC1TO\nqTrpFYRPS5o3Ah069D2rkchSRWOzG8a25bqxpZsxAl0qqFFcjQtvpAphE73y6gI6EJx+oxkZtWLq\nEbhXFSJYjSCwfeycsXw8VsogqT1aulkjMA50DQ4O5tf1PmAt4yR4TYtkgQBXWtLMEdDzBtiJxxKU\nCWETTg4VW/OmmajI2Ct0Csvfsu/YLLHpONBn0fumJU0egQoI3H1akzNrhXfON3ooLFXQGGl7Szd+\nBL4/c0jmBVpS/xGYg+zvnJ5b6p/cBzzCTGR//BCJotA+4MPT1Nufl+xXGEnaM8y5Po66LWrxnU0d\njEUXfG6yn2hOFtoW3EGqjf59b5dLrHHTMZQV2yeP4dWW/qXmIfuJ5mSh7fQBUowN1gC7WddBIksT\nm5ilxkItMhnv0MT7ilRdd7kdvfoKfG6ATocmbjlgxfMncM24bAP2Tb/PRlJG3zXYbdisNlKnucl+\npjlHzSS6ED6LQmhsKk4+3mRqKC9RbG5Rt1z+LUWYH4umSJdi6TZ6wl+Lp2KX2ESwFBLTK7mWND+C\nbbfl2EgqYO/AwSmbFQ11Odlv9A3Y4pLyh2nOdnweRN+k0Tb2xNgfxbypF9m5RLEpqBEr/5YizPQK\nIEkZk8lX18EMZxC0cejbmcSHwLVOtJtW2GQjqW504jkTrOL5CrL/bjvqnZFMc34WKzEOI0UR9Lyp\nF9m5NLExcjQ0JBXfUoTx8H6LdTxtrE0hdDkYhuOmkBi/nmu3I2jhh8JIikCxx1+wWKHAvbwSfQhd\nltdANOemWy2aN/GcZ01rzRvvGMWyslbxLYXNxpsYYb3MTvQeAWIpxmZ/fwk4j2sCm5TEJoUOa9hG\nKC0Uefil42ebnvNZiOZ8Kk9vJUT32HLefPqtfN3IzqU5bxLFmb+lGPSB5QZMgFc143EatJvtLN5m\nbArRNQFat6jGa5oxIrBhNf4ZpKgHKdJl0lkoM6iqevmLNY1MrX2ad5j0vKn4liKcnjTSsND1zUwS\nbcZIaISxKSGXxKe4xnsBmktsJBUxf4VSKMWKfMqk0/uy8bTR+YYlJH7rdr7B0pw3RkENS/m3FGHu\nGrCiaaydTse2fhH7r7JCW27Afb1Dob5tXDNu/AJiObCRVIT22+sGbFaV75t2wiarIvmUj/asGy+w\nNLGpGMacb+gWVqF9WplEUuiyy2wNrC5RbPT5Ro+UpQs16LbKvgNX03SqSu48eu8X53Fte2bq2dlc\notXmNduFarQn7MoLxIYqbbVaNtZ6AeofLTkXMTZdPEOS0cFDlnLpsB3fqgs73B4zx7n1G39zXcpK\nOofAxtAN9DcPfRdRe++VFr7T+0/ay+7t7c0YW/kYIltk5/Mhm0UfY+ulNiBCPC2b6Vfz9XN/Sqmt\ndcKJtoBCwHY60TQk2uA41yBaViD8jnIxUqqwALXS7TNjHKOA/hmeDBrBnXwEkDm0T09nYJyCa+m8\nMdSWiuYwWsRBGr1CIG+sUS34jcmxzoVsFn1OLEZ3QYYQf3HkJnN0paS4JGKWcK83na9zekEXgN8B\n3wJGlSGelYXHgZ9onz5dqFoncm4XfxzabrIQWdg+Ikr+H8Ft6xwCl24GkjkcAv6OZClyGJ1pbGMq\n9yXgAtli7N3HAxl5GLJZ9HkG2KtChLP+AEdVe0liQy+mO8ifHq560vlORmucEoiPOA78TArR2PwI\n2Gcp24lKV6+CttvHH0dh055Gx5uuj1vi/YXOQRzT6OhNJ/K3+edBdKeQAS7Cv+ihMkULsIUH8s5T\nVbPoMwX0Q4bgI2K1EjtPYlNoDJ2v0qOtuyPDL+My4ExV19icRuy3pWyBIce5ysLNHn9/HIVNYHIO\nbHQOYpz/wZdKjsgL0ppGGGXxHjCcFy0Km8yp7EDN3McgtmWLpUKcwg3VyYeDZnwSJSLeGkLnq+Qi\nOTfLjun9xdifjuyRFo0N1daaxqan76UFLezxdvsdTWnA4+SPo7Ch9uAkjr/qu1dmPa5UVPvytSaV\nO8bp+8nUuoEhdJ7c9zUyGOJJSsJ4i7DJihaJTcJmbLhZ9Lmf5s0QGSgELiaNMr6eTXNINmhiJFSi\nZ6AhdL7KJFn0pNQ/ZSI+oQwuNrH/YDn+nD0dSKQ83lUV/+319sVxsRnOGJnEJ4JlMSQ2lANJ0oy+\nY0wV8RC6P472PKIXXsPmyxCl5vUbRIvE5ngQNqKZ+9CChtP4ToczwPPcpYyvZ9PsEuq95KJQIWAn\nijBS5MbvG/5EGWBdPwnm3WvH/7D2QbSllMHFJljABvSbv2ai3PWurnS2x90fx8XmHITN8GTc9rhS\nUWJDOUh5xZg28XWzexJttIHDJUM0ByYQpXmzfoNokdhkGBvZDLyCQDbez3d6Dv3jj/iVfD0Z55B2\nM7L5el47q6fzjZPOIDndQnWbcUqmne5tGXc+42T8BZF3zWCGUwz39Jx0f09PgcsYpK9I8lZDE45z\nFaEm6BoP0r8Z4lzQ0/PJnp5VHAVtORgIlERZu1JFYkM5SLnJfgM4OET7swRv7dsJoEhKrWmiRWAT\nt2TC3AzcZOPAFf2WCEEYkqWCr2e32cXo39bxJXoICo2h81UijI2WXZR3vlsbnHkjX0n/ZS/CZoHC\n2Gjxx3HmjXgldWa1l9YCG/Va/D2/NmhtPViklS/xftRG8jC9RErgvcBoXrQIbJaDsJHNog9dbJ8J\niBCMTWeBfqqXguzSpnrurv4K1fTwrGm0NKCjOEwPrxAHm8f4JkNifxtLV3Ntr+/ZbqUsjsYmnsNt\n9HBYrqMsCWxEDsD9NAfMZ3jedIzQvOmcFNh0F8AnnoOyRWDz27Gxd/fIZtGHLrWKlmQR4m6qzMbX\ny5Cz/tL5hiUkFXJK10l59wL0PCftC/CoDKWxCRH1bf8iMsF/OZJILTQN717AH0djcz6wOW5uoW8g\n/hiMjciBzDQYp2M9vW8QovfNSHsJQTrnrCcznT2vlS3yfQM8DNks+rxkh8iRQ8i9QHtmFr7eH7u8\n5ifeyjMt9661HqbUtdxl0h9KbTeysq6xubn34idih5MT0SKt13RLCxOxcVVd/XEUNtEze7enh/Ov\n486yAIwN58DHzFuR2IWOjPEQ8BROzNO+aTRLixndQzBv/Fi1yLMnjkA2iz5fNpcVIUIAWylApLl8\nvbzF8jtVN56wVSGWVgVSob3Ede4YVwaNzWnT028YezYOsD1g0c/MMk8gXki0+OMobAK0jUmv2HT8\njfQIeEVw25wDaKsW3brXpAPKdot2K4NPAvcNvka+t+f5vHIVGUXL9p+vomtcPv2sbBZ9woN0GhIh\ngJP58tXz9dyrZqEFwJFoySn6C0/rKm3TZxaNjbfVe7onO60RWuYLFLS1Z7nWa1q5vV51LxNSrxiz\nXpe3LVoCdiWxLHjaZdrjBV0o1wZ1LZc+v+FctzpfoITngXF7cekEf7XutXC27iFmD/DVktvGI6+J\nZfWASp7WyTBpu97zlYyUz4OoXkfmDfSq49rkgvNUNiGP6PMlNyqvWNeqqrN4bDFpfR5R1lkXI9Xu\nVWU4EtXryLyBdjiuTS408yEJq6M1D4FYW89Ug+HDBqfoIZrlJa6bvfpFbwWQVK+wzR+It3mLQZxd\nUDOSeZGwuesB4K8UnI7TLrHsx8a7qV1Ymkz1xj6WQaJU50ALS28R9jKyhI1BTAbBAz5ZOsSyH5vn\na85dUL3defA5rq6Bas50sVwgzHQhfQiKv00ZBfN0Xkyp1PzYPFJzwoLqHTUxatc5UM2ZLpYLPMrY\ntKURyFFGzGRKYtlL+PJeQLCxNeUcF1TvdcBH6xyopiwXU2cjw9gEc/gKA8DYOMSyf94wG1uTCKoX\nZyGymi5Tz0A1ZbmYOsfHxvb9NNudj4kPF7ymOcSyHxsvibKgGxBUL73YfsAPQT0DLSi7RdqpvUR/\ncvJLm7PjV7RDLPuxqX0vwFRvfDK+sgGBOMT/hXSWcN3yjLiVcI7+2zRNLPuxYTa2VjmCttW94iJ1\nDlRrooulf/j194qrnpPZxNJwiWWNjeRpBRtbW85E9QaeNcU16huotjQXa2+mUrRobET9WLOxDQuk\n72bpay/LfIL3dhyu02usodywQDXkuMi6lrGTbnbHmo1tWCD3FpZ6aVaW+VizsQ0LtNQR8eQ/C8t8\n7NnYhgXy3NxSKP4PMw1WTFYAYHoAAAAASUVORK5CYII=\n", "text/latex": [ "$$- \\frac{e^{- \\frac{\\left(- 4 t + x\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}}{4 \\nu \\left(t + 1\\right)} \\left(- 8 t + 2 x\\right) - \\frac{1}{4 \\nu \\left(t + 1\\right)} \\left(- 8 t + 2 x - 12.5663706143592\\right) e^{- \\frac{\\left(- 4 t + x - 6.28318530717959\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}$$" ], "text/plain": [ " 2 \n", " -(-4⋅t + x) -(-4⋅t + x - \n", " ───────────── ─────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν\n", " (-8⋅t + 2⋅x)⋅ℯ (-8⋅t + 2⋅x - 12.5663706143592)⋅ℯ \n", "- ─────────────────────────── - ──────────────────────────────────────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) \n", "\n", " 2 \n", "6.28318530717959) \n", "───────────────────\n", "⋅(t + 1) \n", " \n", "───────────────────\n", " " ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phiprime = phi.diff(x)\n", "phiprime" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want to see the unrendered version, just use the Python print command." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 12.5663706143592)*exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))\n" ] } ], "source": [ "print(phiprime)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Now what?\n", "\n", "\n", "Now that we have the Pythonic version of our derivative, we can finish writing out the full initial condition equation and then translate it into a usable Python expression. For this, we'll use the *lambdify* function, which takes a SymPy symbolic equation and turns it into a callable function. " ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-2*nu*(-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 12.5663706143592)*exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)))/(exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1))) + exp(-(-4*t + x)**2/(4*nu*(t + 1)))) + 4\n" ] } ], "source": [ "from sympy.utilities.lambdify import lambdify\n", "\n", "u = -2*nu*(phiprime/phi)+4\n", "print(u)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Lambdify\n", "\n", "To lambdify this expression into a useable function, we tell lambdify which variables to request and the function we want to plug them in to." ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.4917066420644494\n" ] } ], "source": [ "ufunc = lambdify((t, x, nu), u)\n", "print(ufunc(1,4,3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "###Back to Burgers' Equation\n", "\n", "Now that we have the initial conditions set up, we can proceed and finish setting up the problem. We can generate the plot of the initial condition using our lambdify-ed function." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "array([ 4. , 4.06283185, 4.12566371, 4.18849556, 4.25132741,\n", " 4.31415927, 4.37699112, 4.43982297, 4.50265482, 4.56548668,\n", " 4.62831853, 4.69115038, 4.75398224, 4.81681409, 4.87964594,\n", " 4.9424778 , 5.00530965, 5.0681415 , 5.13097336, 5.19380521,\n", " 5.25663706, 5.31946891, 5.38230077, 5.44513262, 5.50796447,\n", " 5.57079633, 5.63362818, 5.69646003, 5.75929189, 5.82212374,\n", " 5.88495559, 5.94778745, 6.0106193 , 6.07345115, 6.136283 ,\n", " 6.19911486, 6.26194671, 6.32477856, 6.38761042, 6.45044227,\n", " 6.51327412, 6.57610598, 6.63893783, 6.70176967, 6.76460125,\n", " 6.82742866, 6.89018589, 6.95176632, 6.99367964, 6.72527549,\n", " 4. , 1.27472451, 1.00632036, 1.04823368, 1.10981411,\n", " 1.17257134, 1.23539875, 1.29823033, 1.36106217, 1.42389402,\n", " 1.48672588, 1.54955773, 1.61238958, 1.67522144, 1.73805329,\n", " 1.80088514, 1.863717 , 1.92654885, 1.9893807 , 2.05221255,\n", " 2.11504441, 2.17787626, 2.24070811, 2.30353997, 2.36637182,\n", " 2.42920367, 2.49203553, 2.55486738, 2.61769923, 2.68053109,\n", " 2.74336294, 2.80619479, 2.86902664, 2.9318585 , 2.99469035,\n", " 3.0575222 , 3.12035406, 3.18318591, 3.24601776, 3.30884962,\n", " 3.37168147, 3.43451332, 3.49734518, 3.56017703, 3.62300888,\n", " 3.68584073, 3.74867259, 3.81150444, 3.87433629, 3.93716815, 4. ])" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from matplotlib import pyplot\n", "%matplotlib inline\n", "\n", "###variable declarations\n", "nx = 101\n", "nt = 100\n", "dx = 2*numpy.pi/(nx-1)\n", "nu = .07\n", "dt = dx*nu\n", "\n", "x = numpy.linspace(0, 2*numpy.pi, nx)\n", "#u = numpy.empty(nx)\n", "un = numpy.empty(nx)\n", "t = 0\n", "\n", "u = numpy.asarray([ufunc(t, x0, nu) for x0 in x])\n", "u" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoMAAAGnCAYAAADMu10zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1w3WWd9/HPt0maUAoJWBoeUmlpaUuAUXwGlxC79BRF\nve+Z+97pOqKMEGV3oSnuLCCESAqN6w4rGHD8Y+9VZn241dGddXaNQw9i08iTAoVFaRUpTy03LaiY\ntClJk3jdf5wE03OuX3LOyTnnd8653q+ZTvPwy69XOSX99Hv9vtfXnHMCAABAmBbEvQAAAADEhzAI\nAAAQMMIgAABAwAiDAAAAASMMAgAABIwwCAAAELBZw6CZfd3MDpjZL2d87EQzu8/MnjGzpJk1FX+Z\nAAAAKIa5KoP3SLok7WOfk3Sfc261pPun3gcAAEAFsrkOnTaz5ZL+yzl37tT7v5Z0kXPugJmdLGnA\nObe22AsFAABA4eXzzGCzc+7A1NsHJDUXcD0AAAAoodr5fLFzzpmZt7QY9XEAAACUnnPOfB/PpzI4\nvT0sMztF0quz/KL8iPnHLbfcEvsaQv/Ba1AeP3gdyuMHr0N5/OB1iP9HqV+D2eQTBv9T0uVTb18u\n6Yd53AMAAABlYK6jZb4j6SFJa8xsr5l9StIXJa03s2ckrZt6HwAAABVo1mcGnXMfi/jUxUVYC4qg\nvb097iUEj9egPPA6lAdeh/LA6xC/cnoN5jxaJu8bm7li3RsAAADZMzO5AjaQAAAAoEoQBgEAAAJG\nGAQAAAgYYRAAACBghEEAAICAEQYBAAACRhgEAAAIGGEQAAAgYIRBAACAgBEGAQAAAkYYBAAACBhh\nEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYBAAACRhgEAAAIGGEQAAAgYIRB\nAACAgBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYB\nAAACRhgEAAAIGGEQAAAgYIRBAACAgBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQA\nAAgYYRAAACBghEEAAICAEQYBAAACRhgEAAAIGGEQAAAgYIRBAACAgBEGAQAAAkYYBAAACBhhEAAA\nIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYBAAACRhgEAAAIGGEQAAAgYIRBAACA\ngBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYBAAAC\nRhgEAAAIGGEQAAAgYHmHQTO70cyeNrNfmtn/NbP6Qi4MAAAAxZdXGDSz5ZI+LekdzrlzJdVI+uvC\nLQsAAAClUJvn1w1LGpe0yMwmJS2S9HLBVgUAAICSyKsy6Jz7g6QvSXpJ0v+T9Efn3E8KuTAAAAAU\nX16VQTNbKelaScslDUn6vpl93Dn37ZnX9fT0vPl2e3u72tvb810nAAAAsjQwMKCBgYGsrjXnXM6/\ngJltlLTeOdcx9f4nJL3POXf1jGtcPvcGAABAYZmZnHPm+1y+3cS/lvQ+MzvGzEzSxZJ25btAAAAA\nxCPfZwb/W9I3JD0m6ampD/9LoRYFAACA0shrmzirG7NNDAAAUBaKsU0MAACAKkAYBAAACBhhEAAA\nIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYBAAACRhgEAAAIGGEQAAAgYIRBAACA\ngBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYBAAAC\nRhgEAAAIGGEQAAAgYIRBAACAgBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgY\nYRAAACBghEEAAICA1ca9AAAohf7+Qd11V1JjY7Wqr59QZ2dCl17aFveyACB2hEEAVcUX+iTpmmu2\n6YUXet+8bs+eLkkiEAIInjnninNjM1esewOAT3//oDZv3qY9e/4c+hobuzQ29rpGR7+acf2GDd26\n997bSrlEAIiFmck5Z77P8cwggKpxxx3Jo4KgJA0N9Wp0dMR7/ehoTSmWBQBljW1iABUnfSt4/fqE\ndu9u0/bt/m9pixeP6dChzI83NEwWeaUAUP4IgwAqim8rOJnsmnprwvs1Z565WMPDXUd9TVPTTdq0\n6ZJiLhUAKgLPDAIoW+kVwE2bErr11qQefXRrxrUrVnTrhhvW6/bbjw6KK1fepL6+VOi7++779Npr\nNdq5c1I1Nev19NNtWrOmZL8dAIjNbM8MEgYBlCVfBbCurkvj4yOSvpxx/UUX9WhgoEf9/YO6++77\nNDpao4aGSW3atD6jY7ijQ/ra16QPflD68Y+L/TsBgPgRBgFUnA0bblYymVkBNNso577nuT77zuBX\nX5XOPFMaHpZ+9CPp0kvnvVwAKGuzhUGeGQQQu5nbwWYTWrkyoZ/9zP/tqbW1SaOjXRlbwbk8/7d0\nqbRli/TZz0rXXitdfLFUXz/v3wYAVCQqgwBi5dsOlrokvS7Jfzbgpk3r59wKnsv4uPS2t0m7dw/q\nzDOTOvVUJpMAqF5sEwMoC+kNIZdfntCWLUk980zmdvAZZ3TIrNnbDFKosHbbbYP6/Oe3SZr5a3Sp\nr28DgRBAVWGbGEDsoo+E8R8IvWxZi667bp3uvrt7RgWwcEFQkh54IKmZQVCS9uzp1d13dxMGAQSD\nMAig4Hzzgf/5nzOng0i9qqvbqPHxzHs0NEzq0kvbihrKxsb83wKZTAIgJIRBAAXlqwA+9FCXDh3y\nVwBXr55/Q0i+6uv9h1QzmQRASAiDAArqrrsyK4CHDvVK2ui9vqVl6VRDSPG2g6N0dib07LNdeu65\n0gdRACgXhEEAeUvfDm5vT+iJJ/zfVlavbtLkpL8CWOzt4CjTv+ZHPtIt52p08cWTuvba0gRRACgX\ndBMDyEtcR8IUw0knSb/7nbR/v9TcHPdqAKDwOFoGwLz4ZgRv2ZLUY49lHglz8skdqq9v1osvFu9I\nmEJbtUras0f6zW+k1avjXg0AFB5HywDIm68CuH379IzgTGvWFP9ImEJrbEz9PDQU7zoAIA6EQQCz\n8jWEjI/3Ts0Izry+FEfCFBphEEDICIMAJPmng7zySpsGB4s3I7hcEAYBhIwwCGCW6SCS5D+LL84j\nYQrt+ONTPw8Px7sOAIgDYRAITC7TQZYs6dbf/E1C3/lOeR0JU2hUBgGEjDAIBMRXAXzwwS6NjPib\nQc4+u0a33dam971PVVEBjEIYBBAywiAQEF8zyMhI9HSQ6bFs1VIBjEIYBBAywiBQpQo1HSQEPDMI\nIGSEQaAKRTeEvO69fsWK6mkGyQeVQQAhIwwCFS5qOoivISQ1HaQrYzpINTWD5IMwCCBkhEGggvkq\ngD/9aZcmJqpnOkgpEAYBhIwwCFQwX0PIxER1TQcpBZ4ZBBAywiBQIWZuB5tNaMWKRBDTQUqByiCA\nkBEGgQrg2w4eGIhuCKmm6SClQBgEEDJzvr2kQtzYzBXr3kA1S28I+eQnE7r11qSeeWZrxrVnnNEh\ns+aMCmBfH8EvF85JNTWpn48ckerq4l4RABSWmck5Z77PURkEykj0kTD+hpBly2gIKQSz1HODQ0PS\nwYPSiSfGvSIAKB3CIFBGomYE19Vt1Ph45vU0hBROY2MqDA4NEQYBhIUwCMRkejt4dLRWhw9PqKEh\noQceiJ4QQkNIcfHcIIBQ5R0GzaxJ0r9KOluSk3SFc+6RQi0MqGb9/YO65ppteuGFmVVAGkLiRBgE\nEKr5VAb7JP3YOfe/zaxW0rEFWhNQNdKbQa65JqFjj23TlVcmdeBA5nbw2Wd3RFYA2Q4urukwyFmD\nAEKTVxg0s0ZJFzrnLpck59yEJP49DcwQPR1Eivpfb8kSGkLiMn3wNJVBAKHJtzK4QtJrZnaPpLdJ\nelzSZufc4YKtDKgg6RXAzs5E5HSQ+vputbQ47dmTeR8aQuLDNjGAUOUbBmslvUPSNc65R83sy5I+\nJ+nzMy/q6el58+329na1t7fn+csB5ctXAXz88S4ND/uPg3nve2t0/fXrtHkzDSHlhDAIoJoMDAxo\nYGAgq2vzDYP7JO1zzj069f4PlAqDR5kZBoFq5asA/v73vZI2eq8/5pjJNyt/bAeXD54ZBFBN0otw\nW7Zsibw2rzDonNtvZnvNbLVz7hlJF0t6Op97AZXENx3kt7/1/290+ulNqq2Nrv6xHVxeeGYQQKjm\n0028SdK3zWyhpD2SPlWYJQHlKXo6iP84mLVrOQ6mkrBNDCBUeYdB59x/S3p3AdcClA1fQ0jUdJDj\nj+/QokVd2r+f42AqGWEQQKiYQAKk8VUAH3ywSyMj/oaQ887jOJhqwDODAEJFGATS+BpCRkaiG0I4\nDqY68MwggFARBhG09O3gtraEdu70/2+xZk2TJiY4DqZasU0MIFSEQQQr14aQ5ctpCKlmhEEAoTLn\nXHFubOaKdW8gV74ZwVu2JPX441szrj355A7V1zfrxRePrgD29RH8qtnEhFRXJ5ml3l6wIO4VAUDh\nmJmcc+b7HJVBVL3oGcH+hpA1a2gICVFtrbRokXT4sDQyIh13XNwrAoDSIAyi6kXNCDbbKF/xmoaQ\ncDU2psLg0BBhEEA4CIOoKjO3g6UJLV+e0OCg/495a2uTRkdpCMGfNTZKr7ySCoMtLXGvBgBKgzCI\nquHbDt6xI7ohpKWFhhAcjSYSACEiDKIipTeEfOITCd16q39CyBlndMjMXwFkOxgzcfA0gBARBlFx\noo+E8TeELFtGQwiyw8HTAEJEGETFuf12fwWwrm6jxsczr6chBNlimxhAiAiDKFvpW8EXXpjQU0+1\naccO/x/b1atpCMH8EAYBhIgwiLIUvRUsSRPer6EhBPPFM4MAQkQYROyipoP4toJXrerWzTcndNtt\nNISg8HhmEECICIOIVa7TQU47rUaXX96mJUtEBRAFxzYxgBARBhGrfKaDSKICiKIgDAIIEWEQJeOb\nDhLVDMJ0EMSBZwYBhIgwiJJgOggqAc8MAgiROd9eXCFubOaKdW+Ut6jpIL/97daMa1PTQZozKoB9\nfQQ/lN4zz0hr1kgrV0rPPhv3agCgcMxMzjnzfY7KIAqK6SCoZDwzCCBEhEHkLb0C2NmZYDoIKhrP\nDAIIEWEQefFVAB98sEsjI/4KINNBUAnq66W6OunIEWl0VGpoiHtFAFB8hEHkxXckzMhIr6SN3utp\nCEElMEtVB3/3u9RWMWEQQAgIg5hT+nbwX/xFQjt3+v/orFnTpIkJpoOgcs0Mg83Nca8GAIqPMIhZ\nRTeE+I+EWb6cCiAqG88NAggNYRBvymVG8Mknd6i+vksvvkgFENWFswYBhIYwCEm5zwhes4YjYVCd\nOF4GQGgIg5CU34xgKoCoRoRBAKEhDAYmfSv4sssS2revjRnBwBSeGQQQGsJgQKKbQSRpwvs1HAmD\n0FAZBBAawmCVymU6yNKl3fq7v0vom9/kSBiABhIAoSEMVqFcp4OcdVaNbrmlTe96l6gAInhUBgGE\nhjBYhXKdDtLQMClJVAAB8cwggPAQBitc+nbw+9+f0OOP5z4dBEAKlUEAoSEMVjCmgwCFxzODAEJj\nzneIXCFubOaKde8Q+aaD9PQktXPn1oxrTzmlQwsXNmdMB+nrI/gBc3n6aemcc6S1a6Xdu+NeDQAU\nhpnJOWe+z1EZrAC5TgdZvZrpIEC+2CYGEBrCYAVgOghQOjSQAAgNYbDMzNwOlib01rcmmA4ClNCx\nx0pm0siINDEh1fJdEkCV49tcGfFtB0vRDSFMBwEKb8GCVBPJ0FCqOnjiiXGvCACKiwaSmPhmBN96\na1LPPpvZELJyZYek5owKIA0hQHGcfrr00kvSc89JK1bEvRoAmD8aSMpM9JEw/oaQlhYaQoBS4rlB\nACEhDMYgakZwXd1GjY9nXk9DCFBanDUIICSEwSKb3g4eHa3VoUMTWrgwoUce8f9nX72ahhCgHHC8\nDICQEAaLqL9/UNdcs00vvEBDCFBJCIMAQkIYLADfdJCFC9t0xRVJvfpq5nbwOed06I03/BVAtoOB\n+PHMIICQEAbnKXo6iBT1n/ctb6EhBChnVAYBhIQwmIP0CmBnZyJyOkhDQ7daWpyefTbzPjSEAOWN\nBhIAISEMZslXAXzssS4ND/uPg3nPe2p0/fXrtHkzDSFApaEyCCAkhMEs+SqAf/hDr6SN3uuPOWby\nzcof28FAZeGZQQAhIQx6pG8Hf/zjCT3zjP8/1emnN6m2Nrr6x3YwUHmoDAIICWEwTfR0EP9xMGvX\nchwMUG14ZhBASIKeTexrCLn99qR27MicD9zY2KFFi5r1yivMBwaq3ZNPSuedJ517rvTUU3GvBgDm\nj9nEHr4K4IMPdmlkxN8Q8va3cxwMEAqeGQQQkmDDoK8hZGQkuiGE42CAcPDMIICQBBEGZ24HL1w4\noQsuSOjxx/2/9bVrmzQ+znEwQMiOOy718/Cw5Jxk3o0VAKgOVR8GfdvB990X3RBy+uk0hAChq6uT\nFi2SDh+WDh36czgEgGpUVQ0k6Q0hV1+d0JYtSe3cmdkQcsopHaqvb9YLL9AQAiDTqadKr7wi7dsn\nnXZa3KsBgPkJooEkekawvyFk9WoaQgBEa2xMhcGhIcIggOpWNWGwr88/I3jBgo36058yr6chBMBs\nOGsQQCgqMgzO3A6WJtTSktCOHf7fyllnNWl0lIYQALmhoxhAKCouDPq2g6XohpCWFhpCAOSOswYB\nhKJsw6BvOsh557Xp2mszt4OlXq1a1SHn/BVAtoMB5IrKIIBQlGUYjJoOcviw5Jx/yaedRkMIgMIh\nDAIIRVmGwejpIN1qbnY6cCDza2gIAVBINJAACEXsYTCX6SDnn1+jrq512ryZhhAAxUVlEEAoYg2D\nuU4HOf74yTcrf2wHAygmGkgAhKJkYdA3HaSnx98MkpoO0pUxHWS6+sd2MIBiozIIIBQlCYO+CuD9\n93dpcpLpIADKE88MAgjFvMKgmdVIekzSPufcR6Ku800HmZxkOgiA8kVlEEAoFszz6zdL2iXJ+T55\n0UU36+MfH5x1OsjKlV1HfSy1Hbx+nssCgPnZuXNQ0s3atatHGzbcrP7+wbiXBAA56+8f1IYNN896\nTd6VQTNrkfQhSb2S/t53zeDgVjEdBECl6e8f1Be/uE1Sr0ZHpWRS2rMn9Q9Xvj8BqBT9/YPq7Nym\n557rVSqu+Zlz3qLenMzs+5K+IOl4Sf+Qvk1sZm66YJiaDtKccRxMXx/BD0D52bDhZiWTWz0f79a9\n994Ww4oAYHbpjbof/WhC//RPSe3dO/29zOScM9/X5lUZNLMPS3rVOfeEmbVHX9kjSRoff1ZXXXWu\nduygAgig/I2N+b81jo7WlHglADA3X6NuMnmZpGeVTdTLd5v4AkkfNbMPSWqQdLyZfcM598mjL+uR\nJK1d260bb9ysG2/M81cDgBKqr5/wfryhYbLEKwGAo6VXADs7E/rCF3xH9X1L9fUbNTbWM/X+lsh7\n5hUGnXM3SbpJkszsIqW2iT/pu5bpIAAqTWdnQnv2MOkIQHnxVQB37OjS2Jj/qL5Vq5o0OtrlCYpH\nK9Q5g94HDzds6GY7GEDFmf6e9bGPdevgwRqdf/6kurr4XgYgXnfdlVkBHBvrlbTRe/3MRt1t26Lv\nm3cDyVzMzBXr3gBQCh/4gDQwIP3kJ9Jf/mXcqwEQivSt4M98JqGDB9vU2dmjgwd7Mq5vbb1KY2NL\nZm3UNStwAwkAhID5xABKzbcVfN99XUrV1/zPMy9bNr+j+giDABCBKSQAisnXDHLnnZlbwc71qrGx\nW5dfntCPftQ1dW5gyvTzzPOZ3EYYBIAIzCcGUCy+CuAjj3Tp4EF/M8jb316jvr42JRIq+LAOwiAA\nRKAyCKBYfM0gw8PRzSDTR1vNpwIYhTAIABF4ZhBAIfimg/zqV/4ItnJlk6TSHm1FGASACFQGAcyX\nfzpIl6TXvdevWjW/ZpB8EAYBIALPDALIRfbTQXq1ZEmHGhq6tG9fYZtB8kEYBIAIVAYBZCvX6SBn\nn92i665bV9IKYBTCIABEIAwCyFau00EaGiZLXgGMQhgEgAg0kADwmbkdXFs7oXPOSejhh/2RqrW1\nSWNj5T3rnDAIABGoDAJI59sOvv/+6IaQ+U4HKQVmEwNAhOHhVCA89ljp0KG4VwOg1Hwzgm+5Jamn\nn96ace3pp3eopqY5YzrIzPnAcWI2MQDkYfFiyUwaGZEmJqRavmMCwYieEexvCFm+vHwaQnLFtzYA\niLBgQep4maEh6eBB6YQT4l4RgFKJmhFcU7NRk5OZ15dTQ0iuCIMAMIvGxlQYHBoiDALVauZ28JEj\nEzrhhIS2b/dHpLVrmzQ6Wt4NIbkiDALALDh4Gqhu/f2D6uzcdtSzflJ0Q0hLS/k3hOSKMAgAs6Cj\nGKgOvukgy5a16W//Nqm9ezMnhKxd26HxcX8FsFK3g6MQBgFgFpw1CFS+6OkgUlQUam6u3IaQXBEG\nAWAWVAaByhc1HaSmplunnuq0d2/m11RyQ0iuCIMAMAueGQQqSy7TQc4/v0af+9w6bd5cXQ0huSIM\nAsAsqAwClSPX6SDHHjv5ZuUvhO3gKIRBAJgFzwwC5Sm9IeTTn06opydzO1jqnZoO0pUxHWS6+hfK\ndnAUwiAAzILKIFB+QpoOUgqEQQCYBc8MAvHyHQkT0nSQUiAMAsAsqAwC8fFVAB95pEsHD/orgNU4\nHaQUCIMAMAueGQTi4zsSZni4V9JG7/XVOB2kFAiDADALKoNAaaRvB3/4wwn98pf+mLJyZZOkMKaD\nlAJhEABmwTODQPH5toOTyegjYVatogJYSOacK86NzVyx7g0ApfL730tLlkhNTdLr/r+XAOTA1xDS\n25vUww9vzbj2pJM61NDQfNTs4JUrb1JfH8EvV2Ym55z5PkdlEABmMV0ZHB6WnJPM+60UQDaiZwT7\nG0JaWzkSphQIgwAwi7o6adEi6fBhaWREWrw47hUBlStqRnBUQwhHwpQGYRAA5nD88akwODREGASy\nlT4juLU1oYce8seO1tYmjY1xJExcCIMAMIfGRmn//lQYPO20uFcDlL9cZwQvW0ZDSJxoIAGAObz3\nvdIvfiE99JB0/vlxrwYoL1Ezgp9+OrMhJDUjuDljRjANIcVHAwkAzAMHTwN+zAiuDoRBAJgDZw0C\nflEzgmtrN2piIvN6GkLKE2EQAObAFBKEzjcdZNeuNm3f7o8Ra9YwI7iSEAYBYA6EQYQsejqIJHnK\nf2JGcKUhDALAHHhmEKGImg6SvhUs9eqtb+3W9dcndOedzAiudIRBAJgDlUGEINfpICtW1Ojqq9u0\nfLmoAFY4wiAAzIEGEoQgn+kgkqgAVgHCIADMgcogqg3TQTATYRAA5sAzg6gmTAdBOiaQAMAcHn1U\nes97pHe+U3rssbhXA2TPNx3klluS2rWL6SChYQIJAMwDzwyiEjEdBNkiDALAHHhmEJWI6SDIFmEQ\nAOZAGES5m7kdfOTIhJqaEkwHQdYIgwAwh4YGqa5OOnJEGhuT6uvjXhHwZ/39g9q0aZuef35mFTC6\nIYTpIEhHAwkAZGHJEun3v5cOHJCWLo17NQhVekPIpZcmdPvtSe3bl9kQctZZHTpypDmjAkhDSJho\nIAGAeWpsTIXBoSHCIOIRPSPY3xCydCkNIcgOYRAAssBzgygl34zgrVv9M4IbGjZqdDTzHjSEIFuE\nQQDIAgdPo1RynRG8ciUNIZgfwiAAZIGzBlEquc4IpiEE80UYBIAssE2MQkvfCu7oSOiPf2zLa0Yw\n28GYD8IgAGSBMIhCip4OIkmeE6HFjGAUD2EQALLAM4PIl68Z5I47/NNBmpq6dcUVCf3wh10ZM4Kp\nAKJYCIMAkAWeGUQ+fBXARx7p0sGD/maQt72tRl/6UpvWrRMVQJQMYRAAssA2MfLhawYZHo5uBmlo\nmJQkKoAoKcIgAGSBMIi5pG8Hf/CDCT31lP+v2VWrmuQcx8GgPBAGASALPDOI2URPB/HPB165kmYQ\nlA/CIABkgWcGMS2X6SAnndShhoYu7d1LMwjKF2EQALLANjEkfwVwYKBLR474G0JaW5kPjPJHGASA\nLBAGIfkbQo4cmb0hhAogyh1hEACyQBgMz8zt4JqaCZ11VkIPPpj7dBCg3BEGASALixdLZtLIiDQ5\nKdXUxL0iFJNvO/inP41uCGE6CCqZudTsm8Lf2MwV694AEIfGxlQ38R/+IJ1wQtyrQaH4ZgT39CS1\na9fWjGtPP71DNTXNGdNB+voIfihvZibnnPk+R2UQALI0HQaHhgiD1SJ6RrC/IWT5chpCUH0IgwCQ\npcZGae9enhusJlEzgmtrN2piIvN6GkJQjQiDAJCl6bMGOXi6Ms3cDh4bm1BjY0IDA/6/BtesadLo\nKA0hCENeYdDMlkn6hqSlkpykf3HO3VXIhQFAuaGjuHL19w9q06Ztev75mVXA6IaQlhYaQhCOfCuD\n45I+65x70swWS3rczO5zzu0u4NoAoKwQBsufbzrIySe36aqrknr55cwJIWed1aEjR/wVQLaDEYq8\nwqBzbr+k/VNvHzKz3ZJOlUQYBFC1CIPlLXo6iBT1193SpTSEAPN+ZtDMlks6T9LP53svAChn02GQ\nZwbLU9R0kNrabp16qtNLL2V+DQ0hwDzD4NQW8Q8kbXbOHUr/fE9Pz5tvt7e3q729fT6/HADEarqB\nhMpg/NKng6xdGz0d5Pzza3TDDeu0eTMNIQjHwMCABgYGsro27zBoZnWS/l3St5xzP/RdMzMMAkCl\nY5u4POQ6HWTRosk3K39sByMU6UW4LVu2RF6bbzexSfqapF3OuS/ncw8AqDSEwdJLbwi58sqEtmzJ\n3A6Weqemg3RlTAeZrv6xHQz45VsZfL+kyyQ9ZWZPTH3sRufcvYVZFgCUH54ZLC2mgwClkW838QOS\nFhR4LQBQ1nhmsHh8R8IwHQQoDSaQAECW2CYuDl8F8OGHu3TokL8CyHQQoLAIgwCQJcJgcfiOhDl4\nsFfSRu+cmw2HAAAJXUlEQVT1TAcBCoswCABZ4pnB+UvfDr7kkoSeesr/V9GZZzbpT39iOghQbIRB\nAMjS9DODw8OSc5JZvOupNL7t4GQy+kiYM86gAgiUgjnninNjM1esewNAXBYtkt54Qzp4UFq8OO7V\nlK/0CuCmTQlt3ZrUz3++NePak07qUENDs/buPboC2NdH8AMKxczknPP+E5bKIADkoLExFQaHhgiD\nUaJnBPsbQlpbORIGiBNhEABy0Ngo7d+fCoOnnRb3aspT1IzgqIYQjoQB4kUYBIAczHxuEJkzgtes\niZ4R3NrapLExjoQByg1hEABywPEyf5brjOBly2gIAcoRYRAAchBiGPRNB7nwwjZdd13uM4LZDgbK\nD2EQAHIQWhiMmg4yNiYdOeL/K4QZwUBlIQwCQA5Ce2YwejpIt044wel1z44wDSFAZSEMAkAOqrky\nmMt0kHe/u0a33LJOmzfTEAJUOsIgAOSgWsNgrtNBTjxx8s3KH9vBQGUjDAJADqohDEZNB/E1g5x0\nUoeOOaZLL73kr/6xHQxUPsIgAOSg0p8ZZDoIgHSEQQDIQaVXBpkOAiAdYRAAclBJYZDpIACyQRgE\ngBxUShhkOgiAbJlzrjg3NnPFujcAxOWVV6RTT5WWLpUOHIh7NSnpDSFXXplQT09Su3dvzbg2NR2k\nOWM6SF8fwQ+oZmYm55z5PkdlEABy8NBDg5KSeu21Wm3YkBrNFmeIij4Sxt8QwnQQAOkIgwCQpf7+\nQd1wwzZJvXJOSialPXu6JCm2MHXHHf4jYWprN2piIvN6GkIApCMMAkCWfJ24e/b06u67u0sSrqa3\ng0dHazU2NqHjjkto+3b/t/E1a5o0OkpDCIC5EQYBIEtjY/5vmaOjNUX/tfv7B7Vp0zY9//zMMBrd\nENLSQkMIgOwQBgEgS/X1nn1XpbZeCym9ISSRSOjOO5N6+eXM7eDW1o7II2HYDgaQDcIgAGSpszOh\nPXuODl7HHHOTrrmmcFuvuTaEnHQSDSEA5ocwCABZmg5Yd9/drYMHa/SLX0zqjTcu0cREfsErlxnB\nDQ0bNTqaeQ8aQgDMF+cMAkCevvIVadMmaflyadcu6Zhjsv9aXwVw4cLpGcFfzrj+7LOv0ujokozt\nYM4HBJCN2c4ZJAwCQJ4mJqTzzpN+9Svpttukm2/O/ms3bLhZyWTmodCpGcHf81zfPdUQct+M7eD1\nBEEAWSEMAkCRbN8urVuXqgr+5jfSsmVHfz59K/jyyxM6cKBNXV09euONnoz7tbZepbExKoAACosJ\nJABQJB/4gPRXfyV9//uDeve7k1q7NhX6OjsTkhTRDCJJ/s5kZgQDKDUqgwAwT/fcM6grrkhNJpn2\n1rd2ybnXtXfvVzOuX7KkW5/5zHp973vbqAACKAkqgwBQRN/9blIzg6AkvfRSr6TLvdeffXaNenvb\ndMEFogIIIHaEQQCYp6jJJHV1Yxofz/z49CHVHAkDoBwsiHsBAFDpoiaTnHPOYq1c2XXUx1LTQdaX\nYlkAkBWeGQSAefKdGTj9/J8kjoMBEDuOlgGAIuvvHyT0AShbhEEAAICAzRYGeWYQAAAgYIRBAACA\ngBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgYYRAAACBghEEAAICAEQYBAAAC\nRhgEAAAIGGEQAAAgYIRBAACAgBEGAQAAAkYYBAAACBhhEAAAIGCEQQAAgIARBgEAAAJGGAQAAAgY\nYRAAACBghEEAAICAEQYBAAACRhgEAAAIGGEQAAAgYIRBAACAgBEGAQAAAkYYBAAACBhhEAAAIGCE\nQQAAgIARBgEAAAKWdxg0s0vM7Ndm9lszu6GQi0LhDAwMxL2E4PEalAdeh/LA61AeeB3iV06vQV5h\n0MxqJH1F0iWSWiV9zMzOKuTCUBjl9IctVLwG5YHXoTzwOpQHXof4ldNrkG9l8D2SnnXOveCcG5f0\nXUn/o3DLAgAAQCnkGwZPk7R3xvv7pj4GAACACmLOudy/yOx/SbrEOffpqfcvk/Re59ymGdfkfmMA\nAAAUhXPOfB+vzfN+L0taNuP9ZUpVB+f8BQEAAFA+8t0mfkzSmWa23MwWStoo6T8LtywAAACUQl6V\nQefchJldI2mbpBpJX3PO7S7oygAAAFB0eT0zCAAAgOpQlAkkHEgdPzP7upkdMLNfxr2WUJnZMjPb\nbmZPm9mvzKwz7jWFyMwazOznZvakme0ys3+Me02hMrMaM3vCzP4r7rWEysxeMLOnpl6HX8S9nlCZ\nWZOZ/cDMdk99X3pfrOspdGVw6kDq30i6WKlGk0clfYxt5NIyswslHZL0DefcuXGvJ0RmdrKkk51z\nT5rZYkmPS/qf/L9Qema2yDl32MxqJT0g6R+ccw/Eva7QmNnfS3qnpOOccx+Nez0hMrPnJb3TOfeH\nuNcSMjP7N0k7nHNfn/q+dKxzbiiu9RSjMsiB1GXAOfczSa/HvY6QOef2O+eenHr7kKTdkk6Nd1Vh\ncs4dnnpzoVLPOfMXYYmZWYukD0n6V0mcNhEv/vvHyMwaJV3onPu6lOrDiDMISsUJgxxIDaQxs+WS\nzpP083hXEiYzW2BmT0o6IGm7c25X3GsK0J2SrpP0p7gXEjgn6Sdm9piZfTruxQRqhaTXzOweM9tp\nZv/HzBbFuaBihEE6UoAZpraIfyBp81SFECXmnPuTc+7tkloktZlZe8xLCoqZfVjSq865J0RVKm7v\nd86dJ+mDkq6eeqQIpVUr6R2Svuqce4ekEUmfi3NBxQiDcx5IDYTCzOok/bukbznnfhj3ekI3tRXT\nL+ldca8lMBdI+ujU82rfkbTOzL4R85qC5Jx7Zern1yT9h1KPdqG09kna55x7dOr9HygVDmNTjDDI\ngdSAJDMzSV+TtMs59+W41xMqM1tiZk1Tbx8jab2kJ+JdVVicczc555Y551ZI+mtJP3XOfTLudYXG\nzBaZ2XFTbx8rKSGJEydKzDm3X9JeM1s99aGLJT0d45LyHkcXiQOpy4OZfUfSRZLeYmZ7JX3eOXdP\nzMsKzfslXSbpKTObDh83OufujXFNITpF0r+Z2QKl/gH8Tefc/TGvKXQ8ThSPZkn/kfp3qmolfds5\nl4x3ScHaJOnbU0WzPZI+FediOHQaAAAgYEU5dBoAAACVgTAIAAAQMMIgAABAwAiDAAAAASMMAgAA\nBIwwCAAAEDDCIAAAQMD+P0+PNg5bVEI3AAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pyplot.figure(figsize=(11,7), dpi=100)\n", "pyplot.plot(x,u, marker='o', lw=2)\n", "pyplot.xlim([0,2*numpy.pi])\n", "pyplot.ylim([0,10]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is definitely not the hat function we've been dealing with until now. We call it a \"saw-tooth function\". Let's proceed forward and see what happens. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Periodic Boundary Conditions\n", "\n", "One of the big differences between Step 4 and the previous lessons is the use of *periodic* boundary conditions. If you experiment with Steps 1 and 2 and make the simulation run longer (by increasing `nt`) you will notice that the wave will keep moving to the right until it no longer even shows up in the plot. \n", "\n", "With periodic boundary conditions, when a point gets to the right-hand side of the frame, it *wraps around* back to the front of the frame. \n", "\n", "Recall the discretization that we worked out at the beginning of this notebook:\n", "\n", "$$u_i^{n+1} = u_i^n - u_i^n \\frac{\\Delta t}{\\Delta x} (u_i^n - u_{i-1}^n) + \\nu \\frac{\\Delta t}{\\Delta x^2}(u_{i+1}^n - 2u_i^n + u_{i-1}^n)$$\n", "\n", "What does $u_{i+1}^n$ *mean* when $i$ is already at the end of the frame?\n", "\n", "Think about this for a minute before proceeding. \n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "collapsed": false }, "outputs": [], "source": [ "for n in range(nt):\n", " un = u.copy()\n", " for i in range(1, nx-1):\n", " u[i] = un[i] - un[i] * dt/dx *(un[i] - un[i-1]) + nu*dt/dx**2*\\\n", " (un[i+1]-2*un[i]+un[i-1])\n", " u[0] = un[0] - un[0] * dt/dx * (un[0] - un[-2]) + nu*dt/dx**2*\\\n", " (un[1]-2*un[0]+un[-2])\n", " u[-1] = un[-1] - un[-1] * dt/dx * (un[-1] - un[-2]) + nu*dt/dx**2*\\\n", " (un[0]-2*un[-1]+un[-2])\n", " \n", "u_analytical = numpy.asarray([ufunc(nt*dt, xi, nu) for xi in x])" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "collapsed": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoMAAAGnCAYAAADMu10zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlYlWX+x/HPAwiIIC4IgiIq7mhlmdlYRpbgpDYtY1mW\ntjtlLjn1mya1GNM229RpqpmmmWrabWoqZwQ1EbPSFnNNjeOGgpZLAiog+Pz+OHAEz0EFzsI5z/t1\nXVzBWe77Rhz5zPd5vvdtmKYpAAAAWFOQrxcAAAAA3yEMAgAAWBhhEAAAwMIIgwAAABZGGAQAALAw\nwiAAAICFnTIMGobxqmEYew3DWFftsVaGYSwyDGOLYRhZhmG08PwyAQAA4Amnqwz+Q9LQkx57UNIi\n0zS7SVpS+TUAAAD8kHG6TacNw+go6RPTNPtUfr1J0iWmae41DKOtpGzTNHt4eqEAAABwv/rcMxhn\nmubeys/3Sopz43oAAADgRSENebNpmqZhGC5Li7U9DgAAAO8zTdNw9Xh9KoNVl4dlGEa8pJ9OMSkf\nPv545JFHfL4Gq3/wM2gcH/wcGscHP4fG8cHPwfcf3v4ZnEp9wuDHksZWfj5W0kf1GAMAAACNwOm2\nlnlb0heSuhuGkWcYxq2SnpA0xDCMLZIGV34NAAAAP3TKewZN07yhlqcu98Ba4AGpqam+XoLl8TNo\nHPg5NA78HBoHfg6+15h+BqfdWqbeAxuG6amxAQAAcOYMw5BZSwNJg7qJAQBA42AYLn/Pw4LqWowj\nDAIAECC4Iof6/J+C+nQTAwAAIEAQBgEAACyMMAgAAGBhhEEAAAALIwwCAAC4wfLly9WjRw+Pz5OR\nkaGbb77ZbeMRBgEAgMe99dZb6tevn6KiopSQkKArrrhCK1as8PWyHLKzs5WYmFin9wQFBWnr1q2O\nry+++GJt2rTJ3Utz4u5thNhaBgCAALZgQY7mzs1SaWmIwsLKNXFimoYNG+TVMZ599lk9+eSTevnl\nl5Wenq7Q0FAtXLhQH3/8sQYOHFjXb6lR8cV2Pm6f0zRNj3zYhwYAAN7g6vfup58uM5OTHzIl0/GR\nnPyQ+emny8543IaO8csvv5iRkZHm/PnzXT5fUlJiTpo0yUxISDATEhLMyZMnm6WlpaZpmubSpUvN\ndu3amU899ZTZpk0bMz4+3vzwww/NBQsWmF27djVbtWplPv74446xHnnkEfPaa681r7/+ejMqKso8\n99xzzTVr1jieNwzDtNlsjq/Hjh1rTps2zTx8+LAZHh5uBgUFmZGRkWZUVJRZUFBgrly50hwwYIDZ\nokULMz4+3rz33nvNsrIy0zRN8+KLLzYNwzCbNWtmRkZGmu+99565dOlSs3379o7xN27caF5yySVm\nixYtzJSUFPPjjz+uMfc999xjDhs2zIyKijIvuOCCGmubOHGimZiYaDZv3tw877zzzOXLl9f4Pm+6\n6SaXf5615a/Kx11mNi4TAwAQoObOzZLNNqvGYzbbLA0fvkiGoTP6GD7c9Rjz5i06ozV8+eWXKikp\n0dVXX+3y+VmzZmnVqlVas2aN1qxZo1WrVmnmzJmO5/fu3avS0lIVFBRoxowZuuOOO/Tmm29q9erV\nWr58uWbMmKEdO3Y4Xv/xxx/ruuuu08GDB3XjjTfqqquuUkVFhcu5DcOQYRiKiIjQwoULlZCQoKKi\nIhUWFqpt27YKCQnRnDlztH//fn355ZdasmSJ/vKXv0iScnJyJElr165VUVGRRo4cWWPsY8eOacSI\nERo6dKh+/vlnzZs3T6NHj9aWLVscr3n33XeVkZGhgwcPqkuXLpo6darjuf79+2vNmjWO72PkyJEq\nKys7oz/zuiIMAgAQoEpLa7sbLLgOo7geo6TkzMbYv3+/YmJiFBTkOnK89dZbevjhhxUTE6OYmBg9\n8sgjeuONNxzPN2nSRFOnTlVwcLCuv/56HThwQJMnT1azZs3Uq1cv9erVS2vWrHG8vl+/frrmmmsU\nHBysKVOmqKSkRF999VWt6zMrL7lW/be6c889V/3791dQUJCSkpJ01113admyZWf0fX/11Vc6fPiw\nHnzwQYWEhOjSSy/V8OHD9fbbbztec80116hfv34KDg7W6NGj9f333zueGz16tFq2bKmgoCBNmTJF\npaWl2rx58xnNXVeEQQAAAlRYWLnLx9PTK6pd9D31R1qa6zHCw11X207WunVr7du3T8ePH3f5fH5+\nvpKSkhxfd+jQQfn5+TXeX9Uw0bRpU0lSXFyc4/mmTZuquLjY8XX79u0dnxuGofbt29cYry62bNmi\n4cOHKz4+XtHR0Zo6dar2799/Ru/Nz893akhJSkpyrMUwjFN+H08//bR69eqlFi1aqGXLljp06JD2\n7dtXr+/jdAiDAAAEqIkT05ScPLXGY8nJD2nChCFeG+PCCy9UWFiYPvzwQ5fPJyQkaPv27Y6vd+7c\nqYSEhDNe38ny8vIcnx8/fly7du1yjBcREaEjR444ni8oKHAETVcdunfffbd69eql3NxcHTp0SLNm\nzao11J4sISFBeXl5NSqOO3bsULt27U773uXLl2v27Nl6//339csvv+jgwYOKjo72WLMK3cQAAASo\nqo7fefOmq6QkWOHhFZowYWidOoEbOkZ0dLRmzJih8ePHKyQkREOGDFGTJk20ePFiZWdn64YbbtDM\nmTN1/vnnS5JmzJjRoD30vv32W3344YcaMWKE5s6dq/DwcA0YMECSdM455+jNN9/UzJkztWjRIuXk\n5Kh///6S7NXG/fv3q7CwUM2bN5ckFRcXKyoqShEREdq0aZNefPFFxcbGOuaKi4uTzWZT586dndZx\nwQUXKCIiQk899ZSmTJmiFStW6NNPP1VGRoakU3cEFxUVKSQkRDExMSorK9MTTzyhwsLCev+ZnA5h\nEACAADZs2KA6byXj7jGmTJmitm3baubMmRo9erSioqLUr18/TZ06VX379lVhYaHOOussSdJ1112n\nadOmOd57csXuVHvsGYah3/zmN3r33Xc1duxYde3aVf/+978VHGy/v3HOnDkaO3asXnjhBV111VU1\nmlp69OihG264QZ07d9bx48e1ceNGPf3007rrrrv01FNPqW/fvho1apSWLl3qeE9GRobGjh2ro0eP\n6m9/+5vatGnjWF9oaKg++eQT3XPPPXr88cfVvn17vfHGG+rWrZtjrbV9b0OHDtXQoUPVrVs3NWvW\nTPfdd586dOhQ43Xu3GvQ8FTJ0TAM01NjAwCAmgzD8Mmed43Jn/70J+Xm5tZoQLGa2v4eVD7uMkFy\nzyAAAAgIVg/D9UUYBAAAAcHdl0+tgsvEAAAEAC4TQ+IyMQAAAOqIMAgAAGBhhEEAAAALIwwCAABY\nGGEQAADAwgiDAADAL2VkZNT76Lrly5erR48eDV5Dx44dtWTJkgaP40uEQQAA4BWpqalq1aqVysrK\n3DJeXfYUDAoK0tatWx1fX3zxxdq0aZNb1uDvexsSBgEAgMdt375dq1atUmxsrD7++GO3jFnXfRXZ\nh9E1wiAAAPC4119/XZdffrluvvlmvfbaa47Hb7nlFo0fP17Dhw9X8+bNNWDAgBoVvEmTJqlDhw6K\njo5Wv3799Pnnn9cYt6oqN2zYMP35z3+u8dxZZ52ljz76SJdccokk6eyzz1ZUVJTef/99ZWdnKzEx\n0fHavLw8XXPNNYqNjVVMTIwmTJggSbLZbBo8eLBiYmLUpk0b3XTTTTp06JB7/3B8jDAIAAA87vXX\nX9f111+v6667TpmZmfr5558dz7377rvKyMjQwYMH1aVLF02dOtXxXP/+/bVmzRodPHhQN954o0aO\nHFnjMnNVte+WW27Rv/71L8fja9asUX5+voYPH65ly5ZJktauXauioiKNHDmyxtoqKio0fPhwderU\nSTt27NDu3bs1atQox/NTp05VQUGBfvjhB+Xl5SkjI8Otfza+FuLrBQAAAM8z/uSe+9rMR+p+qfXz\nzz/X7t27deWVVyoqKkq9evXSm2++qcmTJ0uSrrnmGvXr10+SNHr0aE2ZMsXx3tGjRzs+nzJlimbO\nnKnNmzerT58+NeYYMWKExo0bJ5vNpuTkZL3xxhsaNWqUQkJOH3VWrVqlgoICzZ49W0FB9jrZwIED\nJUnJyclKTk6WJMXExOi+++7TjBkz6vxn0JgRBgEAsID6hDh3ee2115SWlqaoqChJ0siRI/Xaa685\nwmBcXJzjtU2bNlVxcbHj66efflqvvvqq8vPzZRiGCgsLtW/fPqc5wsPDdd111+mNN97QI488onfe\neUcffPDBGa0vLy9PSUlJjiBY3d69ezVp0iR9/vnnKioq0vHjx9WqVas6ff+NHWEQAAB4zNGjR/Xe\ne+/p+PHjio+PlySVlpbq0KFDWrt27Sk7cZcvX67Zs2frs88+U0pKiiSpVatWtTaCjB07VmPGjNHA\ngQMVERGhCy644IzWmJiYqJ07d6qiokLBwcE1nnvooYcUHBys9evXq0WLFvroo48c9xMGCu4ZBAAA\nHvPRRx8pJCREP/zwg9asWaM1a9bohx9+0EUXXaTXX3/9lO8tKipSSEiIYmJiVFZWphkzZqiwsLDW\n11944YUyDEP333+/xowZU+O5uLg42Ww2l+/r37+/4uPj9eCDD+rIkSMqKSnRF198IUkqLi5Ws2bN\n1Lx5c+3evVuzZ8+u459A40cYBAAAHvP666/rtttuU/v27RUbG6vY2FjFxcXp3nvv1ZtvvqmKigqn\n6mDV10OHDtXQoUPVrVs3dezYUU2bNlWHDh1qvO7k944ZM0br1q3TTTfdVOPxjIwMjR07Vi1bttT8\n+fNrvDc4OFiffPKJcnNz1aFDByUmJuq9996TJD3yyCP67rvvFB0drREjRujaa6/1+30FT2Z4as8d\nwzBM9vMBAMA7DMNgHz1Jb7zxhv72t78pJyfH10vxidr+HlQ+7jLFUhkEAAAB4ciRI3rhhRd01113\n+XopfoUwCAAA/F5mZqZiY2MVHx+vG2+80dfL8StcJgYAIABwmRgSl4kBAABQR4RBAAAACyMMAgAA\nWBgnkAAAECACbf87eAdhEACAAEDzCOqLy8QAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGE\nQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsj\nDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgY\nYRAAAMDC6h0GDcP4o2EYGwzDWGcYxluGYYS5c2EAAADwvHqFQcMwOkq6U9K5pmn2kRQsaZT7lgUA\nAABvCKnn+wolHZMUYRhGhaQISbvdtioAAAB4Rb0qg6ZpHpD0jKSdkvIl/WKa5mJ3LgwAAACeV6/K\noGEYyZImS+oo6ZCk9w3DGG2a5pvVX5eRkeH4PDU1VampqfVdJwAAAM5Qdna2srOzz+i1hmmadZ7A\nMIzrJQ0xTfOOyq9vljTANM3x1V5j1mdsAAAAuJdhGDJN03D1XH27iTdJGmAYRlPDMAxJl0vaWN8F\nAgAAwDfqe8/gGkmvS/pG0trKh//qrkUBAADAO+p1mfiMBuYyMQAAQKPgicvEAAAACACEQQAAAAsj\nDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgY\nYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDC\nCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAW\nRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACw\nMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACA\nhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAA\nLIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAA\nYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWVu8waBhGC8Mw5huG8YNhGBsNwxjgzoUBAADA80Ia\n8N45kv5rmuZvDcMIkdTMTWsCAACAlximadb9TYYRLWm1aZqdT/Easz5jAwAQ6MqPl2vDTxu0avcq\n5R7IVbfW3XRO23OUEpui8JBwXy8PAcgwDJmmabh8rp5h8BxJL0vaKOlsSd9KmmSa5pFqryEMAgAg\naU/xHuXsyNHKXSu1Kn+VVhesVmJ0ovq366+urbpqy/4tWrN3jX7c/6OSWyXr7Liz9avEX+l3/X6n\nIIPb+9FwpwqD9b1MHCLpXEn3mqb5tWEYz0t6UNLD1V+UkZHh+Dw1NVWpqan1nA4AAP9xuOywcnbk\naNHWRVq8dbHyCvM0KGmQBrQboIxLMtQvoZ+iw6Od3ldaXqqNP2/U93u+12OfP6burbvrss6X+eA7\ngL/Lzs5Wdnb2Gb22vpXBtpK+NE2zU+XXF0l60DTN4dVeQ2UQAGAJx83jWl2wWpm2TGXZsvRN/jc6\nL+E8Dek8REM6D9F5CecpJKhu9Zd5K+fpq91f6c1r3vTQqmElbr9MXDlojqQ7TNPcYhhGhqSmpmn+\nodrzhEEAQMAqKCpQli1LmbZMLdq6SDERMUpPTldacpoGJQ1SZGhkg8bff2S/kucma9ukbWrZtKWb\nVg2r8lQYPFvSK5JCJdkk3Wqa5qFqzxMGAQABo7S8VJ/v/FyZtkxl2jKVdyhPl3W+zBEAO0R3cPuc\n18+/XqlJqbr7/LvdPjasxSNh8AwmJQwCAPyWaZrasn+LI/wt37FcKbEpSk9OV3pyus5vd36dL/3W\nVWZupqYtnaav7/zao/Mg8BEGAQA4A4dKDmnJtiXKzLUHwAqzwhH+Lut8mVo1beXV9VQcr1DHOR21\n4MYFOivuLK/OjcBCGAQAwIWK4xX6tuBbR/hbs3eNBiYOtAfALunqGdNThuHy96fXTP9suorLivXc\n0Od8ug74N8IgAACV8ovyHeFv8dbFiouMc1T/BiUNUtMmTX29xBpsB2y68O8XateUXQoNDvX1cuCn\nCIMAAMsqKS+xN35UBsBdhbt0eefLHY0fidGJvl7iaaX+M1UT+k/Qtb2u9fVS4Kc8sek0AACNzoIF\nOZozN1O/BB1SYdwWRZ5dqC0lG9U7trc6He+hpovPU5/9V+lQ2HG1ndhViec2/iAoSbf1vU2vfv8q\nYRAeQRgEAPidBQtyNHdulkpLQxQWVq7b7rlQa4q+058X/ltFvQ5KRoVkS1f0+wc05lfTZe6K0rvv\nZmr//lmOMWy2qY7Pq481cWKahg0b5Itvq1bX9rxWkxdOVn5RvhKiEny9HAQYLhMDABqtk0PfxIlp\nkqSJk/+nrUevlpIzpS6ZUtuvFbw7ThWb75Ns6dLPPSVVXRGbLsmUNNNp/JiYOxQcHKe9e0+ExOTk\nqZozJ73RBcK7PrlLnVt21oMXPejrpcAPcc8gAMDvLFiQo0mTMmWzVQa1qHw1P/culbbfrNJ2B6Xi\nOHvwyx0q7bhYKv+dpNecxmnTJkPBwdKePRkuZhkl6R2nR9PTp2vhwkfd+e002Fe7vtKYD8do872b\nfd7hDP/DPYMAgEbNVQXwyWf+K5t5mZR2v70C2Hy3CrdeLm1sI33ymVRY836/Zs1Kdfiw89jnnlsh\n0zS1Z4/zc+HhYSopcX68pCTYTd+Z+1zQ7gKFBIVoRd4KXdThIl8vBwGEMAgA8JraLvvaK4AzpZjN\nUnKmFr1zm8wBeVLnHCk3Xfrkb9Lu8yUzWGFho1Ra6tz40a1bpAoLp56oJEpKTn5IEyYMlWS/R/Dk\n55o3b6bVq53XWVZW4ebvvOEMw7A3kqx+lTAIt+IyMQDAK5wu+0pq2/H3Ohq/XodikuzVP8O0hz9b\nurT1TankQ6dx+va9Q4WFcU7Bbs4ce+ibN2+RSkqCFR5eoQkThjju/VuwIMfpOUlOa5IeUmjoUD3w\ngPT1142rsWRv8V71eKGHdk7eqaiwKJ+uBf6FewYBAF7lqgI4d26Wshb9SUr4xt70kZwpxa2VdjaX\nbA/YQ+C+Hqpq/OjVa5xKS2PqHPrqs9aqsUJDK1RRMUSffSZJmZIaX2PJr/7+Kz1x+RMalNS4GlzQ\nuBEGAQBe46rxIzzlLpUm/iiz036puO2J6t+OixUWfKtKS103cUyYMMRtoe9MmabUvfs0/fijc/dx\nY2gsGfn+SP225291fe/rfboO+BcaSAAAHnFyBfDuu9P06OP/lc28/ETjR1S+SrZeLuU2lzKXSIXt\na4zRq0/t9/oNGzbI65U4w5ASEkL044/OzzWGxpL4yHgVFBf4ehkIIIRBAEC9nKgAnmj8yJp/m3RZ\nnvRTZePHx69I+f0kM1jJyeOkNi/KVlgz9D366BhJ0rx506tVAIf69HJsWFi5y8fDw33fWBIfGa+C\nIsIg3IcwCAA4pZOrf7ffnqZS4yxNnjtXB3q1kkZ0PNH4sfoJGR++KfOIc+NHly6xlZd9XYc+X9+L\nV93EiWlO3ceJiSc6k30pISpBG/dt9PUyEEAIgwAASafZ9mXrDEfjR1bmaCluv9QxQcodL301ucaJ\nH71SFqmkpPFc9q2PqjXOmzdda9cGq6CgQh07+rZaWSU+isog3IsGEgCAy21fWnWcqKMJ63Q0IVbq\nvFgqjnc0frQ49JliWoQoN9d1k4UvGj88JS9P6tJFOnZMWr9e6tXLt+tZt3edRn0wShvu2eDbhcCv\n0E0MAHBwVQF8/vksLc6eJnVYfmLbl6h8aWuUZHtYsqXVaPy45JIMPfDAYKcAWbX1i78Gv9rcc4/0\n4ovSqFHS22/7di37juxTt3nddOAPB3y7EPgVwiAAQNLJFUBTitmkJj0n6liHH6UO+6Wf+tirf7lD\npfx+aho+WkeP1n52r6uNnAMtCEr26mByslRe7vvqoGmaCp8VrkMPHlJ4SLjvFgK/QhgEAAs6uQJ4\n661pmvn0x9pwdIC98pecVe3Ejw3S1gVSScsaY5zqtI9ADH2ncvfd0ksvSTfcIL31lm/X0uG5Dsq5\nNUcdW3T07ULgNwiDABDAzrTxQ13+KsX+JO28zL7hc7UTP3r0GKdjxzx/2oc/27nTfu9gebm0YYPU\ns6fv1nLBKxfo+fTndWHihb5bBPwKm04DQIBy1fjxzZbKxo++sdK1sVJRgj38Lf2Hgne/rIrS+U7j\nJCX5z7YvvtKhg3TbbdLLL+do8OAsde/uuzOLE6ISlF+U79U5EbgIgwDgJ1xVAJ97Lku2HdOlzosc\njR8HovKlrc2l3DFS5rNSUTvHGD1S5vv9ti++dMEFOXr55Uzt2TNLe/bYH7PZpkrybmDmFBK4E2EQ\nAPyAq8aPpU9WNn48MFfa28de/as88cPe+HGr0zjt25+6AohTe+edLEmzajxms83SvHnTvR8G2WsQ\nbkIYBIBG5uQK4C23pGnm0/+RLexCacSd9gqgpGO5Q6XvjkgffC+VtKgxRo8ejeu830BRWur616a3\nzyyOj4rXF3lfeHVOBC7CIAD4yGkbP9p9bT/vd8mNUvpP0s6N9urfl1NqNn60m+0U+hrjeb+BoLGc\nWcxlYrgTYRAAfIDGD//k6sziqoqrN3EkHdyJrWUAwINcVf+GDRukyy+fpiXLqk78WGi/9Bu5x37i\nR27liR/VGj9SUsappMT11i8EPu9asCBHM2Ys0qpVwYqMrNA773h/q509xXt01otn6acHfvLqvPBf\n7DMIAD7gXP0zFd3ld4roE6qCiIVSh73S3rMc5/0q/zxFNB2tI0dcn/gRSOf9+rvSUqllS+noUWnv\nXik21rvzVxyvUPiscB156IiaBDfx7uTwS+wzCAAe5qoC+PTTWbLtvl/qNd9+4keXTB0yDR2ytZa+\nO1v64BWnxo/u3Wn88AdhYdJFF0mLFknZ2dJ113l3/uCgYLWJaKO9h/eqffP2p38DcAqEQQBooBoV\nQKNCave1ljz3f6pIzpUunCPtvKiy8eP30r7u6tnzT5o2bbAefpjGD382eLA9DH72mffDoHTivkHC\nIBqKMAgAdXByBXDkyDQ9+eJ82ZqfLY28Tuq0RCpKUIUtXVpqSDszpfLwGmN06FChG28cpOjo2kMf\n4a/xGzzY/t/PPvPN/HQUw124ZxAAXKht25cJEzK1LW+6lJRTeen3dSmyWNr6m8p7/040fvTqNU6l\npTR9BKrycql1a6mw0H5ucWKid+e/8+M71S+hn8b1G+fdieGXuGcQAOrAVePHyq3jdLTdBpUNjJIS\nY080fvxngUJ+nq3yMuemj8RETvsIZCEh0iWXSJ98Yq8Ojh3r3fnjo+I5nxhuQRgEYGmnbvx433He\n7yEzSMoNk779vTT/nRqNH91TWnHer0UNHuy7MJgQlaDvCr7z7qQISIRBAJblsvHj2arGj7n2xo/c\ndOmL+6V93dWs2Q06fPgap3E479e6qt83aJqS4fIinGdwzyDchTAIwBLOrPGjnSpy06WlQdLOhU6N\nH926se0LaurdW4qJkXbtknJzpa5dvTc3p5DAXQiDAALKqRs/pklJy+3n/a6+Vrq8svHjx19LC5+r\n2fiR9CjbvuC0goKkSy+V3n/fXh30ahikMgg3IQwCCBi1Nn4kVG/8ONvR+NFk32wdK6174wfhD9UN\nHnwiDI7zYmNvXGScfj78syqOVyg4KNh7EyPgsLUMAL9T23m/qanTtGzl76XOix2NHzKDJFuYlPuE\ntG1wjcYPzvuFO2zZInXvLrVpI+3ZY68Wekub2W20/u71iouM896k8EtsLQMgYDhV/4wKfZF3m0Je\n+Yd+6bJY+hWNH/Curl2ldu2k3bulDRukPn28N3d8pH17GcIgGoIwCKDRclUBfOKJLNl+vlvq+3d7\n9a/TEhUdD5VNAAAgAElEQVQXtZNym0lrB0k7/07jB7zKMOyXit94w36p2JthMCEqQQXFBeqrvt6b\nFAGHMAjA52pr+nBUAEOOSkk5WvTnB2SelSv96iXJNkT68Qpp4fNSUYL69MnQn54erAceoPED3lc9\nDE6a5L156SiGOxAGAfiUc9OHtGHjQyptvlX7YvtLF6ZLiV9Ie8+WmTtc+s8yKT9TMmveMJ+QUKGr\nrx6k0FDO+4X3XXqp/b/Z2fZj6kK89NuVjmK4A2EQgNe4qgDOnZtlD4JND9gbP5IztTs5SzIPSrnR\n0rfjpPnvOho/UlJ2qST8YZeXfCVx2Rc+kZQkJSdLNpu0erV0/vnemTc+Ml6b9m3yzmQIWIRBAF7h\nVAEMKtcXebfraMJm6fYlUuwGacfFki1dWvF/alL0sI6Vvew0Dk0faKy6dMmRzZal668PUdeuJ7rc\nPSk+Kl5Lty/16BwIfIRBAG53cgVwwgQXjR+dF6v4UKJkK5U+e8HeAVwR5hijd98omj7gNxYsyNH3\n32dKmqVt26Rt2ySbbaokz96eUNVNDDQE+wwCcKsaFcDKxg+j68MyO+dKkYa98cOWLtnSpKIEJSWN\nU3BwjLZudd7rT5LmzVtUrQI4hBCIRik9fZqysma6eHy6Fi581GPzbju4TamvpWrH5B0emwOBgX0G\nAXhE9QpgaGi50tOHaM5bbykvtvuJxo8958i0jZA+WiYVLHRq/OjRg9M+4P9KS13/Oi0p8ezJIPFR\n8dpTvEemacowXP6eB06LMAjgtGrb+mX8+Ezt+On3jsaPRbtHSIMk5d7govEjTyVNXTd+cNkX/i4s\nrNzl4+HhFR6dNzwkXBFNInTg6AG1jmjt0bkQuAiDAE7JVePHip23q6T9BlUMCZHazKvR+BFaPF1l\npTR+wFomTkyTzeb6HldPq9pehjCI+iIMApDkuvp3xRWD9PjjWbL9/Dvp3FfsZ/12XqLDhxIlW4W0\nZI5T40fXlJYqKaHxA9ZS9fd6+vTpWr06WK1bV3jtjOuqjad7x/b2+FwITDSQAHCu/oUcVdOet0rJ\nx3W03VIpQtLWIfbzfm1pUnG8IiNHqbj4Haex0tOnV1YAafyA9WzbJnXuLLVtKxV4aS/om/59k9KS\n0zTm7DHemRB+iQYSAA4nVwDHj0/TrMcyZSu8UbrwWXv1L/ELHd1zjmQLVvB3I1SR94pkBtUYp2tX\nzvsFTpaUJDVrJu3ZI+3bJ8XEeH5OtpdBQxEGgQB12vN+K0/8yHp1ojTQJg14y37f3ze/k95/TyqN\n1rnnZijjz4N1333TOe8XOANBQVLv3tLKldK6dSeOqfOkhKgEbf9lu+cnQsAiDAIByNV5v2vX/1Gl\nMZt1MLGPdOkAqc1GaccgKfd2GSv+J3PfAkk1ryC0aVOhESMGKSiI836BM9Wnj3fDYHxUvL7c9aXn\nJ0LAIgwCfu6U5/02z7Of9pGcqT2dl0iHJOV2lZY8Ju0c6Gj86JWyXiXR0zjvF3CDPn3s/123zjvz\nVXUTA/VFGAT8mKvGj+UFt6qk/WZp/L+liJ/tjR9bhkv/m6uwY/eptPRJp3HY9gVwn6owuH69d+ar\n6iYG6otuYsBP1Nb4sWrbDfamjy6ZUvsvpT19pdyDku0fUsG5NRo/+va9Q4WFcU4VQG9tgQFYwc8/\nS7GxUmSkdOiQ/T5CTyoqLVLbZ9qq+I/FnEKCWp2qm5gwCPiBGhXAysYPdXlMSrZJFTH2xo/cdGnb\nYKk0Wp07j5NhxLgMfRLn/QKe1rattHevtHWr1KmT5+eLfCxSu6fsVnR4tOcng19iaxnAj1SvADZp\nUq6BFw/WS5++rr2JSWfc+NG1K+f9Ar7Up489DK5b550wmBCVoILiAsIg6oUwCPiIq8YP05TuuSdT\neYXjHJd+Fx95RjovUsq9rU6NHzR9AL7Tp4+0eLE9DF55pefnq7pvsEdMD89PhoBDGAR8wKnxo8kR\nLS+4TaWJG3X8ymNSs5ftJ31sHiH9d57CyyerpITGD8Bf9K48Gc5rTSR0FKMBCIOAB7mq/g0dOkgz\nHs2UregG6cJnHI0fR/f0lXKDpQ9fdWr8SE5pwXm/gB/xyfYydBSjngiDgIc4Vf+aHtCyfber4h9P\nqHzg5ydO/Pj6Hum9+VJpc0VFjVJRUT+nsagAAv4lJUUyDGnzZqmsTAoN9ex88VFUBlF/hEHADU6u\nAP7ud2ma+dhC2cqGS6mP2Kt/bTaqdMcgKfe4wr68SaX5L+jkxo8uXTjvFwgEERFScrKUmytt2iSd\ndZZn54uPjNf3e7737CQIWIRBoA5Oe95v9E4pOVNZ79wlXbZdOvRf+5Yv1Ro/+vfP0MN/HaxJk5wb\nPzjvFwgcffrYw+C6dV4Ig1QG0QCEQeAMuTrv9/sN/6ejcRtU1LWrNLSnFLHPfuLH5ocUlDlfxw99\n7DROy5YVjnDH1i9A4OrdW/rwQ+/cN5gQlcA9g6g3wiDggqsK4Jw5WbLZZkqx6yu3fVmon9p/Je2J\nkHJ/JX34Ro3Gj54pK1QS4/qSr8R5v0Cg8+axdHQToyEaFAYNwwiW9I2kXaZpjnDPkgDfcm782K9l\n++5QaeIWaco/pYpQ+6Xfr8dL732gcOMulZRMdRqHpg/A2rzZUdwivIXKKsp05NgRRTSJ8PyECCgN\nrQxOkrRRUpQb1gJ4Xe2NH8OkSx+2VwBjNtkbP2xh0vKl0oEuqt740bMvTR8AnHXpIoWFSTt32s8o\njvbg4SCGYahtZFsVFBUouVWy5yZCQKp3GDQMo72kKyTNkjTFbSsCPKD+jR+POxo/unYdp+Mt/ynb\nAZo+AJxeSIjUq5e0erX9UvHAgZ6dr+pSMWEQddWQyuBzkh6Q1NxNawE84owbP2xp9saPhfN1vNC5\n8aNzZ877BVA3vXvbw+C6dV4Ig1Hxyi/K9+wkCEj1CoOGYQyX9JNpmqsNw0it7XUZGRmOz1NTU5Wa\nWutLAbdwVQGcO7eq8WOdfb+/5MzTN35w2gcAN/BmE0lCJB3FOCE7O1vZ2dln9FrDNM06T2AYxmOS\nbpZULilc9urgB6Zpjqn2GrM+YwP15arxI6zX7fbGj86HpIow+6VfW7q0bXBl48c7TuOkp0+vrAAu\nqlYBHEIIBFBnCxdKv/61NGiQtGyZZ+d6bPljKiwt1BOXP+HZieCXDMOQaZqGy+caGtgMw7hE0v0n\ndxMTBuFJrho/Zj2+UN/uGeao/ilmk7RjkGTbLeW+69T40bfvHSosjHOqAM6Zw/1+ANxj926pfXup\nZUtp/377EXWe8o/V/1D2jmy9dtVrnpsEfutUYdBd+wyS+uA1NSqA0TukLpWNH4O3S7/8z175W/yE\nlPcrGj8A+FRCgj0IHjwo5edL7dp5bq74qHguE6NeGhwGTdNcJsnDxW9Y0cnVv9tvT5MZ0k/jZ/9F\n+7vGSb/uITU9INnsJ34Y/5svs4jGDwCNh2HYm0iWL7c3kXg0DLLxNOqJE0jgc6fe9qXyxI8uC5W1\nYKzUbq/UI1ay3SV98Ka0p6+j8aMXjR8AGqE+fexhcP16aehQz81DNzHqizAIn3K17ct3P9ynownr\ndbh3e+mqdlJ5uJQ7VFr5vKIPrlBMVKhstoecxuLEDwCNkbdOIomJiFFRaZFKy0sVFhLm2ckQUAiD\n8BpXFcDnn8+SbdufpMQV9saPLgu1r/VmaUdzKfcqafnUysYPu3MuWa0HHhisSZOoAALwD4WFOZKy\n9MEHIdqzx/5vnyf+nQoyghQXGac9xXuU1CLJ7eMjcBEG4RVOFcDoHcqedbfKErdID7wg/dLR3vix\n6Ckp71dqGjpGR4+OdxonPLzC8Y8oFUAAjd2CBTl66aVMSbN0+LCUlSXZbPazzD3xb1bVfYOEQdQF\nYRBud3IF8M470zTzyQWyBV0iDZ1krwA2PaAyW5q0qbX038+l4rY1xujRo/bzfiVRAQTgF+bOzdK2\nbbNqPGazzdK8edM9EwbpKEY9EAZRb2fW+JGprI/HSkN2SwVf2at/1Ro/uncfp/K4ebIVs+0LgMBT\nWur612xJSbBH5qOjGPVBGES9uGr8WL3pPh2Jr2z8+E17e+OHLV1a9byC5r+m40f/7TROx45s+wIg\ncIWFlbt8PDy8wiPzxUfSUYy6IwzitFxVAOfMqWz86PC5/bSPLpn6+RSNHz1T/su2LwAsZ+LENNls\ntd/y4m4JUQn6atdXHhkbgYswiFNyqgC22K7sx+5WWeKPJzV+PCnlDay18YNtXwBYUdW/cY89Nl1f\nfBGspk0rPHrkZXwUl4lRdw0+m7jWgTmb2O+4avx49MlPtfZQamX1b6HU9KBkS5Nyf5S2/sep8YPz\nfgHAWVmZ1KyZVFEhFRXZP/eEb/O/1R2f3KHV41Z7ZgL4rVOdTUwYhKTqFcATJ36oy1+kdrulggvt\n1b/codKecyQzSN26jVNFRYzL0CdJ8+YtqlYBHEIQBGB5ffrYTyFZuVLq398zc+QX5evcl8/Vnvv3\neGYC+K1ThUEuE1vMydW/W25JU/HxFP3h5ed0sHeLysaPplJuurRyjoJ2/tNl40enTjR+AEBdVIXB\ndes8FwZjm8Vq/9H9Kj9erpAgfsXjzPA3JUCdctuXbX+S2n8lJWcqa+koqXWh1L6dlDvRRePHAho/\nAMAN+vSR3n5bWrvWc3OEBIWoddPW2lu8V+2at/PcRAgohMEA5Grbl29yx+tIwnqV9Gst/baNvfEj\nd6i06C21LM5Sm1ZB2rKFxg8A8BRvnVGcEJWgguICwiDOGGHQz7mqAD73XJZsOx+Sui6wn/aRnKkD\nTQ9KtkjphzukBX+p0fhx1iXZnPcLAB521ln2/65bJ5mmZLi8e6vhOIUEdUUY9GM1K4CmFLdOn82+\nT+VJP0r3z5EKzrPf+/fB29KecxTR9EYdOXKz0zic9wsAnpeYKEVHS/v2SXv3Sm3bnv499cEpJKgr\nwqCfcNX4MfPZj2Rrer501S1ScpZ0LELluenSylJp2wapLKrGGN27c94vAPiKYUi9e0srVtirgx4N\ng1QGUQeEwUbm1I0fGVL7lScaPy7fL20fat/2JWeao/GjR49xOpb4hFPo47xfAPCtPn1OhMEhQzwz\nR3xUvNbu9WCXCgIOYbARqbXxo916lfRrJY1sIx3sbL/0u+gtBef/WRVl853GSUpi2xcAaIy80UQS\nHxmvhbkLPTcBAg5h0Efq3vhxp/TpS9LhOMcYPVLeZtsXAPAjVWHQk9vLVHUTA2eKMOgDzo0fa8+g\n8eMmp3HY9gUA/EtVGNy40X40XXCw++egmxh1xXF0HnZyBfDWW9P06DMfaWPJ+ZXVP3vjh3LTJdsa\nadv/nBo/OO8XAAJHhw5SXp60aZPUvbv7xy+rKFPkY5EqmVaiICPI/RPAL3EcnYe5uuQ7bNigExXA\nbRn2Ez+6ZCrrs5MaP5ZNlw4mS6LxAwCsoE8fexhct84zYTA0OFTNw5pr35F9im0W6/4JEHAIgw3k\nqulj8+apem9RgT5c+08V9Quv1vgxVMp6W8EF82j8AACL6tNH+u9/7WHwt7/1zBxVl4oJgzgThME6\ncFUBnDs3yx4EQ4uljtlScqZ2dMnU603mSC2SpR8ecNH48RaNHwBgUd7oKE6ISlB+Ub7Obnu25yZB\nwCAMniHnCqCplTvu1JGETdKYy6R2q6T8fvZLv++/q1bHPlJcbIV++IHGDwDACd7aXoaOYpwpwqAL\nriqAzzyTJVvBfVLvtx3bvhwqi5Rskr56VtqeWqPx4/z0f2vChDTO+wUA1NCjhxQSItls0uHDUrNm\n7p+DU0hQF4TBk9SoAAYdk9p/pSVzHlRFl1zpgnnS9kvtnb/LHpYOdlb79uMUEvKVtpeNcIxRPfBJ\nNH4AAE4IDbUHwvXrpQ0bpP793T9HfFS8tuzf4v6BEZAsHQZPrgCOGpWmx158T7aWvaXrr5Y6LZUO\nJKvCli5lhUq7MqWK0BpjpKScvumD8AcAqK5PH3sYXLfOQ2EwMl7Ldixz/8AISJYIg7Wd9ztxYqa2\n5j1kb/zokqms1VdLlxyVbNdKP1xbo/GjV69xKu34Jy75AgAarE8f6e23PXffIBtPoy4CPgzW1vhx\ntN0GlV0UIbVrW63xY7FC9j+u8mNvOI2TmEjTBwDAPTzdRFLVTQyciYAKg7U3fkyW+rwlJdtP/LA3\nfgRLXz1kvwewLNIxRveUlmz7AgDwKE+HwfjIeO0p3iPTNGUYLg+dABwC5jg6p8aPxC8V3O2PquiU\nK7UqsXf72tLtzR8HkxURMUpHjrzjNE56+vTKCuCiahXAIYRAAIDbmKbUooVUWCjt2SPFxZ3+PXXV\n4okW2jppq1o1beX+weF3Au44ujo1fmS6bvzo3j1ShYVUAAEA3mcY9urgihXS2rXSkCHun6PqvkHC\nIE6n0YbBU53362j86LRUSq5q/CixN35s/K306cvSYfsRPD17jlOZi8YPzvsFAPhSVRhct85DYbBy\n4+mU2BT3D46A0ijDoKvzfjf+8JD6DftRWbY3deRiSQltpfzz7Zd931+ikP2PqfzY605jdejAeb8A\ngMbH4/cN0lGMM+TzMHjK834jfpaSF0nJmdqVnKVdZaVSky7Sl9Np/AAA+DWPdxRH0lGMM+PTMOhU\nAQw6pi9236aSdpuluxZKrX60N37kDpWWPaLYJq8rLq5c67aMcBqL834BAP4kPz9HUpa++y5EaWnl\nmjQpza2/s+Kj4rXz0E63jYfA5bUw6KoC+NRTWbLtv0Pq95J925dOS1V8IFmyHZEyX5byLpSON3GM\n0Te9gvN+AQB+b8GCHE2dmilplkxTWrRI2rp1qiT33b4UHxmvlbtXumUsBDavhMEaFcDQYqljtha/\nOEXHz7ZJ/V+RbGk1Gj86dBin4OBMbTt+4n8QnPcLAAgUjtuhqrHZZmnevOnuC4PcM4gz5NEwmJ4+\nTVdeOUSz33hDO9p2lS4aLCV8LeWfr+O5I6X3l0h7F0pmUI339ezJeb8AgMBVWur6129JSbDb5qjq\nJoa1VV2ZPRWPhsGsZjuVtX2Y1D9Est0kfXlfjcaPXr22qrTZdC75AgAsJSys3OXj4eEVbpujqjLI\nKSTW4Op2PEn2K7P775A0q9b3evYycd6FUnaGQg//UWVlf3Z6mvN+AQBWNHFimmy2mve/JyXZiyHu\nEhUaJUkqKitS87DmbhsXjY9TQ25osVYevE2liT+o5IoSKezv0tO1v9+zYfCbuyVJXVNasO0LAACV\nqt///vXXwTpwoEJXX+3eYohhGI7qIGEwcJxcARw/Pk2PPZ4p2+GR0kVP2BtyE77RofzzpdxQacW/\npL1nSQqqdUyvNJCw7QsAADVVFUPmzpUmTZJ++sn9c1TdN9g9prv7B4fX1agANvtJ6rxIWf8cL120\nVTr/fftBHF/+3r4tX1mkwsNHqaTk7NOO6/EwSAUQAIDapdlv7dKiRdLx41JQ7QWcOqOj2D+dXP0b\nOzZNZtCFmvz037SvY6J0WT+pVa607VLJNl7K/kQ6uMBpnJ49I1VYONWpc/1kHu4mnk4FEACAU+je\nXUpMlPLypDVrpL593Tc2HcWNW21NHxMnZmrr1llSy61ScqayPrlT6pgvnRcl5d4qZT5bYy/mlJTV\nKmnlfDveo4+OkWS/HSEzs/Z1eDQMLlz4qCeHBwDA7xmGvTr4979LWVkeCINUBhslp6YPSd+ue0BH\n4jboaI/O0rCuUliR/dLvhkfUYvk3ahMRoR9/nOk01uluxxs2bJAMw/l9VXx+NjEAAFY3ZIg9DC5a\nJP3hD+4bNyEqQWt/Wuu+AVEvriqAzz6bJZttptT2e3vTR5dM7U/4WsqPlHIvkd6bX9n4Yd8W6OxL\ntuiBBwZ75BQ2wiAAAD522WX2CuHy5dKRI1JEhHvG5Z5B33OqADb7SUv33aljSVuk+1+RSptXNn5M\nkbanKiLkDh054vz/CMLDKzx2ChthEAAAH4uJkc47T/rmGyknRxrqpu0GuWfQu06uAI4enaYnZ/9P\ntvKh0mUP2SuArWw6tj1Vyo2Qlv5X+qVTjTG6pzg3fVRV/yTPnMJGGAQAoBFIS7OHwawsN4ZBKoMe\ncfrGD5vUJVNZ/7tDGr5DOrDIXv3LfM7R+NGz5ziVtX5Ftl9qb/rw1nZ8hmmanhnYMExPjQ0AQKBZ\ntkxKTZVSUqT1690zpmmaajqrqQ7+4aCaNmnqnkEtzlXjR6v4+3U0bqOOJnS2V/9CiyVbmmRLV8jO\nd1V+6D9O49h3XBmiefMWVQt9QzwW+gzDkGmaLs8lJAwCANAIlJVJrVpJhw9Lu3dLCQnuGbfTnE5a\nMmaJOrfs7J4BLaS2xo/Pls6Q4tZIXTIdJ35od6Rku89eAazW+JGSMk4lJTFOl33nzPHu1nunCoNc\nJgYAoBEIDbVXBhcssHcVjx3rnnHjI+OVX5RPGKwjV40fn+27U+VJW6Tf/00qaSHZ0qUv7rc3fjS5\nXUeO/J/TOP5wChthEACARiItzR4Gs7LcGAa5b/C0Tq4A3nhjmp582rnxo3zbpZKtlsaP3rU3fjT2\nU9gIgwAANBKeOJquqjJoda4u+Q4bNkgLFuQ4N34sPKnxY+Hz0q4Bja7xw124ZxAAgEbCNKUOHaRd\nu6TVq6Vzzmn4mK+veV3vb3xfn9zwScMH81Oumj4SE6cqNf1iffT9qyqKbWO//6/J4UbX+OEuNJAA\nAOAnbr9devVV6cknpf9zvgWtzg6XHVb759prwz0blBDlpq6URsxVBXDu3CxlZc2UjOMnNX58Ie1u\nJ9l+12gbP9yFMAgAgJ94911p1Cj7qSSLF7tnzDs/vlNdWnXRHy5y41l3jZCrCmB0wmQdid+kYx3a\nSMlZJxo/ctPV5vAKxbUK0vr1zuf2+msFsDaEQQAA/MS+fVKbNjkyjCwNHBiiiIgT97fV15d5X+qW\n/9yiTeM3yTBc5gG/46oC+NRTWcpZ8bCU+IXjvF+13Cpti5Ryp9tD4C8dHWNUBb6TA6Q/VwBrw9Yy\nAAD4iZUrcxQWlqnS0ln6/HP7YzbbVEmqdzgZ0H6AgowgfZH3hQZ2GOiupfqMUwWwVa4Wzxyv451+\nlB6YJ+3vZg9+/5sj7RqgDu3vVXDwTm2rFgSrd/pK/tn44S5UBgEAaETS06fZ729zeny6Fi58tN7j\nzl4xW5v3b9YrV77SkOV53ckVwCuvTNPzf/lEuRUXn6j+NTlceel3s7T1Y+lImxpjBNol3/qgMggA\ngJ8oLXX9q7mkJLhB49589s3q+UJPPT/0eUWGRjZoLE+o7bzfe+/N1PYdj1Y2fixU1jc3S78pkHav\ntjd9vPtvaW8fSYZ69Rqn0vjn/XKvP18iDAIA0IiEhZW7fDw8vKJB47aNbKtBSYM0f+N83XLOLQ0a\ny91cNX6sXD9ZRxPWq6xvvPTb+BONHyteVJP8v+vY4Q+cxklMbPynfTRGXCYGAKARcRWMYmIe0j//\n2fBQ89Gmj/Tsl88q59achi6z3lxVAJ98MkvLv3DV+BEt5T7k1PgRaNu+eAPdxAAA+JEFC3I0b94i\n5eUFa+PGCkVGDtH27YPUunXDxj1WcUztn2uv5bcuV7fW3dyz2Dpw1fgR1K2y8SNp/4nGj9x0adcA\nRUbcrOLid5zG4R7AuiMMAgDgh0xTGjJEWrJEGjdOeumlho95f9b9Cg0O1WOXPdbwwU7h5Arg8OFp\nmvPiJ7Idv+jEps9NjlRr/PhEOhJTY4y+fe9QYWEcFUA3IAwCAOCnNm6Uzj5bqqiQvv5aOu+8ho23\n4acNSvtXmnZM3qGQoIa3Dpy28aPt6srw95IUv0faPche+bOl12z8KHV92VcSFUA3cHsYNAwjUdLr\nkmIlmZL+aprm3JNeQxgEAMAN7r9feuYZacAAacUKKSioYeMNeGWAHr7kYV3R9YoGjePq/sbmCZN0\nNGG9jnVoKyUvko62coS/2ho/uOzreZ4Ig20ltTVN83vDMCIlfSvpKtM0f6j2GsIgAABuUFgode8u\n7dmTo5SULMXEnKjC1Scw/fXbvyrLlqX5180/o9e7qv5dccUgXXTRNH2xcrrUYcWJxo/oHdK25pJt\nqmRLo/GjkfD4ZWLDMD6SNM80zSXVHiMMAgDgJvffn6NnnsmUVD1ITdWcOel1DlKHSg4p6fkk/Tjh\nR7Vp1uaUr3Wq/hnH1bTjXVLnUB2Nz5IS90o/p1RW/9Kk3RcoMuImGj8aGY+GQcMwOkpaJinFNM3i\nao8TBgEAcBN3n0wy5sMxatakmW4/93b1ju2t8JBwpwrgrbem6bGn/qt1P10ttVtl3/ql85LKPf+a\nK2RnR5XnvmL/uhoaPxofj51AUnmJeL6kSdWDYJWMjAzH56mpqUpNTW3IdAAAWFZtJ5MUFQW7vIx7\nqtC1YEGOtv4zUrkJK/SvnH+rLPIXtQ1rp/0bInS4+BbpSGsp4VtlfTZK+vXP0v5MaXd/aesQafET\n0qEk9euXoUceH6zJk2c7hb5HHx0jydrn/fpadna2srOzz+i19a4MGobRRNKnkv5nmubzLp6nMggA\ngJvUVhkMCrpDLVvGaf9+58vHklx2+p7c9NEs+g8qbb5V5a2HSW2/lyL2SfnnSbv7K/jnZ1VR4rrp\nY+HCRx17InLZt3HzRAOJIek1SftN07yvltcQBgEAcBNXnbtNmz6ko0d/kfQXp9d36nSHSkriVFBw\n4vXR0VNlmgdVWOj8emms7L/aa6LpIzB44jLxQEk3SVprGMbqysf+aJrmwnqOBwAATqEqeFW/9Hrv\nvUP1xz9+pvXrnV+/bVuxpFdqPHbo0CzZQ5+zqKhSFRU5P96+Pef9Bjo2nQYAwI/Vdvk4JGSsysud\nK33Nmo3S4cPOnb40fQQ2jzWQAAAA35o4MU0221SnENe8eTOtXu38+m7dIlVY6Px6mj6si8ogAAB+\nzlUTh+TcKMIRb9bF2cQAAFgQnb6oQhgEAACwsFOFwQYedQ0AAAB/RhgEAACwMMIgAACAhREGAQAA\nLAXC6PUAAAWZSURBVIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACA\nhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAA\nLIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAA\nYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAA\nAAsjDAIAAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIA\nAFgYYRAAAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAA\nAMDCCIMAAAAWRhgEAACwMMIgAACAhREGAQAALIwwCAAAYGGEQQAAAAsjDAIAAFgYYRAAAMDCCIMA\nAAAWVu8waBjGUMMwNhmG8aNhGH9w56LgPtnZ2b5eguXxM2gc+Dk0DvwcGgd+Dr7XmH4G9QqDhmEE\nS/qzpKGSekm6wTCMnu5cGNyjMf1lsyp+Bo0DP4fGgZ9D48DPwfca08+gvpXB/pJyTdPcbprmMUnv\nSPqN+5YFAAAAb6hvGGwnKa/a17sqHwMAAIAfMUzTrPubDONaSUNN07yz8uubJF1gmuaEaq+p+8AA\nAADwCNM0DVePh9RzvN2SEqt9nSh7dfC0EwIAAKDxqO9l4m8kdTUMo6NhGKGSrpf0sfuWBQAAAG+o\nV2XQNM1ywzDu1f+3dwehUpVhGMf/j97CrkmCQhYKuWkXpLYQRY2woAhpWRBBi1ZSQRSUi7Yta9VG\nvaIpbm4YRFBRSeTGrO7F8l5aRIFCmlAE5tKnxZwghAjEmVd9n99mvpnVwxzmzHu+833vgU+ApcB+\n24vXNVlEREREjN01rRmMiIiIiFvDWJ5AkobU9STNSLog6fvqLF1JWifpuKQzkn6Q9FJ1po4kLZN0\nUtK8pAVJb1Vn6krSUklzkj6sztKVpF8knR6Ow9fVebqStFLSrKTF4by0uTTP9Z4ZHBpS/wjsZLTR\n5BTwTG4jT5akbcAl4JDtB6rzdCRpDbDG9rykO4FvgafyW5g8SdO2L0uaAk4Ar9o+UZ2rG0mvAJuA\nFbZ3VefpSNLPwCbbv1dn6UzSQeBL2zPDeWm57T+r8oxjZjANqW8Atr8C/qjO0Znt87bnh/ElYBG4\ntzZVT7YvD8PbGa1zzh/hhElaCzwB7APSbaJWvv9Cku4CttmegdE+jMpCEMZTDKYhdcRVJN0HbABO\n1ibpSdISSfPABeC47YXqTA29DbwGXKkO0pyBzyR9I+mF6jBNrQcuSjog6TtJeyVNVwYaRzGYHSkR\n/zLcIp4FXh5mCGPCbF+x/SCwFtgu6eHiSK1IehL4zfYcmZWqttX2BuBxYPewpCgmawrYCLxreyPw\nF/B6ZaBxFIP/25A6ogtJtwHvA4dtf1Cdp7vhVsxHwEPVWZrZAuwa1qsdBR6RdKg4U0u2fx1eLwLH\nGC3tisk6B5yzfWp4P8uoOCwzjmIwDakjAEkC9gMLtt+pztOVpNWSVg7jO4BHgbnaVL3Y3mN7ne31\nwNPAF7afq87VjaRpSSuG8XLgMSAdJybM9nngrKT7h492AmcKI13z4+j+UxpS3xgkHQV2AKsknQXe\ntH2gOFY3W4FngdOS/ik+3rD9cWGmju4BDkpawugC+D3bnxdn6i7LiWrcDRwbXacyBRyx/WltpLZe\nBI4Mk2Y/Ac9XhknT6YiIiIjGxtJ0OiIiIiJuDikGIyIiIhpLMRgRERHRWIrBiIiIiMZSDEZEREQ0\nlmIwIiIiorEUgxERERGN/Q1nM7PU4/GDkQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "pyplot.figure(figsize=(11,7), dpi=100)\n", "pyplot.plot(x,u, marker='o', lw=2, label='Computational')\n", "pyplot.plot(x, u_analytical, label='Analytical')\n", "pyplot.xlim([0,2*numpy.pi])\n", "pyplot.ylim([0,10])\n", "pyplot.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "\n", "What next?\n", "----\n", "\n", "The subsequent steps, from 5 to 12, will be in two dimensions. But it is easy to extend the 1D finite-difference formulas to the partial derivatives in 2D or 3D. Just apply the definition — a partial derivative with respect to $x$ is the variation in the $x$ direction *while keeping $y$ constant*.\n", "\n", "Before moving on to [Step 5](http://nbviewer.ipython.org/urls/github.com/barbagroup/CFDPython/blob/master/lessons/07_Step_5.ipynb), make sure you have completed your own code for steps 1 through 4 and you have experimented with the parameters and thought about what is happening. Also, we recommend that you take a slight break to learn about [array operations with NumPy](http://nbviewer.ipython.org/urls/github.com/barbagroup/CFDPython/blob/master/lessons/07_Step_5.ipynb)." ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.core.display import HTML\n", "def css_styling():\n", " styles = open(\"../styles/custom.css\", \"r\").read()\n", " return HTML(styles)\n", "css_styling()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.4.3" } }, "nbformat": 4, "nbformat_minor": 0 }