{ "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": [ "## Burgers' Equation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Hi there! We have reached the final lesson 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", "We have learned about the finite-difference solution for the linear and non-linear convection equations and the diffusion equation. It's time to combine all these into one: *Burgers' equation*. The wonders of *code reuse*!\n", "\n", "Before you continue, make sure you have completed the previous lessons of this series, it will make your life easier. You should have written your own versions of the codes in separate, clean Jupyter Notebooks or Python scripts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can read about Burgers' Equation on its [wikipedia page](http://en.wikipedia.org/wiki/Burgers'_equation).\n", "Burgers' equation in one spatial dimension looks like this:\n", "\n", "$$\n", "\\begin{equation}\n", "\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} = \\nu \\frac{\\partial ^2u}{\\partial x^2}\n", "\\end{equation}\n", "$$\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 the previous notebooks of this module. Using forward difference for time, backward difference for space and our 2nd-order method for the second derivatives yields:\n", "\n", "$$\n", "\\begin{equation}\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", "\\end{equation}\n", "$$\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", "$$\n", "\\begin{equation}\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", "\\end{equation}\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", "The initial condition for this problem is going to be:\n", "\n", "$$\n", "\\begin{eqnarray}\n", "u &=& -\\frac{2 \\nu}{\\phi} \\frac{\\partial \\phi}{\\partial x} + 4 \\\\\\\n", "\\phi(t=0) = \\phi_0 &=& \\exp \\bigg(\\frac{-x^2}{4 \\nu} \\bigg) + \\exp \\bigg(\\frac{-(x-2 \\pi)^2}{4 \\nu} \\bigg)\n", "\\end{eqnarray}\n", "$$\n", "\n", "This has an analytical solution, given by:\n", "\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", "\n", "The boundary condition will be:\n", "\n", "$$\n", "\\begin{equation}\n", "u(0) = u(2\\pi)\n", "\\end{equation}\n", "$$\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": {}, "outputs": [], "source": [ "import numpy\n", "import sympy\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": "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": 3, "metadata": {}, "outputs": [], "source": [ "sympy.init_printing()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Start by setting up symbolic variables for the three variables in our initial condition. It's important to recognize that once we've defined these symbolic variables, they function differently than \"regular\" Python variables. \n", "\n", "If we type `x` into a code block, we'll get an error:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "ename": "NameError", "evalue": "name 'x' is not defined", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;31mNameError\u001b[0m: name 'x' is not defined" ] } ], "source": [ "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "`x` is not defined, so this shouldn't be a surprise. Now, let's set up `x` as a *symbolic* variable:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "x = sympy.symbols('x')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's see what happens when we type `x` into a code cell:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAsAAAAJCAYAAADkZNYtAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAwUlEQVQYGVWQ3Q3CMAyE04oBKrEBbFCJDWADJDaADUB9St+QGKGMwAjtBghG6AZI2SB8l8b8WDqf7Vxd20WM0cnatq2gfUqcW8EnsAA1mINXiTNhwwcXgVKfUed8S76b4WQNOKdocktIXbtcO8Cj0xje+0psIO8Fy41tjJA7GK0JNMqfJfFvhRm1kGyY6OvTzAi0wADrD+rqiJ9iGfER6koCPd5AEsEb8DHetWiAQ6HlSK7gnhW6gO6tizxUQ5iu8gZXOlF/Vp9rRgAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle x$" ], "text/plain": [ "x" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The value of `x` is $x$. Sympy is also referred to as a computer algebra system -- normally the value of `5*x` will return the product of `5` and whatever value `x` is pointing to. But, if we define `x` as a symbol, then something else happens:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAABQAAAAOCAYAAAAvxDzwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABZklEQVQ4Ea2U303DQAyH04oBKtigbFCJDWADEBvABqA+JW9Vu0HDBrQbtBtAO0I2ALJB+L7rXSB5ACph6Xf+c7brs50OmqbJiqLYZVm2BM/INZgg34MN8hr+ZxrEhB9EjHpRC5I99my/qifRo4JvwRi8gDXJtB1PVpjn+Ur+HxgeX8LPEenJGU+8i6728hwsse2/h6N7l/wukO2xbXKIZ+AtVaijEy7BAnkOdsiX8EAx2dT76LPhQkyifo18GyrEcBWi4oFeAYfkKlmtNAWzIB0O7VZXRptrVqUKo63DnPKYxAZJM+T6IIbTZ26TzQJANeRYARe7Tyk4JMQn6cnPdvjkDllh26fOzdeiv/bsGcmtTrItHbKHDsP398kf2qfK4DY9PTEUga3dAuQH7ksrnKPY/JbQXY1TcKMR3QQrEBLB+0O0Lf4H1Olb1uBO2Sfld3Ud4CYcwZ6An6XkZP1RJx36j0+Y9icZucR0ZBkKdgAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle 5 x$" ], "text/plain": [ "5⋅x" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "5 * x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This will let us manipulate an equation with unknowns using Python! Let's start by defining symbols for $x$, $\\nu$ and $t$ 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": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOgAAAAfCAYAAADtPoHZAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAGvElEQVR4Ae2b7XHcNhCGTxoVoJE7iDtQkgoidyBbFdjuQBn/kv5l5A5kV5CRO5A7sOMO4g6SqIPL+8AAB8evIwARJChiBkcQwC4W+wEsiL2D7Xa7WVM5HLi+vj4VtReWYspvVfe9nBmslPoc2CfPI7/zWi6CAzcS6gsotcK9V/F5EZSvRLZxoFeeh20Qa92sOXDiUcfO+ZP3vhbL40CvPFcDLUyg2jV/9kg+U3l1bz2GlFbcJ8/sBiqCzqxrlsRL4bhXPg5FIpibUJiY/syRsWyG1qCdTv0vlTvnZ9veiTbj7sbQWCKM5p2sP8KRXXc0ZpQ8sxqoiERJUdxvKcoh+DeCbyg8+JXv9uDuVPo9cKHNGOfvZAGSOSsOToJ7r86ti4namMNH5d9UfjI7qOaarD/CMYnuaNwoeWY1UBTOEqpiXBI8bt1X5YcWDLR9aamfoqr3bDGQIFb6c7+v3jFODPc19XpvNWLaFpiS9Ee8mlp3guV5lEuIVrH+TRnP4jjW87PyDiq9o6goM21vlD/sdBjwIhhWaHA8U/5HmUXghepf6hmUBNN5VlQbinKrjBHDE/f85o+l8idlPIJPyi79pQJ0shOQ2EHZoRedxAcWpmj9sfCT6k6MPLMZqJj7i/LfiVqE4eEqNJLqcSfPld/6jXrHcBGuS5xhMA6XWNWcAeB+v1fGgMCHcSZdYQiesetnRRTlubJzgzmftM5LsD7tG+Ac4U/smao/c9GdIHnmNFBW/cZ5SQoHwX1umjEg9TtVv2d6XlrFBN87vf+pzM4DngfbVj1Uv7O76P1WeceIXWfVV4aqOnNm9OFVHkSrh4/+7qxY0cY4FhftpD6j813lH72f5m9Df4bKQ/0m1x1PZEHyzGmgGCeM2kliHorbajB+R/XjwxLZJL2zK/1h4alj1+NsulEdq2WMi2t2NKHgy6gxbB+XyoNotTRgfCw81VlR8P5i8Upt7sMRtHelaLeuC2Gh9Q39GSoP9ZtcdzyeB8nz0AMcu4jx/Jo6iJiNe+jcVgyKlZVkjFfv7LCfTU34D7g522E4FxaXMfpwVBvOiuD6z2bOtn5iEXB04gE0Fi/bmUVhTT8W3yT9EY/noDth8iQWN1e+urq6yzVW1zii4barbW71ovWcPDe6pqJnav1J1Z0YeebcQdkF+PDizpBT7Qr+B6KpaBg6Lh+p3Ll4KMyS+02tP6m6EyzPA1bDnEkKZ85bejr3LufwxYxlF7IPeoa5RMXMMI7QUvUnVp5BBmqZ43/o6OLyS1+xVM67CnRRtfB68flgzlNc9SdCOlOdJx5zXPn298rHoTgFc6Z8GgF3sw9GeO/qNOn9sl63D8/aPvwbiXi7OD3IfQaNWEL6QbQqt8ZW9kOZqxi+/jbigoUvKZ5X8CZIXrgbVydq64zH3Efv2O2ijSCPqb8PRE9TtC9SD4o2UAnF3X3GnNO4ommL3gFndDyvcHJlwjGg676rEY8ZrZUroOHAkvWgWAOVUI4lHe61zP2nr6uqYzfY2j5OiC4oYGPrGwakeu5XMS7C8FiRHz0JL19lLx4d8RNFaGW5WD048uVqJ0uEDoHiJELQ2GkaIXqmddofooXadsCN6gmleyBDop71IIDWuE714zN+ajzvEK6wuCwuiXfMK7f+LFoPKgMVc3HtuOfhrsYYpJ6c0wjDG/LlVt3yJNGFwfXF5dLuX+MQD+v/I4V5NRYd4UXBjFHrWSXV78xf753xvBVQf+Gkv7m8VvEku/5ozMXrgTFQTRTF5G9Nr1V2xkkdipl6OSsUj5tEI25t5drqnVW7Hpf7xc4Lt5VrH98gKSPcekLJTGif+rMyB8fz1hF2vDfc645+RVRbPmfXH427eD0wBiotYIchnWjS7kseOwkuX2NHMT1n8GMVA+NkMcEVd+44MZsYIX8tawvExwjb6hE4HgQ8iIrgESyGj6GzS0OP/3c2VZk0GU9FDwsW9NUTPET+befj76r3PZA67KT6I9qgfZl6wD2b7o+2ypPHyea880uZr2Cj43kFO8v4WkvXZYwMBLtVLlJ/UugW7Oh6cOgthb4L6FUvtpgS15ni9gfHYxYigVL1Z9Z64Ay0k7lyH5zLW4ieDCNT82LO3Fm2uXu9SARTnX97O9YaLS93PjjVupT6Wqz+zF0P3BmUHWHn7CHCzblO9ZxZFpk0R/9L7+hz1Hit10KjDzz+AEXrz5z1oAqWF5HslNx78kdjk1Q31ldMN8T6nAkHJGv+UE6YY9QisurPOIJ0O+gmVjDjkLViLY0Dq/6MI7HKQMdBv2ItiANZ3f2C+DIpqf8DY+mXBUAJC10AAAAASUVORK5CYII=\n", "text/latex": [ "$\\displaystyle e^{- \\frac{\\left(- 4 t + x - 2 \\pi\\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 - 2⋅π) -(-4⋅t + x) \n", " ─────────────────── ─────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) \n", "ℯ + ℯ " ] }, "execution_count": 8, "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))) +\n", " sympy.exp(-(x - 4 * t - 2 * sympy.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. To take a derivative with respect to $x$, we can just use:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfwAAAA+CAYAAADOBakYAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAT6klEQVR4Ae2d65EVtxaFDxQBUDgC2xlgiIBxBthEAGSAi1/4HwUZ2ETgizPAjgBwBtgRGE8Gc9fXIzX9fut0nz5LVRp167G1taS9t16n58bV1dXBzggYgToCP//8813FPgopPD9V3N/1nI4xAkbACKyDwBg9dWsdFl2rETgJBF5JmL6H0yBU7/T47UlwbiaNgBE4FwQG66mb54KI22kEJiBwp1CGlf03hXc/GgEjYAS2gMBgPWWDv4XuMg+bRECr+u8KjF3o2dv5BUD8aASMwPoIjNFTNvjr99duONDAuwhb37PaJBrv5G+PJaIyr8aWiflV9llbnSH+ufJm2/uxjEMjYASWQUAydrK6YwwCauddebbg8ei5b0aWn6WnbPDHoO28rQiEgctg/qs104AElX+ibDUhgL782x4SoycJkZ5ov9ZzbcKgeGi+kX+gZ6/wI2AOjcBCCEiukPeT1R0jYcDQ/4RXOTz3ggY7lZulp2zwB0PtjD0IMJAZjJOdyrNt/kH+soEIae8b4peMYsb9MBLUM8aeScBj4vRemxAQb2cEjMAsBPagO4YCMPi8vYPgZD3lW/odqDppGALBMH4elrs5V6BxW+Ef8qVMesfQYohJeyL/aynDgBeVYRUBja/k/5VnUvG94n9QmDk9/y7PLsLvIeqjQsqx64Bjhc+s3M4IGIEFEJC8Mak+ed0xFAq1t/VekNJY1Pwiz6QATGL4l9IW0VM2+ELVbjYC90Th00wqGPLGHQLFswX2UP5psQ69MxFAYUTHOSACEx0z4Wi82TJ8LY9QQQ9j3/QTu5xeS3qk7dAIGIH5COxJdwxGQ7oFPVO9F8SC51v5uO3PeX2jTlTZSXrKBn9wFzljBwKsgmvn22FQd22DZwZZ+e6q/FcKn4U6oPdc77/JM7tlcF+GtDxQfGm1rfdf5EuTgphZ8bnhV1x2blYtH/IWt9xicYdGwAikQeDkdYf0CPqpV89F+EL+eC8o12voqJAGPVzTguQ65Xr1H58Hhzb4g6Fyxg4EMPYY7ZLT4GUwNxrgYkbl46IfPnN6Z+b7MpQnjlU5Z/sHxbETMGVLP5s1iwQ37bOJQgutWduL8GhnBIzAYAROXndIjwzScyCivBhzJgf5vSDFFRcuPyotXuRD77W5SXrqZhs1xxuBEQhgjO+PyN+YFWGQj9v0GGhm/7hsMqB3dgD+yGLG/4E2Z/EI06NAK5tEVEghvHZGwAgcB4E96Y4hiHEvCD30X/DcKyo6FiRRx7G7WVtIhczT9BTf0rc3BnPHwIsXL97OpTG3vHj4ZQ4NlX+In0PDZS1LHgPjxsAedMcx+3yOnvIKvzi38vMcBLgIF8/g59CZU7Z4YW8KHS7y/T6loMsYASMwGYE96I7JjZ9QcLKeSmrwpTz5WEr8SdOEdrlIFQHhyW31uNVdTV7tXTxxFscWVNe5U1L+VHd+D2BsRSrLZKV4ljaWxGbyqy2Wu830xnYY0biw7mjpjjm6o4Vkkui5euoGWxEpnBjjcsJbhZxJLOYC3TiJ4DfVsZ547rFYXXMJBV7j7c1opJnNTjZM8KTy/FYcOrWb8VN5Fi0M9RCD94Py5udHek4zgKY25Ijl1PYbR6xuUFXiKcrD2codQAUcFpe9QZ2wcCa1hXNc/jVz7wXYvqpFw7qjD6SdpavPv+ipVGcPOmf4KP/N0vRFs3ZOq7i38o1nr4r/JL84H33tUp235Uvn2np/JX8lf9FXvitd5aH9sSvP1tLE7zv4HsuXylzI351Q7lVfGdFl3JR40vuzalwfnS2li/ezljv6gv6TTyJ7a/S12oIOK7VnKh8Bm5PSHTPaap1TuaOXZEtfMwpW4B8ULrYCZdIletxo5JZj1fETB37KVXIhP9ubk/igvPzUc+n8pxeRKdFiBc3qmFn2ZCc60OCrc1N5m1z3lILik/EQdzgGk1A5ytS+sU28fB+GrHQbncpCk/6pHT8ovvFb1Y2ENhYp3nchd8CqtmxS9o7d5cIBGV/s2xCid1K6Yyreaqd1TgN4SQy+6sGwxe20hmonR/HTr/wTgxUqTQqebc1Z2+eVOsa8/qjMTZMTjh74idhoA1ip/KXea5OcSp7VX9VOjCo/vUHRjHX8NK/pS1PQfD+WWMwvmtw1YIx+jnGVsPSt6krall8td9e9k1r2jjIGNEbREcjNpAVLB5MnoTs6+O9MEm7WOS0ILW7wwyA9KFx6kNIEfkN9Idoo5KKBx/A1TTDo+LXO9mk/K9Ein4rKXVt8nqHrQXQzRaCQNm7ShbYzualNuhTHCu4q5Mn413P84MQhxNcMsuLpZwwbn6BkFr+4E11u6j9anHBCguI5m0AqPHe5A+WkspewG6ukObcf/ZGpKpHqu2huXndUeR76rrahV61zWgC71RI/J5pt9yRGVp3JNjbKmDr+0zOXWPj84PsQf1CIAWQXAAWIZ/uWn2t9Uti0WlTS8k51Ff9JQrGC+CGFXDErL4M0Gi92MTBoGe8Ks3/20sI7ONPWXrxDHUyM/pXHgRsr6JyPLHbZP09a+D4ons9IXuKpUmHEJXJwTw+17/MrH5cVmSyULjDpnYkAOEbHxJB+j674Xf0Y1xUWaXXl20qa5S70hPr9GLKXtN/VBvRBcfzW6lMetvuR6a6x+qvylWQlENq67qi1d2CEdU4HUCkMPtvofeerHSx1J2nwcksc5c5gRyBYPbJFlTmlMZCZGKAAUfqL3la+rmXaX/HCZARDzj9xiYYOYeW78Rj5g0IEnZUu34UnH0aPvE2TFdKahFnRX5xoUC9Y8fvNzMArhA/KZvUqXNSJPga86/v4pBcnKm/0Xjyugb/aZER0wSvDTmHuFF9qh95bv6ufF+p+uNOdvLlUy11Hl2g8LC17HbXNSwpjnFVqbfxHykqLO5oPFIes4FkM0U5khOdDB43N6g74nuLUVuucHuBu9aRPSWbg1bZipxBqKqNOxZAzoFFwGDE6+R/FMxEoGhDSa1vJilvTwS8r26Jxeq64fMKiZ1beYBi38jDKbYIPzmDR6lQX6UzAHus5Gnvi4AF+kjjVBfY5/nrP2qnwMlSIYmJnBl5QXvRfsZ0807dVR7kPRCr/E/mIUzXf3HewPSVnuevuraVlr7u2ean5AqCJjMY8clHc1Xykdya4f8sjTyx4irLURGazuqOJ2SFxarN1Tg9QKQz+HdUZlXqp+jAY/yxF9r9gqDLDoZDVL6vUuBLkHJeVPgaD3/x/LR/rxjBks1yFnU5lKE/+qkN47igdgao6hCvyUU2rvSsvCgdBrK7Ii/8khnIIM/mydigsTmJILzqEGkXf5Vg542gHWOGgzdZ4xCqLTPFHdYAhxp6QI4R4jHBf7/DPLkwVE0VnRr0pnrHAGKAtg/oXYkWnsmBMf4Md/DRt9yfHpsjTAs8nJ3e0WdifquxlXSb+GdeTdVpGpPBH9FjQ/FaIqj0qT8mwKQPHltHAI1dDJsGb1x21hg+MCH1indOAVwqD31DNdZQ6AiXadr7WWq6QgHL4uvB+EE22vTGK3IjPjLzeEUKUeX4JTM+tTvmLK+48n+IRPi7eNW2n5/n6HlQ+M7QKawZMcVXDQhsa+emrpyWdNrCrMEQJtJCYHh3aR3tKbVJ852SJcvK1ihWHoqrhWMvYESEaUWE29qvSe5VuB/nNJak9m5Q7gBJvpXERwQt9sHnZWwDb2GSwQG/dVzh4Iqu8TF6LOgTjX3zP6U94WFV3TOA3KxLab53TAODNhri5UZ9FgIG7qFMnQpNzrdpgVhwKHCFhlYO7xx/F56tjPS/OE3X0OdWL0LATkRspPaPImJCUnOIQXlzO9/Vr61/aG2f2rZkG5ukqv1YauxBxV2IsD+yoTHXsIAxWulMrWbic5a4CqPowpexValvklck+uoHt+dwrDr3AThhxVXlArxQXNpQfouv2rjsEyyS3a52TwuBjgBhMizoNYgw9q76aoQwVMcijoSyd36sMgt9WLhRfPlC9CCoz9tzYh1rgJ3PwJh8FFIE/6J0JTOb0/KyQHqNjSJvApcu1Tgig3VVw7TTxB++Tvs+vsjmGY9oRMGlcdY6hs0Jey10BdPVjatkr1LbMo3hmJ467LNlndGMo6sh4dhyouHxXSs/oDY45/ydfdPeKLy3Pu9YdLW3ujRamu9Y5t3oRGJ8BRTtn276rRraBOat/IJ8bOj1nN9tDZ1EeQcgMneJ4ZtZ71BWb6kOguCyHEHMUgYMXJkNsu3EUgYEnD+2CPyYquQs0mOTkbc0Trx/AOU5yKkn5KyvdR/mbHkQPPuAp8lVM3tSzeO1r36L8qr5coS5KOD0xy13AWH14DNlL36PlGpDZqkPvcZeoqB94Ro/0yc3udUcVrKHve9Y5KQw+20tztlNb+4WOkH+sDG8UflbI4EYQmAQUBziGjO0vBOKgcA0ljiFH8TStorPJiNI+yGPomZCQD8P/RM/gl32lT89dZ+/M5MGj1dF2+YN8TpPMeq/uOrTScMJJIGC5+9JNx5C9L7UlepKMosfYqUDHXeiddnHBNOoEDHt1N6qoB5Xc6qw7WqHZb0KS/5anAclvPDkHjYbtJBEU/2y9Y4zXmDB0YiaeUAIfFfIzPjsjcNBY2IXc0ZVqy2Zl79SHmrC17jj1TpzI/82J5fqKMTOtzjz7ymwxndlynE1vjb/nYijJTsrWGmp+BiOwF7mjwVuWvcEdstGM1h0b7ZjUbCVZ4cO0ZpFsMXL55KRX+ak7YAr9MEP/U2GquxJT2HKZDSBgudtAJ2yYBeuODXfOEVhLtcKHdc6jvQJN04nxol8a6qZ6yghY7k6599Lzbt2RHuPN1pBshU+LNZvk0hqXTba6Lb7ZjmljTFhyEZHLi945aQPpzOMtd2c+AFqab93RAswZRSc1+GeEo5tqBIyAETACRmDTCKTc0t90w82cETACRsAIGIFzQsAG/5x62201AkbACBiBs0XgxosXL/hN5uz/9qTzoauzRdENNwJCQDJwow0IpS0iZ1X6lrsqIn43AkagDQGf4bch43gjYASMgBEwAjtCwFv6O+pMN8UIGAEjYASMQBsCNvhtyDjeCBgBI2AEjMCOELDB31FnuilGwAgYASNgBNoQsMFvQ2bleF3G4p/2ZP/tb2VWRlcvvh/C/+iCLmAEjEAjAtYHjbA4ciQCNvgjASO7hI//Z5/ss8GizY1u/r3vSX6hUHzzL39fKbTRZ8DY7R4BjfVkOkG0rQ92P4KO08Bbx6lmd7XwPeq/EraKn0nyTfRGJwWw+r9BFQ/8+9I38l/r+bKB0ceKox3+Bz8N4Dhqdwik1AnWB7sbLus0yCv8kbjLuD1TkTsjiw3OLvps439Q2PitfMXH/xPemN5XEeXlacNop3K35d/K829YH8mz8mh0ysMkgG/+T6qrkagjjcAGEQhjPIlOEG3rgw32+amy5BX+iJ6T8LFFjSGbZGwHVvWT8n3fkZe0lLsLrVWr/bQ923kIiojJR5d7qcR/5F93ZXKaEThVBI6gE6wPTnVwbJBvr/DHdcpTCXiyc/WgPA4KuyYUF2L5j3Fsr5Nb7cgmRwrh2c4I7BGBZDpBcpPdgbE+2OOwWadNXuEPxF1Cx9Za60W9YNTY6kZIf9U7M/PchfI/KOxavbNirhnzQJuVNbTx8YLQJ6VtffVMe+C91i7F2RmBk0UgyHRKnWB9cLKjY5uM2+AP6BcJ9m1l4/y6deWtNM6rHyjff/LvG8gyGejbHWAywOWfkoO2IqCPArhQ2DVpKJXdwAsXDJ9ugA+zYAQWQ0AyeAydYH2wWI+ZEAjcNAyDEHguAR+yko5b16XVbDDUKIjW1UDggtX75w6OUACrnN938NSXRHtou50R2BMCx9AJ1gd7GjEbaItX+D2dEIz1bz3ZYjIG+W+VuYwRIXyukG3+1h2CkI+bvtWyISkLmFDwG/dep7rYUYgTkGJ+jO8dpXPLvurgvfXngNXMA99pM4rLzgjsAgHJCDttx9AJ1ge7GDHbaYQNfkdfSLAxjvcVDjKyyouBra7u+VkadEpn+nof5QIvGM53Qwoqf2N9ikdZ8RW/ITsWQ6pyHiNwNggEOVxdJwQ+rA/OZuQt01Ab/G4cMeAYx+pW/N1CfHZxrkkAFUc+VvdtH6dRUsl1bX/fI6do5hMK6pTv2hEoEV/phVVK387GSqy5WiMwGoFj6gTrg9Hd4wJdCNjgd6AjY8rKvra6V/yPiucSXfEyGooAlxlkpbGSxth/p+ehRhnDiIFscqXz+0Cf/Fs/02cVMrT9Te12nBHYDAKSu2PqBOuDzfT8PhixwZ/ej2zTF112fk+ElALn54RjPyuL8W4rQ33ZSll0eWbnoTYZUfyxXLX9bfXSnnxXoi2T443ADhCoysRcnWB9sINBsaUm2OCP6A0ZWAw52/QINj+P4yd07xTyczu23BFQVv8vFXepcKzjfL56fBBpUDf/UCf7D3oKVzmDV73wxy5E3NH4U3Ef9P5RYdPPDsHlsbydEdgdAhrzKXWC9cHuRsy6DbpxdXW1LgeuvYSAFEjSf4wj+hw1HOXSnupiYsRE4NtSI/1iBIzAIASsDwbB5EwDEbg5MJ+zHQ8BVgyNN+wXYoHt9aaV+ELkS2S4w9C2Y1HK6BcjYAQaEbA+aITFkVMQ8Ap/CmqJy2hWz1Ye3+jOzuwTV5eEvHhndc92f9udhCT1mqgR2BsC1gd769H12uMV/nrYd9XMx29OfWXM/YalP+LThZnTjMBeEbA+2GvPHrldXuEfGfCh1WlWz8/ZuBh4rO33oaz15hPPXCzkZ4snu0PR20hnMAJHRMD64Ihg77iq/wNrxyWGkPw9qAAAAABJRU5ErkJggg==\n", "text/latex": [ "$\\displaystyle - \\frac{\\left(- 8 t + 2 x\\right) e^{- \\frac{\\left(- 4 t + x\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}}{4 \\nu \\left(t + 1\\right)} - \\frac{\\left(- 8 t + 2 x - 4 \\pi\\right) e^{- \\frac{\\left(- 4 t + x - 2 \\pi\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}}{4 \\nu \\left(t + 1\\right)}$" ], "text/plain": [ " 2 2 \n", " -(-4⋅t + x) -(-4⋅t + x - 2⋅π) \n", " ───────────── ───────────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) \n", " (-8⋅t + 2⋅x)⋅ℯ (-8⋅t + 2⋅x - 4⋅π)⋅ℯ \n", "- ─────────────────────────── - ───────────────────────────────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) " ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phiprime = phi.diff(x)\n", "phiprime" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you want to see the non-rendered version, just use the Python print command." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "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 - 4*pi)*exp(-(-4*t + x - 2*pi)**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": 11, "metadata": {}, "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 - 4*pi)*exp(-(-4*t + x - 2*pi)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)))/(exp(-(-4*t + x - 2*pi)**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 usable function, we tell lambdify which variables to request and the function we want to plug them into." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The value of u at t=1, x=4, nu=3 is 3.49170664206445\n" ] } ], "source": [ "u_lamb = lambdify((t, x, nu), u)\n", "print('The value of u at t=1, x=4, nu=3 is {}'.format(u_lamb(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": 13, "metadata": {}, "outputs": [], "source": [ "# Set parameters.\n", "nx = 101 # number of spatial grid points\n", "L = 2.0 * numpy.pi # length of the domain\n", "dx = L / (nx - 1) # spatial grid size\n", "nu = 0.07 # viscosity\n", "nt = 100 # number of time steps to compute\n", "sigma = 0.1 # CFL limit\n", "dt = sigma * dx**2 / nu # time-step size\n", "\n", "# Discretize the domain.\n", "x = numpy.linspace(0.0, L, num=nx)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have a function `u_lamb` but we need to create an array `u0` with our initial conditions. `u_lamb` will return the value for any given time $t$, position $x$ and $nu$. We can use a `for`-loop to cycle through values of `x` to generate the `u0` array. That code would look something like this:\n", "\n", "```Python\n", "u0 = numpy.empty(nx)\n", "\n", "for i, x0 in enumerate(x):\n", " u0[i] = u_lamb(t, x0, nu)\n", "```\n", "\n", "But there's a cleaner, more beautiful way to do this -- *list comprehension*. \n", "\n", "We can create a list of all of the appropriate `u` values by typing\n", "\n", "```Python\n", "[u_lamb(t, x0, nu) for x0 in x]\n", "```\n", "\n", "You can see that the syntax is similar to the `for`-loop, but it only takes one line. Using a list comprehension will create... a list. This is different from an *array*, but converting a list to an array is trivial using `numpy.asarray()`. \n", "\n", "With the list comprehension in place, the three lines of code above become one:\n", "\n", "```Python\n", "u = numpy.asarray([u_lamb(t, x0, nu) for x0 in x])\n", "```" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "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,\n", " 4. ])" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set initial conditions.\n", "t = 0.0\n", "u0 = numpy.array([u_lamb(t, xi, nu) for xi in x])\n", "u0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now that we have the initial conditions set up, we can plot it to see what $u(x,0)$ looks like:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEoCAYAAACkdq2MAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxX9Z3v8dcnCSFAwh72JO6gIqAgSa1atGqtrftO0t6209o7nZkud9pHp+3UttO9ffTOTDud22v3O4ni2lqtjnvcEzYRFRAUSQIhhIAsAbJ/7h/nJIT4A35Zz+/3y/v5eOTxg/PbPslJ8s453/P9fszdERERiUda1AWIiEjyUGiIiEjcFBoiIhI3hYaIiMRNoSEiInFTaIiISNwUGtJrZrbWzHaamZtZi5nVmdk3+vmaT5rZG2Y2spfP+62ZbTWzqX1834+H9beYWcpcfx5+TnvDffSJHvddFt73yV6+5llm1mBm3xrQYiWpKDSk19x9HnBu+N+X3H2au3+/ny87GRgPpHffaGbfDn/xLTnK8yYBY4HMvrypu/8/d58GvNSX5yeq8HP6wlHuzgZygHE97wi/1uVHed6o8HkTBqJGSU4ZURcgEloIpLl7ay+fdy2Q6e7Ng1BTSnL3B8xstLs39fJ5y81sXG+fJ6lFoSEJwd3bgfY+PM8BBUYv9fUXvwJDdHpKBpSZ/Ud4Pt3NrNzMFoW3O82sysx+ZmZZ3R4/O3x8Y/fTUGaWbmZ1wJfDhz4QPq7OzL5hZqOPc97+Q2Z2l5m9bWY7zGyXmf3VzBYP4Od6iZk9bmb1YS0bzOwPZnZJj8eZmf1PM1sVPrbezF4wsxt7PO6r4eu0m9kWMzvFzB4Ot203s9+b2fgYdWSY2bfC5+w2s/Vm9hXAYjz21933T8/3Dv97XrevdZ2ZjTSz27vXFuN1M83sa+G41I7wc3zczC7u8bhefX90e16xmVWaWa2ZbTOzV8zsJ2Z22rH3kgw4d9eHPnr9AZwAOFB+lPsdeBsoAyYS/AL7dLj9xzEe/+3wviXxbO92/yfC+z/RY3s58AIwM/x/LnA/0AScE+N1ygkPXOL8/G8DOoCfAKPCbQuAt4A9PR77R6AFuCX8OmQAXwzr/k6M194CNAD3AbPCbR8mOBK7O8bjS8NaPkPwh+Ao4GvA67G+Nt32z3v23bH2abfatvTYNgJ4EtgDfDDcNgr4aVjX3/Tn+wMoDl/num7bLgL2At+O+mdhuH1EXoA+kvMjztBoAab12F4LbI7x+IEOjf8LLOyxbSzQBiyL8TpxhwYwMwyfTYD1uO+q7qEBXBfW97sYr/NM+MuwZ51bwuf03P4Swam4Ed22XRQ+9pEYr//sEIXG/wqfd3uP7UYQok3AjL5+fwAPALti1PId4PNR/ywMtw+dnpLBtNnd63psqwJmDfYbu/tn3X1Vj237CH4pzevny98IjAT+7OFvr24eA67v9v+PhbcPxHid+wl+sRbHuK+pZ/0EX7tMgqOmTp3v9VCM13gyxrbBEPNzDL82fyb4Wt0Q43nxfn/sACaa2X+aWdd97v4td/95vyqXXlNoyGCqj7GtmeB0xqAyszwz+1cze7X7+XmCo4TR/Xz5zvPo23re4e7N7v5Ut02zw9utMV6n8/lzYty3M8a2zgH/7pcXnxLe1sZ4fKxtg6Gvn2O83x+3A08AfwtUh2MbXzOz6X0pVvpHoSGDqSOKNw1/mawGbgI+RzCuMc2DuQs1A/EWA/zYWJMKe/u1601NUejz5+juO939MuAs4HsE80R+AGwys6sGrkSJh0JDUtFNBJMFf+HuL3pwOe9A2hDezux5R3jV1ywzyzzeY7tte7MftbwV3s6IcV+sbYOhs/5B+RzDr6m5++vufru7n0YwP2cE8Iu+vq70jUJDEl1jeDsCwMzmmNlPjvOczrkER/wla8ESJdMGoKZ7CU6jXGNmPf/CvwHYSHD1EARXNkHwS66nzm139qOWznGEq2Pcd2kfXu8A3U4PmdkXe142G0PMzzH82lxN8LW6rw+1dHqK4A+BLu7+Z4KrwzQ7fYgpNCTRvRbezg9vPw588DjPeYTg8s/Pm9nZAGaWA/yKYFC2X9y9Fvg8cDLww855BeEckJ8BP3D3veFj7ye4rPRjZnZTOGcjw8z+AbgY+L67r+xHLU8Dy4BLzew2M0szsywz+xpwah9e8jXgFDMbY8F6Xt8CxhznOT8nuPrsHzsDJvya/JDga/T58GvWH183s86xE8zso8Bc4L/6+brSW1FfvqWP5PsA1hIM1HZeNlkHfCO87/bw/93vOy/8qAu3eedzCAZR6wiOKBzYDazo8X4/IxjUrQOWA4sIBrPrCK7V9/C2rttzzgEeJZjv0AC8SnB5bhXBfIc6ggD6eIy6/j3Or8MlBAO0O8PnvQJ8OsbjjGAQdzXB4O9O4EXg5h6P+2T4Ou3darwZyAv/fSiscSfwq27PG0FwaXJV+PV7C/g+wVySrq8NQWD+Osb+ubnbay0Mv8YNBOM/Pwvr79yv3Wv7crfnjQS+Dqzr9jk+AVza43Ps1fdH+Jzzw7rXh98H24FVwN8D6VH/PAy3Dwt3SiTCAcvfAx9y90QfyBMRGfYiOz1lZtcCLxMcvh7rcSPM7LvhEg2vm9lLZnb+0FQpIiLdRTmm8U8EA3UvHudxvyA4RL/A3ecCvwOeMLMFg1yfiIj0EGVovN/dNx3rAeHA123Aj9x9J4C7/wbYTHDOVkREhlBkoeHubXE87FqCQbhnemx/GrjMzLIHvDARETmqRO+nMY/gWvvqHtvfIaj9DIIrPY5gZrcRHKGQlZW1MD8/f5DLlP7q6OggLU1XgCcD7avk0J/9tHHjxgZ3z411X6KHxmTgoL93Ru++8HZSrCe5+x3AHQCzZ8/2N9/sz4RbGQrl5eUsWbIk6jIkDtpXyaE/+8nMqo52X7L+uaDLc0VEIpDoodEAjDaz9B7bc8LbXUNcj4jIsJboobGWoMa8HttPJGims37IKxIRGcYSPTT+RNi1rcf2i4DH3X3/kFckIjKMJXRouPubBAPaXzOzyQBm9imCWeTfiLI2EZHhKLKrp8zspwQzwvPD/68J71rs7i3dHvoPBCttvmhmrcB+4DJ3X4OIiAypyELD3b8S5+NagX8OP0REJEIJfXpKREQSi0JDRETiptAQEZG4KTRERCRuCg0REYmbQkNEROKm0BARkbgpNEREJG4KDRERiZtCQ0RE4qbQEBGRuCk0REQkbgoNERGJm0JDRETiptAQEZG4KTRERCRuCg0REYmbQkNEROKm0BARkbgpNEREJG4KDRERiZtCQ0RE4qbQEBGRuCk0REQkbgoNERGJm0JDRETiptAQEZG4KTRERCRuCg0REYmbQkNEROKW8KFhZovM7FEzW29mr5nZcjO7Meq6RESGo4QODTM7AXgKaADOcvezgN8B95jZlRGWJiIyLCV0aABXAGOB/+3ubQDu/itgH7A0ysJERIajRA+NtvA2o3ODmRlB3emRVCQiMowlemgsAzYA/2xm2WaWBnwdGAn8KtLKRESGIXP3qGs4JjObAfwe+ADQCOwFPuXuzx7jObcBtwHk5uYuvOeee4aiVOmHxsZGsrOzoy5D4qB9lRz6s58uuuiiVe6+KNZ9CR0aZjabYCD8EeCLQBNwE/BLoMTdHz3ea8yePdvffPPNQa1T+q+8vJwlS5ZEXYbEQfsqOfRnP5nZUUMj0U9PfRcYD3zB3Q+6e4e7LwOeA/5oZhnHfrqIiAykRA+Ns4Ct7n6ox/aNQC5w4tCXJCIyfCV6aNQD02McURQADrw79CWJiAxfiR4avyCYp/Ev4aW2mNlFwHXA3e7eEGVxIiLDTUKPCbj7fWZ2OfBPwDozawc6gG8AP4+0OBGRYSihQwPA3R8DHou6DhERSfzTUyIikkAUGiIiEjeFhoiIxE2hISIicVNoiIhI3BQaIiISN4WGiIjETaEhIiJxU2iIiEjcFBoiIhI3hYaIiMRNoSEiInFTaIgMkh37mti0Y3/UZYgMqIRf5VYkWXR0OBWbd/HUhnpe2NTAm2Fg/O4Ti7h4ztSIqxMZGAoNkX6q39/Efau2smx5DdW7D3Ztz0gz2jqc7z68nvNPySUzQwf2kvwUGiJ94O68/PYuyiqreeyNOto6HICZ40dx1YIZXHhqLvNmjePK/3iBzTsP8MeXtvCZC0+KuGqR/lNoiPTCnoMt3LdqK3dWVrO54QAAaQaXnjGVpYvzufC0XNLTrOvx3/zIGXzyDyv4+VObuPacmUzOHhlV6SIDQqEhchzuzpqaPZRWVPPw2lqa2zoAmDp2JLecm88ti/OYPm5UzOdeNGcKHzgtl2c37uRnj2/kh9edNZSliww4hYbIURxobuPBNbWUVVbxRu2+ru0XnDqZ4sICLjl9Chnpxx+n+OZHT+eFf2vg7hXVfKyogDNmjB3MskUGlUJDpIcNdfsoq6jmT69so7G5DYAJo0dw46I8li7O54TJY3r1eqdMyaGkMJ8/vlxFaWUVP7hWRxuSvBQaIkBTazuPvr6dsopqVla927V9UcEESooKuHzuNLJGpPf59YtOmsQfX65iV2PzQJQrEhmFhgxrWxoOcOfyau5dWcO7B1sByB6ZwXXnzGRpYT5zpg3MqaTsrOBHbX9T24C8nkhUFBoy7LS1d/Dk+h2UVVbz/KaGru1nzhhLcWEBVy+YwZiRA/ujkR2+XufpLpFkpdCQYWP73kMsW17DshXV7NgXnCYamZHGlfNnUFJUwPxZ4zCz47xK3+RkjQCgUUcakuQUGpLSOjqc599qoKyiiqc21NMeTsI7OXcMxYUFXH/OLMaNHjHodeSEp6f2KTQkySk0JCXtamzm3nASXufSHhlpxkfmTae4MJ/3nTRp0I4qYukMjcbm1iF7T5HBoNCQlOHurNjyLmWVVTz6Wh0t7cEkvJnjR3Hr4jxuOjePKTlZkdQ2akQ66WlGU2sHre0djIhjfodIIlJoSNLb19TKn1Zvo6yyio07GgEwg4vnTKGkKJ8PnDbliKU9omBmZI/MYO+hVhqb2pgwJjPSekT6SqEhSev1bXsprajiwTW1HGptB2By9khuOTePWxbnMWvC6IgrPFJnaOxXaEgSU2hIUjnU0s5Da2spq6ji1a17u7afd/IkigsLuOzMqQl76qdzXGO/xjUkiSk0JCm8Vb+fsspq7l+1tesKpLFZGdywMI/ionxOzs2OuMLjy9EEP0kBSREaZnY98AVgDDAB2A38u7v/V6SFyaBqaevgv9+oo6yiisp3dndtn583npLCfK6cP6NfS3sMNc3VkFSQ8KFhZl8CPgZc5e5bzWwE8Efgg4BCIwXV7D7YtbRHQ2MLAKMz07l6wUyKC/OZO3NcxBX2TeescJ2ekmSW0KFhZicAPwLOd/etAO7eamZfBmZEWJoMsA53nly3g7LKKso37sSDOXjMmZZDcWE+15w9s+sv9WTVNVdDRxqSxBI6NAiOMPa4+4ruG929FqiNpiQZSPX7mrh7RQ2/f/4Qu5tWApCZntY1CW9hwYQhnYQ3mLI1K1xSQKKHxnnAlnBM44tALsF4xm/c/XdHe5KZ3QbcBpCbm0t5efkQlCrxcnfW7+7g6epWXqlvpz08qpgy2rgobwTnz8wgJ3MPjVv28OyWSEsdUA21wam2dRs3U25bI66m7xobG/UzlQQGaz8lemjkAScAXwauBeqB64G7zGy6u38/1pPc/Q7gDoDZs2f7kiVLhqRYObYj+2s3AZCeZnzojCmcNWovn7vuYtIinoQ3mKoyt3D/pjeYOHUGS5bMjbqcPisvL0c/U4lvsPZToodGFsEVU19x97pw271mdgvwdTP7V3c/GF15cjzuzurqPZRVVvHw2u20hP21p43N4tbF+dx8bh7TxmVRXl6e0oEB3S+51UC4JK9ED4394e2aHttfAa4DzgBWDmlFEpfG5jYeXLON0opq1m8/3F/7wtNyKSnM5+I58fXXTiXqqSGpINFDYwOwAOj526U9vB1ev3WSwPrt+yirrOLPr9R2/XKcOCaTm8L+2vmTEmtpj6GkgXBJBYkeGg8BtwDzgBe6bZ8LHALeiKIoOVJTazuPvLadsspqVnXrr734hIkUF+Vz+dxpjMxInkl4g2WsJvdJCkj00Lib4Kqp75nZR9290cwuAG4A/sXdD0Rb3vD2TsMB7qys4t5VW9kT9tfOGZnBtefMpLiwgNnTciKuMLFocp+kgoQODXdvN7PLgR8Db5hZE9AM/L27/zra6oan1vYOnlq/g9KKal5463B/7bkzx1JSWMCV8we+v3aq0OQ+SQUJ/9Pt7ruBz0Rdx3BXu+cQy1bUsGx5NfX7g/7aWSPSuGr+DIoLC5g3iP21U0V2twUL3V1fL0lKCR8aEp2ODue5TTsprajm6Q07CNtrc8qUbJYuzh+y/tqpYmRGOpkZabS0ddDc1pFUiy2KdFJoyHs0NDZz78qt3Lm8iprdhwAYkW5cceY0SooKKDxxov5K7qOckRnsamthX1OrQkOSkkJDgGAS3vJ3dlNWWc2jr2+nNVzbY+b4USwtzOemRXnk5oyMuMrkl5OVwa4DLTQ2tTFF1wlIElJoDHN7D7Xyp9VbKausZlN90F87zeCS06dQXFjAhaflRt5fO5VkqxGTJDmFxjC1duseyiqq+curh/tr5+aM5OZFedxamM/M8aMirjA15YwM52poVrgkKYXGMHKwpY2HXq2ltKKa17Yd2V+7pKiAS89I3P7aqUJHGpLsFBrDwKYdYX/t1Vu7flmNGzWCGxfO4tbC5OivnSq0aKEkO4VGiursr11aUcXybv21z84fT0lhAR+ZN11X70QgR4sWSpJTaKSYzv7a96yoYdeBw/21rzk76K995ozk7K+dKjpb1ur0lCQrhUYKaGvv4Jk3d1JWWcWzPftrFxVwzYIZSd9fO1V0jmnoSEOS1YCGhpn9p7t/biBfU45uR9hf+67l1WzfG3TCy8xI46NnTae4KJ9z8lOnv3aq0JiGJLtehYaZffw4D7miH7VIHDo6nJfe3kVpRRVPrN9Be7i2x4mTx1BcGCztMWFMZsRVytF0rXSr01OSpHp7pPGHY9zn/ahDjuPdA2F/7eXVvNMQrAifnmZ8eO40igsLOO/kSSnfLjUV5OiSW0lyvQ2N9bz3aGIMMAdYCvxyIIqSQNBf+13KKqp5+LXD/bWnjzvcX3vq2KyIq5Te6Bxb0piGJKvehsbn3b0qxvZ1ZvYosAx4pv9lDW+NzW386ZVtlFVUsaEuaJNuBh84LZeSogIump077Pprp4rDp6c0piHJqVeh4e5PHeO+Q2Y2p/8lDV/ravdRWlnFg69s40BLsLTHpDGZ3Kj+2ilDjZgk2fV2IPzCWJuBCcA1QNNAFDWcNLW289e12ymtrOKV6j1d2xefOJHiQvXXTjWda09pTEOSVW9PT5UTe8DbgK1ASX8LGi4272zkzspq7lvdrb92VgbXnzOL4sJ8Tp2qdbNTUdc8jZY2OjpcFy9I0ultaLwNfLrHtnagHnjb3dsHpKoU1drewRPrdlBWWcWLb+3q2j5v1jiKC/O5cv4MRmdqvmUqS08zRmemc7ClnYOt7V1jHCLJorffsb9092cHpZIUtm3PIZYtr2bZihp2duuvffX8mSwtzGd+3viIK5ShlJOVwcGWdvY3tSo0JOn0diD83warkFTT3uE8tzFY2uPpDfVH9NcuLsznunNmMW6UlvYYjrJHZrCD5mAwXEuBSZLRnzkDrKGxuWtpj63vHu6v/ZG50ykpzGex+msPe51zNfZpMFySkEJjALg7le/sprSiisfeqOvqr503cRRLFxdw46JZTM5Wf20J5GjRQkliCo1+2HuolQfC/tpvHdFfeyrFRfl84NRcXR0j76FFCyWZKTT64NWaPZRVVvGXV2tpag2W9sjNGcmt5+Zx82L115Zj6xz81gQ/SUYKjTgdbGnjL2tqKas8sr/2+0+ZRElhAZeov7bESY2YJJkpNI5j4479lFVU8cDqbewPz0GPHx32116cz0nqry291LX+lMY0JAkpNGJobmvnv1+vo6yimuVbDvfXXlgwgeLCfK44S/21pe80piHJTKHRTfWuoL/2vSsP99ce09Vfu4AzZoyNuEJJBVq0UJJZ0oWGmT0PnA+c6O5b+vt6be0dPL2hnrLKap7bdLi/9unTx1JSlM/VC2Zq1q4MqGwtWihJLKl+G5rZ9QSB0W91e4P+2stW9OivPW86xYUFnJM/XpPwZFBonoYks6QJDTPLBH4IPEIfe5F3dDgvvt1AWUV1zP7aNyycxfjR6q8tg6tzpVsNhEsySprQAP4OWAlspJehsftAC/etquHOymq27DoIHO6vXVJUwPtOUn9tGTpjNRAuSSwpQsPMJgJfAc4DPtGb5+485BT98Kmu/tozuvXXnqL+2hKBzjENDYRLImpr7zjm/UkRGsDtQKm7b+ntOMOBVmdcewcXzc6luLCAi+ZMIV1HFRKhw5fcKjQkcXQf5z2WhA8NMzsFuAk4vRfPuQ24DSBnSh4/uWAUuaMPQv16nq9fP0iVSn80NjZSXl4edRlDwt0x4FBrO089/UzS/REznPZVMotnP3W4s25XB8/UtPJKfXtXC4djSfjQAH4C/Mjd9x73kSF3vwO4A2D27Nl+4xUXD1ZtMkDKy8tZsmRJ1GUMmZxnH2NfUxsLi96fdBdfDLd9layOtZ+OHOcNrh7NSDMunzuV4sICzv/x0V83oUPDzC4A5gI3R12LyEDKyRrBvqY29je1JV1oSHJyd1ZVvUtpRRWPvFZHS3vfxnkTOjSAS4F0YEW3sYxp4e0jZtYCfN3dH4miOJG+0riGDJX9Ta38+ZVtlFVWs6FuPwBm9HmcN6FDw91vJxgE72Jm3wa+BVwxEDPCRaLQtTy65mrIIKna187XHniNB9ds42BLOwCTszO5aVEety7OJ2/i6D69bkKHhkiq0qKFMhiaWtt5eO12SiuqWFPTBARXQhWdNJHiwgI+dOY0MjP618IhaULDzK4AfkCP01PuviDCskT6JDvsqaEjDRkIb+9spKyimvtXb2XvoeAPkVEZcEvhCRQX5nPKlJwBe6+kCY1w3EJjF5ISOk9P7dOYhvRRS1sHT6zbQWlFFS9v3tW1ff6scRQXFjBu31t86INnDvj7Jk1oiKSSsVoeXfpo67sHWba8hmUramhobAZg1Ih0rjl7BksXF3DWrHEAlJe/PSjvr9AQiUBX9z6NaUgc2jucZzfWU1ZRzTNv1ndNwjttajYlRQVcc/ZMxoanPAebQkMkAloeXeKxc38z96wMJuFt23MIgMz0ND581jSKCws494QJQ97CQaEhEoGugXCdnpIe3J2XN++irLKax16voy08rMibOIriwgJuXDiLSdkjI6tPoSESgc4jDQ2ES6e9B1u5b/VWyiqr2LzzAABpBpeeMZWSogIuOGVyQrRwUGiIRCCna3KfxjSGM3dnTc0eyiqreejVWprDFg5Tx47klnODpT1mjB8VcZVHUmiIRCAnS33Ch7MDzW08uKaWssoq3qjd17X9glMnU1yYzwdPn8qI9P5NwhssCg2RCGRrIHxYerNuP2WVVTywelvXvp8wegQ3Lspj6eJ8Tpg8JuIKj0+hIRIBLVg4fDS3tfPoa3WUVVaxYsu7XdsXFUyguCifD8+dTtaI9Agr7B2FhkgEuhYsVGikrC0NB7hreTX3rtrK7gMtQLDfrz17JksL8zl9+tiIK+wbhYZIBLJGpJOZnkZLewdNre1J9ZemHF1bewdPrq+nrLKK5zc1dG0/Y/pYiovyuXrBzK4/GJJVclcvksSyszLYfaCFxuY2hUaSq9vbxF3Lq1m2opod+4KlPUZmpHHl/BkUF+azIG/8kE/CGywKDZGIZI8MQmN/UxuTI5ysJX3T0eG88FYDpRVVPLWhnvZwEt5JuWMoLizghnNmMW700CztMZQUGiIRydGihUlp94EW7l1Zw53Lq6nadRAI+mt/ZN50igvzed9Jk1LmqCIWhYZIRLRoYfJwd1ZWvUtZj/7aM8eP4tbFedx0bh5Tco7fXzsVKDREItI1wU9zNRLWvs7+2hXVvLnjcH/ti+dMobgwnyWze9dfOxUoNEQiorkaiev1bXspq6ziwTW13fprj+Tmc2dxy7l976+dChQaIhE5PKah01OJ4FBLOw+traWssppXa/Z0bX/fSZMoLsrnsjP63187FSg0RCJyeExDRxpRequ+kTsrq7lvVU3XqsNjszK4fuEsigsLOGVKdsQVJhaFhkhEOsc0tP7U0Gtp6+DxdXWUVlRRsXl31/b5eeMpKczno/NmMCpTc2diUWiIRKRz0UINhA+dmt0HuWt5NfesrKGhMVjao7O/dnFhAXNnjou4wsSn0BCJyFgNhA+J9g7nmQ3B0h7lG3fiYX/t2VNzKC7KH9L+2qlAoSESkcOLFmogfDDU72/inhU13LW85oj+2lecNY2SogIWFgx9f+1UoNAQiYgaMQ08d+flt3dRWlnF42/s6OqvXTBpNEsX53PjojwmjsmMuMrkptAQiUjXkYbGNPptz8EW7lu1lTsrq9ncEPTXTk8zPnTmVIoLCzg/QfprpwKFhkhENLmvfzr7a5dWVPPw2vf21751cT7Txg2PpT2GkkJDJCKHQ0NjGr1xoLmNP68JlvZYt71nf+0CLjl9ChkJ2l87FSg0RCIyptvpKXfXoOxxbKjbR2lFFX9+pfaI/to3Lcrj1iTpr50KFBoiERmRnsaoEekcam3nYEt7V4jIYU2t7Tz6+nZKK6pZVXW4v/a5J0yguLCAy+dOUwOrIabvUpEIZWdlcKi1nf1NbQqNbrY0HODO5dXcu7KGdw8Gp++yR2Zw3TkzKS4sYPa0nIgrHL70XSoSoZysDHbub6axuRUY3oO2re0dPBWjv/aZM8ZSUlTAVfNnKFgTQMLvATNbAPwdcA5BvSOAJ4HvuvvOKGsT6a+c8JfgvmF8BdX2vYe4a3kNd/for33V/BkUFxUwf9Y4jfckkIQPDWAZ8AZwobsfMLOZwFPA5WY2390PRVueSN91LVo4zEKjo8N5btNOyiqreWr9DsI5eJycO4alKdxfO30h5UgAAArbSURBVBUkQ2gAfNXdDwC4+zYz+ynwG+AK4P5IKxPph+G2PPquxmbuDSfhVe8+3F/7w3OnUVyU+v21U0EyhMY8d2/psa02vJ0w1MWIDKSuRkzNqTtXw91ZseVdSiuq+O/Xj+yvvbQwnxsXzRo2/bVTQcKHRozAADgNcOC5WM8xs9uA2wByc3MpLy8ftPpkYDQ2Ng7L/bSnITiH/8rrbzL1wOaIq4lPvPvqYKvzUm0bz9S0sq0xOP9kwPzcdC7Ky2BerpFmW1m3aivrBrfkYWmwfqYSPjR6MrN04FPAb919Y6zHuPsdwB0As2fP9iVLlgxdgdIn5eXlDMf9tLp1I09UbWLqrAKWLDkt6nLicrx99drWw/21D7Ue7q99y7l53LI4j1kThm9/7aE0WD9TSRcawDeBNuBLURci0l85KbJo4aGWdh56tZbSyirWbt3btf19J02ipKiAy86cyggt7ZESkio0zOyTwE3AEndvjLoekf7KTvL1p96q309pRTX3r97aNZg/NiuDGxbmsbQwX/21U1DShIaZfQz4R+Bid6+Puh6RgXB4IDx5jjTaOpy/vFpLWUUVle8c7q+9IG88xYX5XDl/hpb2SGFJERpmVgJ8FbjE3evCbR8FZoTjFyJJKZkuua3ZfZA7l1dT9tJB9rW8AsDozHSuXjCT4sJ89dceJhI+NMysGPg1wVjGJd2u4b4A2B5VXSIDIdG793X21y6trOLZHv21S8L+2jnqrz2sJHxoAL8gWJTnpzHu+84Q1yIyoBK1p0b9vibuXlHDXcurqd3bBEBmRhofOWs6Z2Tu4tPXXKBJeMNUwoeGu0+MugaRwZJIYxodHc7Lm3dRFqO/dnFhPjcsDPprl5eXKzCGsYQPDZFUlghjGu8eaOH+1bH7a5cUFfD+k9VfWw5TaIhEaExmBmZwsKWd9g4nfYh+Obs7r9TsobSiiofXbqcl7K89fVwWt5ybz83n5qm/tsSk0BCJUFqakZ2Zwf7mNhqb2gZ9ZdfG5jYeXLON0opq1of9tc3gA6flUlyYz8Vz1F9bjk2hIRKxnKwgNPY3tw5aaKzf3tlfexsHWoKlPSaOyeTGRbMoXlxA/iQt7SHxUWiIRCw7KwP2Dvy4RlNrO39du52yyipWV+/p2r74hIkUF+Vz+dxpjMzQJDzpHYWGSMS6GjEN0BVU7zQcoKyiivtWb2VP2F87J+yvvVT9taWfFBoiEeu8gqrzF3xftLZ38OS6HZRWVvHiW7u6ts+dOZaSwgKuVH9tGSD6LhKJ2OxpOTy7cSelFVVcesbUXj23ds8hli2vZtmKGur3B705skYE/bVLigqYN2v8YJQsw5hCQyRin73wJO5aXs2zG3fyzIZ6Lpoz5ZiP7+yvXVpRzdMbjuyvXVJUwHVnq7+2DB6FhkjEJmWP5AsfPJXv/XU93/3rOs4/dXLM3hMNjc3cszJY2qNm9yEARqQbV5w5jZKiAgpPnKiZ2jLoFBoiCeDj7zshmJG98wD/9XIVnzr/RCA4qnjp7V3ctbyax9fV0doeHFbMmhD2116YR27OyChLl2FGoSGSADIz0vjGR07nb/64kn97ciNjR42gcvMunt/UQN2+YMHANINLTp9KcVE+Hzg1V0t7SCQUGiIJ4uI5U7jg1Mk8v6mBL9/7atf2GeOyuFlLe0iCUGiIJAgz49tXncnnSlczZexIzj9lMhecmsucaTk6qpCEodAQSSAn52bz2JcujLoMkaPSymQiIhI3hYaIiMRNoSEiInFTaIiISNwUGiIiEjeFhoiIxE2hISIicVNoiIhI3BQaIiISN4WGiIjETaEhIiJxU2iIiEjcFBoiIhI3hYaIiMRNoSEiInFL+NAwsylmVmZmb4Yf95nZrKjrEhEZjhI6NMwsE3gCyATOBM4ADgDPmFl2lLWJiAxHCR0awP8A5gFfdfc2d28HvgqcBPxtpJWJiAxDiR4a1wPV7r65c4O71wHrwvtERGQIJXpozAPeibH9HeCsIa5FRGTYy4i6gOOYDKyKsX0fMNrMRrn7oZ53mtltwG3hf5vN7PVBrFEGxmSgIeoiJC7aV8mhP/up4Gh3JHpoHI0d6053vwO4A8DMVrr7oiGpSvpM+yl5aF8lh8HaT4l+eqoByImxPQc4GOsoQ0REBk+ih8Za4IQY208EXhvaUkREJNFD4wGgwMxO6NxgZlOB04H743yNOwa+LBkE2k/JQ/sqOQzKfjJ3H4zXHRDh5L6VwHqgGOgAfgucD5zt7o0RliciMuwk9JGGu7cAlwLtBHMz1gNjgYsVGCIiQy+hjzRERCSxJPSRhogkHjN73sy8+1ijDB8pGRpaGTfxmdkCM/u1ma0ys1fNbJ2Z/dzMcqOuTY7OzK4nGFOUBGVm15vZc+HP1mYzW2lmHxuo10+50NDKuEljGTARuNDd5xOMXV0GvGhmoyKtTGIKf7Z+CDwSdS0Sm5l9CfgGsNTdFwKzgY3ABwfqPVIuNNDKuMnkq+5+AMDdtwE/BU4Froi0KjmavyO4mnFF1IXIe4WnC38EfNbdtwK4eyvwZeA/Bup9UjE0tDJucpjn7m/12FYb3k4Y6mLk2MxsIvAV4OtR1yJH9TFgj7sfEeruXuvuKwfqTVIxNLQybhIIL6fu6TTAgeeGuBw5vtuBUnffEnUhclTnAVvCMY3nzWyDmb1kZp8ayDdJ1gULj6VPK+NKtMwsHfgU8Ft33xh1PXKYmZ0C3ESwEoMkrjyCZZe+DFwL1BOcXbnLzKa7+/cH4k1S8UjjaI65Mq5E7ptAG/ClqAuR9/gJ8CN33xt1IXJMWcAY4CvuXufuHe5+L/Ag8HUzGz0Qb5KKoaGVcZOMmX2S4C/ZD2umf2IxswuAucD/iboWOa794e2aHttfAUYTXEnab6l4emotMCfGdq2Mm4DC68f/kWBpmPqo65H3uBRIB1aYdR2sTwtvHzGzFuDr7q7LcKO3AVjAew8G2sPbATlISMUjjYFYGVeGgJmVEFwOfUl4hRtm9tGw86IkAHe/3d1PdvcFnR/Ar8K7rwi3KTASw0Ph7bwe2+cCh4A3BuJNUjE0/kBwRPFjM8swszSCa5ffQYfYCcPMioFfE+yvS8ysJAyRK4EZUdYmkqTuJphD873Oiczh6cUbgO93zonqr5RcsDA8svhXYBHBJZyvA19095pIC5MuZrabo8/H+I67f3sIy5E4mNkVwA8ITk9NJVh1uiU8+pAEEM6n+THB6gpNQDPwC3f/9YC9RyqGhoiIDI5UPD0lIiKDRKEhIiJxU2iIiEjcFBoiIhI3hYaIiMRNoSEiInFTaIiISNwUGiIiEjeFhoiIxE2hISIicVNoiIhI3BQaIkPEzMrMbJ+ZdZjZk+G2X5rZu2b2jpl9OuoaRY5HCxaKDCEzuxG4B/iMu//GzAoI+iCcp66FkgwUGiJDzMweAC4h6LL2O+AH7v54tFWJxEehITLEzGwasI6gDefD7v7JiEsSiZvGNESGWNja9jvAZOCZiMsR6RUdaYgMsbAFcTkwCsgHznT3hkiLEomTjjREht4XgErgGiAL+PdoyxGJn440RIaQmZ0M3EdwtdQhM/ss8CvgSnd/ONrqRI5PRxoiQ8TMvg+8AEwDPhVu/lx4W2Zm90VSmEgv6EhDRETipiMNERGJm0JDRETiptAQEZG4KTRERCRuCg0REYmbQkNEROKm0BARkbgpNEREJG4KDRERidv/ByERQVS08GEzAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the initial conditions.\n", "pyplot.figure(figsize=(6.0, 4.0))\n", "pyplot.title('Initial conditions')\n", "pyplot.xlabel('x')\n", "pyplot.ylabel('u')\n", "pyplot.grid()\n", "pyplot.plot(x, u0, color='C0', linestyle='-', linewidth=2)\n", "pyplot.xlim(0.0, L)\n", "pyplot.ylim(0.0, 10.0);" ] }, { "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", "We will implement Burgers' equation with *periodic* boundary conditions. If you experiment with the linear and non-linear convection notebooks 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", "$$\n", "\\begin{equation}\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", "\\end{equation}\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." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# Integrate the Burgers' equation in time.\n", "u = u0.copy()\n", "for n in range(nt):\n", " un = u.copy()\n", " # Update all interior points.\n", " u[1:-1] = (un[1:-1] -\n", " un[1:-1] * dt / dx * (un[1:-1] - un[:-2]) +\n", " nu * dt / dx**2 * (un[2:] - 2 * un[1:-1] + un[:-2]))\n", " # Update boundary points.\n", " u[0] = (un[0] -\n", " un[0] * dt / dx * (un[0] - un[-1]) +\n", " nu * dt / dx**2 * (un[1] - 2 * un[0] + un[-1]))\n", " u[-1] = (un[-1] -\n", " un[-1] * dt / dx * (un[-1] - un[-2]) +\n", " nu * dt / dx**2 * (un[0] - 2 * un[-1] + un[-2]))" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# Compute the analytical solution.\n", "u_analytical = numpy.array([u_lamb(nt * dt, xi, nu) for xi in x])" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEbCAYAAAAmmNiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXxU5dn/8c892fedQAgJSYCwhCRsAmELKmpBFGtdWhVta936q7XW1j5oVVwera2PtVVbt6pVqrhUba2VupCKQGQx7BDIwpZ9gez7XL8/JkTCOkkmmUlyvV+veQ1zzplzruEk88055z73bUQEpZRSyh4WZxeglFKq/9DQUEopZTcNDaWUUnbT0FBKKWU3DQ2llFJ209BQSillN6eGhjFmmDHmY2OMtvtVSql+wGmhYYy5DFgPJJxlOQ9jzEPGmD3GmB3GmHXGmNl9U6VSSqnjOfNI41fAAmDtWZb7I3AVMEdEkoC/AJ8YY1J7uT6llFIncGZozBKRfWdawBiTCNwEPCYiZQAi8iKQBzzS+yUqpZQ6ntNCQ0Ra7VjsMsAAq0+Y/jlwgTHG3+GFKaWUOi1Xbz2VDFiBgydMzwfcgfF9XpFSSg1i7s4u4CzCgXoRaTthenX7c9ip3mSMuQnbaS28vb2nxMTE9F6FyiGsVisWi6v/DaNA91V/0ZP9tHfv3nIRiTjVPFcPjdMxZ5opIs8DzwMkJiZKdnZ2nxSlui8jI4P09HRnl6HsoPuqf+jJfjLGHDjdPFf/c6Ec8DXGuJ0wPaD9uaKP61FKqUHN1UNjG7YaR5wwPQ5oBXb3eUVKKTWIuXpovAcIkH7C9PnAf0Skps8rUkqpQcylQ0NEsrFdm/gfY0w4gDHmB9juIr/HmbUppdRg5LQL4caY32K7Izym/fWW9lnniEjzcYv+BLgfWGuMaQFqgAtEZAtKKaX6lNNCQ0R+YedyLcC97Q+llFJO1F+b3CqlgKqqKsrLy2lubj77wg4SFBTE7t3aBsXVnWo/eXp6Eh4eTlBQULfXq6GhVD/V2NhISUkJ0dHR+Pj4YMwZb19ymJqaGgICAs6+oHKqE/eTiNDQ0MDhw4fx8vLC29u7W+t16QvhSqnTKysrIyIiAl9f3z4LDNV/GWPw9fUlPDycsrKybq9HQ0OpfqqxsRF/f+2zU3VNQEAAjY2N3X6/hoZS/VRrayvu7nqGWXWNu7s7ra32dDJ+ahoaSvVjelpKdVVPf2Y0NJRSStlNQ0MppQaA1NRUQkNDGTlyZK9uR0NDKeUSsrOzO774PD092bBhw0nLXHPNNYwaNQp/f39SU1P59NNPnVBp9z377LOMHDmS+vp6h697y5YtXHLJJQ5f74k0NJRSLiExMbHji6+lpYXvfve7VFdXd1pmxYoVvPjii0ydOpUtW7Zw/vnnO6na7gkNDSUmJgY3txNHe+g/NDSUUi5nyZIl5OXlceuttzq7FIe6+uqr+eKLL/Dy8nJ2Kd2moaGUcjmXXnopt99+O3/729945ZVXTrvck08+yahRozDGkJGRAcCnn37K+PHjMcZ0eu/x5/w/+ugj5s+fz9ChQ7nsssuorq5m7dq1XHTRRQwfPpwrrriCqqqqTttqaGjgrrvuIi4ujsTERJKTk3nttdc65m/atInU1FQ8PT254YYbeOKJJ0hLSyMgIIDU1FQeeeSRk2o95sMPP2TatGmMGTOGiRMncsEFF3Ra9+rVq1m8eDGTJ08mJSWF6dOn89FHH3X7/7dHRGRAP8aMGSPK9a1evdrZJfQ7u3btcsp2q6ure3X9119/vbz88svS1NQkU6dOFT8/P8nOzu6Yv3r1apk3b16n10Cnn6H8/HwB5OWXXz5p3YGBgfLggw+KiEhxcbEEBwfLNddcI48//riIiBQVFUlgYKDcc889nd67cOFCSUhIkIKCAhERWbNmjXh5ecmrr77aabnY2FiJjIyU1157TURENm/eLCkpKaet9e233xY3Nzd57733RESkra1N7rzzTgkKCupY5uabb5Zly5aJ1WoVEZG1a9eKj4+PbNy48aTPFxsbKyJn3k9n+9kBNslpvlP1ziClBpCRv/qXs0sAYP9ji3q8Dk9PT1auXMmkSZO4+uqryczMxNPTs8frra2t5Sc/+QkAkZGRzJ49mzfeeIOnn34agKFDhzJnzhxWr17d8Z5PP/2Ujz76iBdeeIGoqCgAZs+ezZIlS7j//vtZunRpp22EhYVx7bXXAjB58mTefffdU9YiIvz85z8nPT2dJUuWAGCxWFi+fDkrV67sWG7ZsmVERER03GORlpZGcnIyL730ElOnTu3x/0lX6OkppZTLio+P58UXXyQrK4tf/vKXDllnWFgYwcHBHa9DQ0NPmhYWFkZxcXHH62OttGbNmtVpXUlJSezfv5/9+/d3mj5u3LhOrxMSEk5ZS3Z2NgcPHmTatGmdpvv7+3P48OGO135+ftx7771MmTKF5ORkUlNT2bFjB3l5eXZ8YsfSIw2lBhBH/IV/Nn3dy+0VV1zBrbfeylNPPcWCBQvw8/Pr0fp8fX07vT7Wkd+J09ra2jpel5eXA3DllVd2avlUX19PZGQkFRUVne6PsLdPsGPrDQ0NPe0yVquVxYsXU1VVxapVq4iOjgYgPT2dpqYmu7bjSBoaSimX9+STT7Ju3TpuuOEG/vCHP3Sad+xL3HYq3qaurs6h2w8PDwfgX//6FzExMQ5f75EjR067TE5ODuvXr+eJJ57oCAxn0tNTSimX5+XlxVtvvUVjYyM/+9nPOs0bMmQI0PmLNzs726HbX7BgAQBbt27tNL2goICrrrqq24NgJSYmEhMTw6ZNmzpNr6ysZMaMGRw9erTjaOLEPqOOP33WlzQ0lFL9wpgxY3juuecoKSnpND0hIYHo6Gjef/99wNY0dsWKFQ7d9nnnncfixYu57777Or6s6+rquOOOO4iMjOz2BXpjDE888QSrV6/mn//8J2DrvXjZsmUkJCQQHBzM2LFjiY+P5+WXX+4IxrffftvhwWi30zWrGigPbXLbP2iT264baE1u9+zZIykpKRISEiIjRoyQ2bNnn3K5G2+8sVOTWxHbz09SUpKMGjVKLr74Yvn0008FkBEjRsjNN98sIiLp6ekSEhIiHh4ekpKSIuXl5bJkyZKzTsvJyRERkcbGRrn77rtl5MiRkpSUJKmpqbJ8+XJpbW0VEZGcnBxJSUkRDw8PCQkJkZSUFMnMzOyo8eGHH5aEhAQBJCEhQZYvX94x7x//+IdMnTpVRo0aJUlJSfKTn/xE6urqOubv2LFD5s+fL5GRkTJv3jy54447ZMqUKeLn5ycpKSnS1NTU8X93rO4vv/zytP/XPWlya+S484ADUWJiojgtkZXdMjIySE9Pd3YZ/cru3btPaqXTF3S41/7hTPvpbD87xpjNInLKtrx6ekoppZTdNDSUUkrZTUNDKaWU3TQ0lFJK2U1DQymllN00NJRSStlNQ0MppZTdNDSUUkrZTUNDKaWU3TQ0lFJK2c3lQ8MYM9UY829jzG5jzHZjzAZjzBXOrksppQYjlw4NY8xI4DOgHJgoIhOBvwBvGWMWO7E0pdQAUVhYSGpqKv7+/v2m/7Mbb7yRmJgYjDEnjRrY21w6NICFQCDwfyLSCiAifwaqge85szClVO8qLCzEzc2NW2+9tVe3ExUVxZYtW7o91vb777/P73//+5OmZ2VlERoayoYNG3pa4klefPFFHnzwQYev1x6uHhqt7c8dIwwa20gkFsDtlO9QSg0Ir732GhaLhZUrVzplWFN7nS40/Pz8iI2N7fHwtK7G1UPjTWAPcK8xxt8YYwGWAV7An51amVKqV7355ps8+uijHDlyhH/84x/OLqfLxowZQ1ZWFhMmTHB2KQ7l0mOEi0i1MeY84GVs1zVqgSpggYj893TvM8bcBNwEEBERQUZGRh9Uq3qitrZW91MXBQUFUVNT0+fbbWtr6/Xtbt68mYSEBJYuXcojjzzCSy+9xEUXXdQx//LLL2fr1q2UlpayevVqHnroIXJychgyZAjPPPNMp7EiPvjgA1588UWqq6tpbm4mJCSEBx54gHPOOeekzwW2cShWrlzJfffdR1FREePGjePee+9l8eLFPP/88zz99NMcPXqU2267jczMTL7++mvq6upITk7uqC05OZn777+f7du386tf/Yply5Z1bGfPnj3cd9997Ny5k6CgIDw9Pbnkkku46aab8Pf3Jz8/n9/+9rds3boVYwwiwtVXX81tt93WMR46QGNjI2D73TnV/jjTfmpsbOz+79vpRmdyhQeQCBwGngd8sR0ZXQ1UAN+yZx06cl//oCP3dd1AG7nveD/+8Y87fiZuv/12cXd3l+Li4k7L3H///QLI7bffLm1tbdLS0iJz586VmTNndlruwgsvlOeee67j9TvvvCN+fn5y8ODBTsvNmzev04iAa9asEUDee++9TsstW7ZMnn766Y7X119/vcTGxp7ycwBy//33d7zOycmR4OBgufPOO8VqtYqIyN///ncxxkhWVpaIiLzxxhsyf/58aWhoEBGRoqIiGT16tDzxxBOd1v3yyy8LIPn5+afc9pn2U09G7nPpIw3gISAY+KmINLRPe9MYcxXwqjEmStovkCul2j0QdPp5F/8epn7f9u9NL8OHd5xhPVXf/Pu5uVC0FYCTxoKbfD1c8gfbvwuzIGpSl0s+XnNzM1lZWTz99NMA3HLLLfzhD39gxYoV3HnnnSct//3vfx+LxYLFYmHx4sX84he/oKmpCS8vLwD++Mc/Eh8f37H85Zdfzo9//GP+9re/cffdd5+2jlmzZhEfH89rr73GkiVLANsf2W+99RaZmZnd+mwPPPAAbW1tPPTQQ9guz8Jll13G7NmzsVhsVwsuvPBCzj33XLy9vQEYOnQo3/72t3nhhRdO+fn7mquHxkTg8HGBccxeYAkQB+zr86qUUr3mww8/5Dvf+U7H63HjxjFnzhxeffXVU35pjhkzpuPfoaGhAJSWljJixAgAvL29ufXWW9m4cSNWqxVjDJWVleTl5Z2xDmMM1157LY899hhHjhwhJCSEjIwMkpKSCAsL69Zn++STT5gwYQK+vr6dpn/xxRcd/w4ICOBPf/oTb775JlVVVbi7u1NcXMyRI0e6tU1Hc/XQKAVSjTHuJxxRxAICuMb/olKu5PgjhDOZ+v1vjjrO5uZvvtTOOEZ4D48yAF5//XWys7N59dVXO6YdOXKEgwcPsmXLFlJTUzstf/wX8LG/1o9dn6irq2P+/PlER0fz+eefExISAsDIkSPtapF13XXX8eCDD7Jy5UpuueUWXn31VZYuXdrtz1ZeXs6UKVPOuMy9997LU089xWeffUZaWhpgO0JZvnx5t7frSK7eeuqP2O7TeLC9qS3GmPnAt4GVIlLuzOKUUo5VVlaG1Wpl586dbNmypeOxY8cOPD09OwWJPdauXUtubi633357R2B0xahRo5gxYwavvfYa9fX1ZGRksGjRoi6v55jw8PCzHjH89a9/ZcGCBR2B4WpcOjRE5B3gImAmsMsYswN4CrgHuMGJpSmlesGKFSu48MILT5oeEBDAnDlzWLFiBS0tLXav79jRxLHrBwBWq5WysjK717F06VLWrVvH448/zqJFi/D09Ow038PD41jDHerq6s7YPHjBggXs3LmThobOZ9yvvPLKjtZMTU1NneoFKC4utrve3ubSoQEgIqtEZL6IjBORJBFJFpHfiojr3u2jlOqWv/71r1x88cWnnHfxxRdTVlbGv//9b7vXl5aWRnBwMM8++2xHE9UnnniC+vp6u9dx1VVX4enpycMPP3zKU1NxcXGUl5fT1NTEunXruOOO0zcueOCBB7BYLDzwwAMdQfP666+TlZXF9OnTAVi0aBGffPIJ27dvB2Dv3r2sXLnS7np73emaVQ2Uhza57R+0yW3XDaQmt+Xl5ZKcnCxubm6SkpIie/bs6TT/qaeekoSEBAFkyJAhsnDhQomMjBRAUlJSZNu2bfLoo4/KiBEjBJBx48bJ22+/LSK2prPTpk2TqKgoSU9Pl+XLl8vw4cMlJCREzjvvPCkoKJCUlBTx8/MTPz8/SUlJkZKSkk7bX7JkiZzuu6SkpETS09Nl9OjRMmHCBPnggw/k448/lpSUFAEkMjJSFi9e3LH8zp07ZdGiRRITEyMpKSly6aWXSm5ubsf8yspKWbp0qURGRsqMGTPkyiuvlKVLl3Z81rVr18oPf/jDTp/1hRdeOKmu3mpya6Q97QaqxMREyc7OdnYZ6iwyMjL6TWdxrmL37t2dbmLrK2e8ED5APfjgg7i5uXHPPfc4uxS7nWk/ne1nxxizWURO2RmXq7eeUkopp3vvvfd4//33nV2GS3D5axpKKeUM8+bNo6mpiTVr1jBs2DBiY2OdXZJL0CMNpZQ6BWMMY8eOJSIigtdee83Z5bgMDQ2llDoF7UDz1PT0lFJKKbtpaCillLKbhoZS/dhAbzKvHK+nPzMaGkr1Ux4eHid1R6HU2TQ0NODh4dHt92toKNVPDRkyhIKCAurr6/WIQ52ViFBfX09BQQFDhgzp9nq09ZRS/VRgYCAAhYWFXerEr6caGxs7BghSrutU+8nDw4PIyMiOn53u0NBQqh8LDAzs0RdAd2RkZDBpUs/HzVC9q7f2k56eUkopZTcNDaWUUnbT0FBKKWU3DQ2llFJ209BQSillNw0NpZRSdtPQUEopZTcNDaWUUnbT0FBKKWU3DQ2llFJ209BQSillNw0NpZRSdtPQUEopZTcNDaWUUnbT0FBKKWU3DQ2llFJ209BQSillt34RGsaYy40xXxhjNhtj8owxm4wx1zm7LqWUGmxcPjSMMT8D7gG+JyJTgERgL3CeUwtTSqlByKXHCDfGjAQeA2aLyGEAEWkxxtwFRDmxNKWUGpRcOjSA64CjIrLx+IkiUggUOqckpZQavFz99FQasL/9msYaY8weY8w6Y8wPnF2YUkoNRkZEnF3DaRljdgAjge3AZUApcDnwBnC/iDxymvfdBNwEEBERMeWtt97qk3pV99XW1uLv7+/sMpQddF/1Dz3ZT/Pnz98sIlNPNc/VQyMHSADmiMiXx01/F7gIiBCR+jOtIzExUbKzs3u3UNVjGRkZpKenO7sMZQfdV/1DT/aTMea0oeHqp6dq2p+3nDA9C/AFxvdtOUqpQaWlwfY4ZvMr8Oc58PkjUD04L6u6emjsaX8+sc6200xXSqnua66H3NXw6XJ46UJ4dATs+ddxCxgo3gZfPA5PJsHK66A8x2nlOoOrt576J3A1kAx8edz0JKAB2OmMopRSA4gIrPmdLSwOb4S25uNmGqjM++bluMXgHwlb34A9H8Luf0BFDty2vs/LdhZXD42VwB3Aw8aYi0Wk1hgzB/gO8KCI1Dm3PKVUv2K1Qsl2OLAept8Mxtge2f+Ggs2AgWEpEDcXYmdDzAzwCf7m/b6hkHiR7VFdBH+aCaW7oGwvRIxx2sfqSy4dGiLSZoy5CPgNsNMY0wg0Af9PRF5wbnVKKZcnAhW5kP/f9scaaKi0zYubC5Htl0Xn/gKsrRA7yxYM9ggcBqMvhG1v2o46Iu7snc/gYlw6NABEpBL4kbPrUEr1M0f2w8sLobqg8/TAaIhPB8txX3+J3+reNs75EYy5AEYt6GaR/Y/Lh4ZSSp1RfSXs/9J2JNHaBJc+bZseGA1NNeAbBiPnQPw8iJsHofG2U1KOED3V9hhENDSUUv1Lcx0cXA/5X2DN/S+meCsG2/1mVosnu1N/TVR4CCF+nnDLlxA0Aiza0NJRNDSUUq6trQVpquVQgxdbDx/FbctfWZj/KGBrc98k7myRUXzZlsRaaxJb/rQBi8WN+WOHcPnkaM4NAM/ezIyjhyDjMVurq8sH/qVWDQ2llGuxWmkr3k7p1v/QmrOaiMrNvG/O5Vf11wIQbSIZ7hHPOusEviKJwsBUQoODsBgDLW2MaW5jX2ktn+wq4ZNdJYT7e/LsNVM4J87OC9xd5eEDW/8Gxg0WPQHegb2zHRehoaGUcjqrVTi87k1at71LRPkGAqxVDDtufmhrCSG+HkyKCSFp+GhKhy/kkqhAbgr0xs1y8vWJ0ppGPsgq5O3Nh9hbUssPXtnI3340neTo4JOW7TG/cIiZCQfWQs4nkHS547fhQjQ0lFJ9TqoLKd26ii9aJ/DZYTcy8yv4SfMH/ND9EwAKJIwt7ilURc7Eb9y5TEgcy9cRfhg7L2APCfDmR3Pj+cHsOH76ZhYfbiti6V82sPKmmSQODXD8Bxq7yBYae/6loaGUUj1WX0n5js84svMTAorWM7T5IJHA5pYb+bjtXADWBZ5LYNh4/Medx8SkVBaF+vV4s24Ww5NXpdLQ3MZne0q59qWvePvmmYwM7/m6O0lcCKuWwd7/2FpwuXs5dv0uRENDKdUrymqaWJdTRvIn3yW2fjvhCOHt82rFmy2W8cTFJfBY8kRmJoQRE+pr95FEV3i4WXjmmsl8/+WNrM+r4JfvbOOtW2Y6diOhcRA2ytalSPleGDrRset3IRoaSqmea2uhNi+T4i2r4NAGfsLd7C5tBOB1jzaiLG5sM2MoDDkHj9HzGZ06l1nDQpjdCyFxKt4ebjy3dAqzHv2cDfsr2X64ionRQY7dSHCMLTSqizQ0lFKqE6uVxoJtFGV9jDUvg6ijWfjTyKj22T5N2/DxGMfUkSEcGP4bgseOZnJMJNNOcdG6rwR6e3DltBG89GU+L6/N5/+uSnXsBiZeabsgHhrn2PW6GA0NpdTZidBSf5Rt5cIHOc18uHMlvyu7heO/HnMkihy/KTTHzGHZpItIThiBp7tr3VR3Q9pIXl6bzz+3FfKrhWMZEuDtuJWnftdx63JhDg0NY8yzInKbI9eplHIOa1UhRVmrqMv+nLDS9RS1BnB508PtcwO4zjOeUu94GqJnET7xfFLGj2eUl2v/HToi1JcF4yNZtbOE1zMPcueCwdEzrSN1aQ8bY5aeZZGFPahFKeVEIkLR7kxqMl8hqHgdQ5sPMvy4+VaaSAz3YLiPcMXcZGLivyLFz9Np9XbX92fFsWpnCX/76gC3pSfg7eHmmBU3HIWDmWBxg9EDtwPDrv5Z8MoZ5rnuYONKqZM113NkzxdkVVj4qDyS9bkVTKr+nKc93wSgTrzY6jaB8ogZ+Caey4TJaawK9rONPT1x2FlW7rqmx4Uyflggu4qq+efWQq6YOsIxK67MhTeusl0E19DosJuTjyb8gLHA94BnHFGUUqoXtLVQk2tr4eRxcA3Da3cQQivlrfN4p/VmADx9Uvkw5AZbC6dJ85gZGdwrzWCdyRjDD2bHcdfbW/nL2v18Z0q0Yz5jQJTtubqo5+tyYV0NjdtF5MAppu8yxvwbeBNY3fOylFKOUN/cysb9R/Bc/QApRe8SQAPH7oe2imEH8VjDE7ln8jjSRoUxbmggFst3nFpzX1icMoxH/rWL3UXV5JfXER/h3/OV+kWAsUB9ObQ2g3v/O3Vnjy6Fhoh8doZ5DcaYsT0vSSnVLSI0l+VS8PXHtOSs5s/yHf5ZHERLm3Cney0z3RvIlShy/CbTEjOX4ZMWMCFhJEku1sKpL3i5uzEzIYyPthfzVX6lY0LDzd02fnhNEdQW2+7bGIC6eiF87qkmAyHAEqDREUUppezTVlVIYdYq6vd8RnhZJmFtZR3NYENaImm1LiI5OghG/IjM2LtIGT+eBE8HXfjt56bHtYdGXgXfPcdBX/ABw2yhUaOhcUwGp77gbYDDwLU9LUgpdXrSVEtuFazLLWf9vjIey7uUEdR1zK8Uf3Z4JFM1LI25yYu4fcJEgnw9nFix65oeb+sq/av8SkTEQdc12hsIVBf2fF0uqquhkQvceMK0NqAUyBWRNodUpZSyaa6nfPd/Kd/2H3wLvmRI434ubnyORmwd4i3wmMQwj3oqhszAL/FcJkyexdwgXycX3T+MGRJAsK8HRVWNHKpsICbMAf9vge2hUVfW83W5qK6GxjMi8t9eqUQpBUBFyUHKMp7D8+CXRNftIJzWjo7+msWN6X7FBI2aSVpCGNNGvc2IUA2J7rBYDOeMDOU/u0rIzK9wTGic+2tY8CB4OrgXXRfS1Qvhv++tQpQalKxWag5uITvvAB/WjmF9bgU1Jfms837aNlsMO008h4KmYRLmETfpfF6JHjLgmsE6y/T4MP6zq4Sv8iq50hH3a/j0wiBPLsa17/lXaqARoaFkH4c2f0xbbgZRRzYSJNUEW6N4pfl3AHh7DOG9wGvwGTGJ6EkXMC4+hglO7OhvIJsed+y6RoWTK+k/NDSU6mVNrW1kHTxKeeYbTM95ighrKcf3eFQkoRz2Hc/PZsUxY3QkqTHBeLl/y2n1DibjhgUS4O3O4SMNFBxtYHiwT89WWF0If78J3Dzhur87pkgXo6GhlIO11FZy8OtV1O35nIzGMTxTmkRTq5V0y1Eu9iylUvzZ7ZVCTVQaYUkXMGHiJNK9PEh3duGDkJvFMG1kKJ/vKeWrvAq+PTm6Zyv08IH9a8DTAfd9uCgNDaV6qK2pjoNZn1O161MCi9cR27SPBGNrmX6o7RyaWsczdmgAo0d+i/URsxg/KY1ZvgN3OND+ZnrcsdCo7HloeAeDuw8010JjNXgHOqZIF6KhoVQXWVua2VNSx/r9R1mfW8HS/LuYS1bH/Gbc2OE2lrLwGQSNu4DNU+cT5q8h4aqmx4cBDrquYYyt2W1lnu0GPw0NpQYfsbZxYNcmSretwvvwlyTUb+W+pl+ySWy95sS5jSPSs4qisOl4jp5PwpTzSQkPc3LVyl5JUYH4ebqxv6KekupGIgN7ODBTwLHQKISIgTdeh4aGUicQEXKKj1Kx5iU8Dq4hvvZrRlLNyOOWmeV7mJgx5zEzIYyZ8elEh/qR6KR6Vc+4u1mYMjKUL/aW8VV+JZekRPVshcfuCq8p7nlxLkhDQw16IsL+/Xns3/4l79ZOJDOvkvLaRjK9nmaoOQJACWHsD5yKxM0jevKF3BGToPdKDCBTYkL4Ym8ZOwqqHBAaQ23PA7QrkX4XGsaYNcBsIE5E9ju5HNUPiQj5hwo4mPUJJv+/jDi6kXgOEyuGnzX9maMEEBHgzZqQ64gND2TYpAuITmvOwikAABe3SURBVJhIpGXw9QY7WIwbZuswfldhdc9XFpsGzXUwNLnn63JB/So0jDGXYwsMpewmIuSW1bI+r5IDuzdx2YFHGCt5xJtv+t6sx5t8vxSWz40iKXky8eF+GHO+E6tWfWncMNsF691F1T3vvHDsIttjgOo3oWGM8QQeBT5CxyJXZyAi7Cs6Qu6WL2jNyWD/0WaeqLf9EgfTwjKvPNqMG/k+46kbPouwiQsYPmE2E9y9mODk2pVzRIf4EODtTkVdM2U1TQzp6cXwAazfhAbwY2ATsBcNDXUcq1XYW1LF3q2ZNO/LYGjFV6TKLsYY2/AuRRLKX/0vY0ZCODPiQyn0/jvDx57DKK+BewOW6hpjDOOGBbIhv5KdRdU9Cw1rG5TsgPpKSJjvuCJdRL8IDWNMKPALIA24wbnVKGezWoXdRVVsyC0j80AVG/Ir+W7T2/zS461vFjJQ4jmCqqFpBI4/nw3T5mPcjv24xzqlbuXaxreHxu6iauYnDun+iqyt8NxcsLjDvWUwwK6F9YvQAO4DXheR/dpiZfBpswq7i6rZunsPDdmfE172FVNlO0WtC1jVthiAvIBkjsrnHIlMw3/ceYRPXEBk0HAinVy76j/Gt1/X6PHFcHcv8AmFhkrbeOH+PQggF+TyoWGMGQVcCYzrwntuAm4CiIiIICMjo3eKUw5TW1vbsZ+sIhystrK70opP6WZG123mHHZyjaXgmzcYON83l9ZYT8aGuhHhPYEt5kXbHbmtQNY+YJ8zPsqAd/y+Gkjqq2xjyG3OLe7x55tqCcSfSjat/ie1AQkOqK7rems/uXxoAI8Dj4lIlb1vEJHngecBEhMTJT09vZdKU47Q2mblrx+upkTCObLnS94qHkZpk+1H81mPz1jotgGAJosPlWFT8Ek8l+AJ53NO5ETOGWCH/v1BRkYGA/F3qrGljYe+WkVJvTA9bQ4+PRlL/fAoyNnP1DHDITHdYTV2RW/tJ5cODWPMHCAJuMrZtSjHaW2zsrOwmg05xZTuWU9g8TqmyXYmmX14mVZ2tPwS79A0ZsSHEuxzHVWWuQSNPx+v4VMY5u7p7PLVAOXt4UZChB97S2rJLqkhdUQPBlTquCt84N3g59KhASwA3ICNx13LaL/dko+MMc3AMhH5yBnFKfscC4nMvAoy8yrYsr+U31p/x3ctu/E3jWAAA1YMR4LG88S8JMImH2t1kuLM0tUgM25YIHtLatldVO2Y0BiAd4W7dGiIyH3YLoJ3MMY8ANwPLNQ7wl1Tm1XYVVjN+twycrO341Owjti2Azzaen37EoYxPqX4SyM1/nFYEuaR3xpJ0qJbCPENdWrtanAbPyyQD7YU9vxieFB7F+tVBWderh9y6dBQ/YPVKuwqsh1J7N67F49DXzK5dRsL3XYSbcptRxLuUDzhR4wdO44Z8WEMq/krBEYREGjr56c8IwM0MJSTHX9neM9WtBhiZ0HQcAdU5Vr6TWgYYxYC/8sJp6dEJNWJZQ1KVquwt7SG9bkVZOaUkbn/KFUNLZxjdvOW10O2hdp/spo8grCOnIPPmPn8T9JU8AmxzQia6pzilTqDY6Gxp7gGq1WwdHdsdt/QAftHUL8JjfbrFnrtwglEhLzyOtblVrB5XwFN+etIbt5CmmUnQyScVS13MDzYh4SRs2nODcQaNQXvMedC/Dy8IicOuJub1MAVEeBFRIAXZTVNHDpST2yYn7NLcjn9JjRU3zpUWc+63HLW51ZQlLOFaQ1rmWXZyZWWvXiZ1o6fnDFeVay5KZ0Rx365rPvB0oOmiko52bhhgZTVlLG7qLpnofHRL6FkJ1zx8oC6wU9DQwFQXNXI+rxy1u8royQ3i31VbhQSDsCNbhu5y+NtAARDU8REPEfPx8TPwydmJiM8j/vF0sBQ/dz4YYF8sbeMXYXVXJQ0rPsrOpQJRVvh6EENDdX/VdQ2kZlXybqcMg7k7CS6ahOzLDv4pWUX4aaa57wv5+uE25gZH0Z66FAk1wcTPxczcg5eA/RcrVJw3NgaPb0YHjTCFhpVhyB64FzD09AYJGoaW9iQX8nanArW5Zazp7iGu9xXcqvbWlsLJ49vlm3xG8qPzhmPZd6xH/Q4GDfFKXUr1dcmRAUBsL3A7k4oTq2j2e3hHlbkWjQ0BqjGljY27T/Cutxytu47QEBxJtPNTt5vvYxKAvFytzAxsJ7o+nJavYKxxM3FkjAP4tLxCEuw9eGk1CAUH+5HgJc7JdVNlFQ3EtndbtKDRtieNTSUK2pps7Lt8FHW5lSwMacAt0NfMZ3tXGjZyc9NPm4etlHqIsbPJWz6BUyKCcb7SBy0/Rp3beGkVAeLxZA0PIj1eRVsO1zFgvHdDQ090lAu5NgNdetzK1iXU8aG/Ueoa27DmyayvG7Gx735m2WNO23Dp+KWkM7FEy+E8DDbjCFjnVS9Uq4tecSx0DjKgvHd7GC/40jjkOMKcwEaGv2EiJBfXsfa3Aoyc0opy93CxOYtzLTs5FxTzHnNvyMhwp+0hFiaDozFwxPc2083WWJmgI5Sp5TdUqJt/U5tPdyD6xohsTD6Ahgy3kFVuQYNDRdWeLSBdbkVrMspZ1dOHqn1a5ll2cHy9hZOx1+83nRrHOGx7SNct34G2husUt2WHG27GL7t8FFEhG4N/uYXDte87eDKnE9Dw4VU1jWzPreCtbnl7NmXQ/WRUnLEdl50sjnAY14vdizb6jcMt1HzMHHpEDeX8OP7uNHAUKpHhgf7EObnSUVdM4cqG4gJ83V2SS5DQ8OJapta2ZhfydqccrbsO0BI2QbSLDu5wbKDMZYC1ntO5MW4J0kbFU7ayBlI5g5M7EyIS8ddWzgp1WuMMSRHB7E6u4yth492PzQajthu7guIAv8IxxbpJBoafaiptY2vDxxlXW4563Ir2HroKN9iLT90/4j/Mfm4eUrHsm3uPkxPiGHm1VO/CYfvvOSkypUafJKjg1mdXca2w0dZnBLVvZWsuge2rIDFf4Ap1599+X5AQ6MXtVmFHQVVrM0tJ3NfKU0HNjJNtrPOOoGvZQwWA+PCDKm1eVgtHliHT8ESnw7x83AbPlVPMynlRCkjbNc1enQxvKPZ7cBpQaWh4UAiQk5pLWtzylmXU0Z5/lZSW7aSZtnBdZY9BLg3AJA5xEJt+vc4Jz6UwNZpUJSOJXYmeGqPmkq5iuT2FlQ7Cqposwpu3ekmfQDeq6Gh0UOHKus7Ll6vy62grKYJgJWeDzLdsqdTC6fWkATcE9KZMf4SiD/W9nsIjD6/7wtXSp1RuL8Xw4N9KDjaQG5ZLWMiA7q+kgF4V7iGRheV1zZ1NIPdk5NDbNVG0iy7+IXbDq5uvpeIgBGkJYQR2jiZtrKjuCWkQ/w8iJuH+wAcxUupgSw5OoiCow1sPXS0h6Ghp6cGjZrGFr7Kq2RdbgUbcooYVvYlaZad/KC9hRPHXXZ496I2hsw9z9amu+X/wN1bWzgp1Y8lRwfz7x3FbDtcxRVTR3R9Bcf+UKwqAKt1QHTXo6FxgsaWNr4+cIS1ueVsyimktjCbnW0xAHjRzLtef8TLtAC2Fk6W2FmY+LkQN5fIoSnfhISHj7M+glLKQVKOu8mvWzx8wDcc6suhrhQChp79PS5u0IdGa5uVbQVVrM+tYP2+EpoPbmKa7GCWZQe3W/bR7O7OD4e/zYxRQ5iZEI5b9g/BJ1hbOCk1CCS1h8buohqaWtvwcu/GIGPfe8v2neEb7uDqnGPQhYaIsLekvYVTbjlf5VUS05zDz9zfZulxLZzANkqde+Q43vpe/DeHmQm/cVLlSqm+FujtwdihAewprmF9bgXpid0YgS96YI1FMyhC42CFbbzrtTnlHMzdyfjGLGrFh0+taQAMCfHl/IYsANpCEjouXpuRc3DTUeqUGtQuShrKnuIaPtxW1L3QGGAGfGgU17TwuyceZpZlJ3e77egYpW6/93jmnn8LaaPCGR7kDdt8YeRs3LSFk1LqOBcnR/H7T/examcxj1yW1PVTVAczYdNfIGoyzLild4rsQwM+NEZzgD94PtPxus0rGEv8HEYmnMvI41tDpFzlhOqUUq5u1BB/xg0LZHdRNV/sLe/6+Bq1JbBtJTTVuHRoNDQ0sOlwLetyK8643IAPDTEWJOE8TPw8iJuL29BksHTjYpZSatBanDKM3UXVfLitsOuhceyu8Mo8xxfWAy21lRz8ehW1uz8jvDSTHS3DuLn5jrO+b8CHRq1/HOa6vzu7DKVUP3bxxCge/zibT3aV0NDcho9nF/7wjBgHXkFQtgeKtsKwlN4r9AysViF/+5fUbn6HwOJ1xDTtI8F800kqpp6JUYGkjQpn2Rna+wz40AC9uU4p1TMxYb6kRAex9XAVq7NLWThxmP1v9vSFSddA5rPw1fOw5Jmzv8cR2loo3J7BV2UefFLkzfrcCi5vep97PVYA0IwbO9zGUhExA9+x5zF2Sjr/DLD1f7fsDKsdBKGhlFI9tzgliq2Hq/hwW2HXQgNg2o2Q+SfY/jYsWG4b1c/RrFYq8rMo+PrfuB9Yw7TaLfjRyKHWJXzUeiUAOwKmszrQisfodEZPXUBKeFiXN6OhoZRSdlg4cRgP/2s3n+0upbapFX+vLnx9hiXAmAth78ew+RWYe5dDajpa30xmXgWBGb9mfMUqwqSa42Mgj+FEDR3OI9OSSEsIZ2SYL8b0bFwPDQ2llLJDVLAP00aGsHH/Ed7PKuDaGbFdW8Gsn0JsGky6rts1NFQcYv/Gj2nOyeA3bd9lfbFBBB53LyXNvZoiCSPXfzKtsXMpswzl8m9fQXx3unQ/Aw0NpZSy0zXTY9m4/wj/+9FupseFMrorPd/GptkeXdBcU8H+TR9Tn72a8LJMotsOMa59XkhzPB6WNCbHBtM47HZ2xN5H4vhJDGu/jyQjIwOLgwMD+kFoGGNSgR8Dk7HV6wF8CjwkImXOrE0pNbhcmhrFf/eW8V5WAbeu+JoPfjwLv66cpjqmtgy8/E/q2LStsYYdZW2sy61g875D/Onwtxlj2jrm14kXuz0nUj0sjesnfZvfTUjuWksuB3D50ADeBHYCc0WkzhgzHPgMuMgYkyIiDWd+u1JKOYYxhkcuS2JHQRX7SmtZ9t52fn9Vqm04BHutfxb++xtorEKiUjkSkEhd2UECqrJpbWvl0sZnOdbqc6tnAh7unlQMmYn/2HQSp8xnqr9zR/jsD6EBcLeI1AGISIEx5rfAi8BC4F2nVqaUGlR8Pd3507WTueTptXywpZDoEB9unpdAoLfHWd/b0NzGoZZwQt2GEEo1lsIsQsniWA93teLN1OB6Ro8Zy8yEcGLi/suQQN/e/UBd1B9CI1lEmk+YVtj+HNLXxSil1KghATz67Yn89M0tPLM6l5fX7ufbk4dzcXIU/l7uuLvZLlAfqqwnv7yO/PI6th2uIrukhjarH/AAfjSQbMljhl8RgUMTGDZ6ChMnJvNOiHOPJM7G5UPjFIEBMAYQ4Is+LkcppQC4NHU4gd4evLAmj3W5FbyeeZDXMw+e8T0WA+OGBZI6IpjpcaGcE7eIqOD+NWCbEZGzL+VCjDFuwNfABhH50WmWuQm4CSAiImLKW2+91YcVqu6ora3F39/f2WUoO+i+OllBrZXPD7aQX2WlTaDNKliBcG8LkX6GSF8LIwIsjAy04OXeN71U9GQ/zZ8/f7OITD3VvP4YGg8Ai4F5IlJ7tuUTExMlOzu71+tSPZORkUF6erqzy1B20H3VP/RkPxljThsaLn966njGmO8DVwLp9gSGUkopx7I4uwB7GWOuA34OnCsipc6uRymlBqN+ERrGmGuBu4HzRaS4fdrF7dculFJK9RGXPz1ljLkGeAH4NXD+cTfRzAGKnFWXUkoNRi4fGsAfAW/gt6eYt7yPa1FKqUHN5UNDRELPvpRSSqm+0C+uaSillHINGhpKKaXspqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaGUUspuGhpKKaXspqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaGUUspuGhpKKaXspqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaGUUspuGhpKKaXspqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaGUUspuGhpKKaXspqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaGUUspuLh8axpghxpgVxpjs9sc7xphoZ9ellFKDkUuHhjHGE/gE8AQmAOOBOmC1McbfmbUppdRg5NKhAVwPJAN3i0iriLQBdwPxwK1OrUwppQYhVw+Ny4GDIpJ3bIKIFAO72ucppZTqQ64eGslA/imm5wMT+7gWpZQa9NydXcBZhAObTzG9GvA1xviISMOJM40xNwE3tb9sMsbs6MUalWOEA+XOLkLZRfdV/9CT/RR7uhmuHhqnY840U0SeB54HMMZsEpGpfVKV6jbdT/2H7qv+obf2k6ufnioHAk4xPQCoP9VRhlJKqd7j6qGxDRh5iulxwPa+LUUppZSrh8bfgVhjzMhjE4wxkcA44F071/G848tSvUD3U/+h+6p/6JX9ZESkN9brEO03920CdgPXAFbgJWA2MElEap1YnlJKDToufaQhIs3AAqAN270Zu4FA4FwNDKWU6nsufaShlHI9xpg12I7240Rkv5PLUX3MpY80uks7OXR9xphUY8wLxpjNxpitxphdxpg/GGMinF2bOj1jzOXYAkO5KGPM5caYL9p/t/KMMZuMMdc5av0DLjS0k8N+400gFJgrIinYTkNeAKw1xvg4tTJ1Su2/W48CHzm7FnVqxpifAfcA3xORKUAisBc4z1HbGHChgXZy2J/cLSJ1ACJSAPwWGA0sdGpV6nR+jK1hykZnF6JO1t7K9DHgZhE5DCAiLcBdwNOO2s5ADA3t5LB/SBaRnBOmFbY/h/R1MerMjDGhwC+AZc6uRZ3WdcBREekU6iJSKCKbHLWRgRga2slhP9DeMu5EYwABvujjctTZ3Qe8rhe+XVoasL/9msYaY8weY8w6Y8wPHLmR/tr31Jl0q5ND5VzGGDfgB8BLIrLX2fWobxhjRgFXYrupVrmuEdh60LgLuAwoxXZ25Q1jzDARecQRGxmIRxqnc8ZODpXT/RpoBX7m7ELUSR4HHhORKmcXos7IG/ADfiEixSJiFZG3gQ+AZcYYX0dsZCCGhnZy2M8YY76P7S/Zb+lNm67FGDMHSAL+5Oxa1FnVtD9vOWF6FuCLrSVpjw3E01PbgLGnmK6dHLqg9vbjP8d2l3+ps+tRJ1kAuAEbjek4WB/a/vyRMaYZWCYi2gzX+fYAqZx8MNDW/uyQg4SBeKThiE4OVR8wxlyLrTn0+e0t3DDGXNw+iJZyASJyn4gkiEjqsQfw5/bZC9unaWC4hn+2PyefMD0JaAB2OmIjAzE0XsF2RPEbY4y7McaCre1yPnqI7TKMMdcAL2DbX+cbY65tD5HFQJQza1Oqn1qJ7R6ah4/dyNx+evE7wCPH7onqqQHZ91T7kcWTwFRsTTh3AHeIyCGnFqY6GGMqOf39GMtF5IE+LEfZwRizEPhfbKenIrF1INrcfvShXED7/TS/wda7QiPQBPxRRF5w2DYGYmgopZTqHQPx9JRSSqleoqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaGUUspuGhpKKaXspqGhlFLKbhoaSiml7KahoZRSym4aGkoppeymoaFUHzHGrDDGVBtjrMaYT9unPWOMOWKMyTfG3OjsGpU6G+2wUKk+ZIy5AngL+JGIvGiMicU2DkKajlqo+gMNDaX6mDHm78D52EZZ+wvwvyLyH+dWpZR9NDSU6mPGmKHALmzDcH4oIt93cklK2U2vaSjVx9qHtl0OhAOrnVyOUl2iRxpK9bH2IYgzAB8gBpggIuVOLUopO+mRhlJ976fAV8ASwBt4yrnlKGU/PdJQqg8ZYxKAd7C1lmowxtwM/BlYLCIfOrc6pc5OjzSU6iPGmEeAL4GhwA/aJ9/W/rzCGPOOUwpTqgv0SEMppZTd9EhDKaWU3TQ0lFJK2U1DQymllN00NJRSStlNQ0MppZTdNDSUUkrZTUNDKaWU3TQ0lFJK2U1DQymllN00NJRSStnt/wO7gu38eOns4AAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the numerical solution along with the analytical solution.\n", "pyplot.figure(figsize=(6.0, 4.0))\n", "pyplot.xlabel('x')\n", "pyplot.ylabel('u')\n", "pyplot.grid()\n", "pyplot.plot(x, u, label='Numerical',\n", " color='C0', linestyle='-', linewidth=2)\n", "pyplot.plot(x, u_analytical, label='Analytical',\n", " color='C1', linestyle='--', linewidth=2)\n", "pyplot.legend()\n", "pyplot.xlim(0.0, L)\n", "pyplot.ylim(0.0, 10.0);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's now create an animation with the `animation` module of Matplotlib to observe how the numerical solution changes over time compared to the analytical solution.\n", "We start by importing the module from Matplotlib as well as the special `HTML` display method." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": [ "from matplotlib import animation\n", "from IPython.display import HTML" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We create a function `burgers` to computes the numerical solution of the 1D Burgers' equation over time.\n", "(The function returns the history of the solution: a list with `nt` elements, each one being the solution in the domain at a time step.)" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "def burgers(u0, dx, dt, nu, nt=20):\n", " \"\"\"\n", " Computes the numerical solution of the 1D Burgers' equation\n", " over the time steps.\n", " \n", " Parameters\n", " ----------\n", " u0 : numpy.ndarray\n", " The initial conditions as a 1D array of floats.\n", " dx : float\n", " The grid spacing.\n", " dt : float\n", " The time-step size.\n", " nu : float\n", " The viscosity.\n", " nt : integer, optional\n", " The number of time steps to compute;\n", " default: 20.\n", " \n", " Returns\n", " -------\n", " u_hist : list of numpy.ndarray objects\n", " The history of the numerical solution.\n", " \"\"\"\n", " u_hist = [u0.copy()]\n", " u = u0.copy()\n", " for n in range(nt):\n", " un = u.copy()\n", " # Update all interior points.\n", " u[1:-1] = (un[1:-1] -\n", " un[1:-1] * dt / dx * (un[1:-1] - un[:-2]) +\n", " nu * dt / dx**2 * (un[2:] - 2 * un[1:-1] + un[:-2]))\n", " # Update boundary points.\n", " u[0] = (un[0] -\n", " un[0] * dt / dx * (un[0] - un[-1]) +\n", " nu * dt / dx**2 * (un[1] - 2 * un[0] + un[-1]))\n", " u[-1] = (un[-1] -\n", " un[-1] * dt / dx * (un[-1] - un[-2]) +\n", " nu * dt / dx**2 * (un[0] - 2 * un[-1] + un[-2]))\n", " u_hist.append(u.copy())\n", " return u_hist" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# Compute the history of the numerical solution.\n", "u_hist = burgers(u0, dx, dt, nu, nt=nt)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "# Compute the history of the analytical solution.\n", "u_analytical = [numpy.array([u_lamb(n * dt, xi, nu) for xi in x])\n", " for n in range(nt)]" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZsAAAELCAYAAAAP/iu7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd3hUZfr/8feTXkkjhJoAAUIJJKEoIiUIWGhSLLgqlnUtP1fXdYv7VVcBy6J++bquukXdRUUURKXooq4CWRFEQUBKKEsVAgECIb1N5v79MUOWQICUmZyZ4X5dV66QM2fOuR+SzCdn7nOeY0QEpZRSyp38rC5AKaWU79OwUUop5XYaNkoppdxOw0YppZTbadgopZRyOw0bpZRSbmdp2Bhj2hhjPjPG6PnXSinlwywLG2PMROAbIPkC6wUaY54yxmw3xmwxxqw2xgxuniqVUkq5gpVHNr8DRgGrLrDey8CNwBARSQX+AXxhjEl3c31KKaVcxMqwuVxE/nO+FYwxKcDdwEwROQYgIm8Ae4Bn3F+iUkopV7AsbETEVo/VJgIGWHHG8uXAlcaYCJcXppRSyuU8/Wy0PoAd+PGM5XuBAKBns1eklFKqwQKsLuACWgKlIlJ9xvJC5+e4up5kjLkbx9tvhISE9EtMTHRfhRay2+34+Xn63wuN58vj8+WxgW+Pz5fHBk0f386dO/NEJP7M5Z4eNudizvegiLwGvAaQkpIiO3bsaJaimltWVhaZmZlWl+E2vjw+Xx4b+Pb4fHls0PTxGWP217Xc0+M5DwgzxvifsTzS+fl4M9ejlFKqETw9bDbhqLHDGcs7ATZgW7NXpJRSqsE8PWwWAgJknrF8OPAvESlq9oqUUko1mEeHjYjswNF7+R9jTEsAY8ydOGYdeMzK2pRSStWfZScIGGNewDGDQKLz643Ohy4RkcrTVn0AeBJYZYypAoqAK0VkI0oppbyCZWEjIr+p53pVwOPOD6WUUl7IW099Vko1QWFhIUePHqWqqsot24+KimLbNt88f8eXxwZ1jy8gIICQkBDi4+MJCQlp1HY1bJS6yBQWFnLkyBHatWtHaGgoxpz3srVGKSoqIjIy8sIreiFfHhucPT4RwWazUVxczI8//khCQgJRUVEN3q6GjVIXmaNHj9KuXTvCwsKsLkV5AWMMgYGBxMTEEBwcTG5ubqPCxqPPRlNKuV5VVRWhoaFWl6G8UGhoKBUVFY16roaNUhchd7x1pnxfU35uNGyUUkq5nYaNUkopt9OwUUqpi1h6ejqxsbF07NjRrfvRsFFKebUdO3bUvGAGBQXx3XffnbXOzTffTJcuXYiIiCA9PZ0vv/zSgkob789//jMdO3aktLTU5dveuHEj48ePd/l2z6Rho5TyaikpKTUvmFVVVdx0000UFhbWWmfu3Lm88cYb9O/fn40bNzJy5EiLqm2c2NhYEhMT8fc/824r3kPDRinlMyZMmMCePXu47777rC7FpaZMmcJXX31FcHCw1aU0moaNUspnXHvttTz44IO8++67vPnmm+dc78UXX6RLly4YY8jKygLgyy+/pGfPnhhjaj339J7G0qVLGTNmDK1bt2bixIkUFhayatUqrr76atq1a8f1119PQUFBrX2VlZXx61//mk6dOpGSkkKfPn2YM2dOzePr1q0jPT2doKAgbr/9dmbNmsWgQYOIjIwkPT2dZ5555qxaT/nkk08YMGAA3bp1o3fv3lx55ZW1tr1ixQrGjRtH3759SUtL49JLL2Xp0qWN/v9tEhHx6Y9u3bqJr1qxYoXVJbiVL4/PyrFlZ2e7fR+FhYVu38eZbrvtNpk9e7ZUVFRI//79JTw8XHbs2FHz+IoVK2TYsGG1vgZqfS/27t0rgMyePfusbbdo0UJmzJghhYWFkpubK9HR0XLzzTfL888/LyIihw8flhYtWshjjz1W67mjR4+W5ORkycnJERGRlStXSnBwsLz11lu11ktKSpKEhASZM2eOiIh8//33kpaWds5aFyxYIP7+/rJw4UIREamurpaHH35YoqKiata555575NFHHxW73S4iIqtWrZLQ0FBZu3btWeNLSkoSkQt/7y708wOskzpei3W6GqUUAB1/90+rSwBg38wxTXp+UFAQ8+fPJyMjgylTprBmzRqCgoKaXFdxcTEPPPAAAAkJCQwePJj33nuPV155BYDWrVszZMgQVqxYUfOcL7/8kqVLl/L666/Ttm1bAAYPHsyECRN48sknmTp1aq19xMXFccsttwDQt29fPvzwwzprERF+9atfkZmZyYQJEwDw8/Nj+vTpzJ8/v2a9Rx99lPj4+JqLMQcNGkSfPn34+9//Tv/+/Zv8f9IQ+jaaUsrndO7cmTfeeIMNGzbw29/+1iXbjIuLIzo6uubr2NjYs5bFxcWRm5tb8/Wps94uv/zyWttKTU1l37597Nu3r9byHj161Po6OTm5zlp27NjBjz/+yIABA2otj4iI4ODBgzVfh4eH8/jjj9OvXz/69OlDeno6W7ZsYc+ePfUYsWvpkY1SCmj6EcXpPGFm5Ouvv5777ruPl156iVGjRhEeHt6k7Z05cakxps5l1dXVNV/n5eUBcMMNN9Q6k6y0tJSEhASOHz9e6/qWiIiIetVyaruxsbHnXMdutzNu3DgKCgr4/PPPad++PQCZmZmNnt+sKTRslFI+68UXX2T16tXcfvvt/OlPf6r12KkXf0ebwaGkpMSl+2/ZsiUA//znP0lMTHT5dvPz88+5zq5du/jmm2+YNWtWTdBYSd9GU0r5rODgYN5//33Ky8v55S9/WeuxVq1aAbVfsHfs2OHS/Y8aNQqAH374odbynJwcbrzxRiorKxu13ZSUFBITE1m3bl2t5SdOnGDgwIGcPHmy5ujlzMkzT3+brzlp2CilfFq3bt3429/+xpEjR2otT05Opn379ixatAhwnKI8d+5cl+57xIgRjBs3jieeeKLmRb6kpISHHnqIhISERp+4YIxh1qxZrFixgo8//hgAm83Go48+SnJyMtHR0XTv3p3OnTsze/bsmkBdsGCBywO13uo6Rc2XPvTUZ+/ly+PTU59dZ/v27ZKWliYxMTHSoUMHGTx4cJ3r3XXXXbVOfRZxfB9SU1OlS5cuMnbsWPnyyy8FkA4dOsg999wjIiKZmZkSExMjgYGBkpaWJnv37pUJEybUWpaXl3fWsl27domISHl5uTzyyCPSsWNHSU1NlfT0dJk+fbrYbDYREdm1a5ekpaVJYGCgxMTESFpamqxZs6amxqefflqSk5MFkOTkZJk+fXrNY0uWLJH+/ftLly5dJDU1VR544AEpKSmpeXzLli0yfPhwSUhIkGHDhslDDz0k/fr1k/DwcElLS5OKioqa/7tTdX/99dfn/f9u7KnPRk57v9IXpaSkiGVJ7mZZWVlkZmZaXYbb+PL4rBzbtm3bzjrrydU84QQBd/HlscGFx3ehnx9jzPcictZ51fo2mlJKKbfTsFFKKeV2GjZKKaXcTsNGKaWU22nYKKWUcjsNG6WUUm6nYaOUUsrtPD5sjDH9jTGfGmO2GWM2G2O+M8Zcb3VdSiml6s+jw8YY0xFYBuQBvUWkN/AP4H1jzDgLS1NKKdUAHh02wGigBfB/ImIDEJG/AoXAT6wsTCmlVP15etjYnJ9rboVgHFOY+gH+dT5DKaWUx/HoudGMMS2Ab4GdwM1AKfA/wJPA1SKy/BzPuxu4GyA+Pr7f+++/3zwFN7Pi4uJ632zJG/ny+KwcW1RUFF26dHHrPqqrq2vdLMyXnBrb4cOHue6669izZw8ZGRksXbrU6tIu6Oc//znLly/n4MGDbN68maSkpLPWudD3bteuXRQUFJzz8eHDh9c5N5rlszJf6ANoC3wOlOPo3ewGhtX3+Trrs/fy5fHprM/uk5OTI35+fnLvvfe6Zftnjm3YsGFnzSZdHwsXLpQXX3zxrOXr16+XmJgY+fbbbxtb4nnNnj1bANm7d2+dj1/oe9fYWZ89+m00Y0wK8B2wH4gFWgGPAR8ZY66xsjallGeaM2cOfn5+zJ8/35LbH9fXokWL+OMf/3jW8vDwcJKSkpp8G2tP49FhAzwFRAO/EJFSEbGLyDzgK+AtY4ze1lopVcu8efP4wx/+QH5+PkuWLLG6nAbr1q0bGzZsoFevXlaX4lKeHja9gYMiUnbG8p1APNCp+UtSSnmqtWvX0r17d+6//35iYmJ46623aj0+evRoWrdujTGGtWvXctVVV9GxY0cGDhzI1q1ba6374YcfMmLECPr370+fPn3IzMxkzZo1593/3LlzadeuHcYYUlNTWbhwIQCvvvoqnTt3JjY2lhkzZnDVVVexZMkSDh06RHp6Ounp6cycOZPPP/+c9PR0jDFMmzat1razs7MZN24cSUlJpKWlcckllzBz5kyKi4sB2LNnD3feeSfp6elkZGSQnp7OrFmzqK6ubuL/qovU9d6ap3wA/wYKgIAzls8D7EDLC21Dezbey5fHpz0b97j//vtr/m8ffPBBCQgIkNzc3FrrPPnkkwLIgw8+KNXV1VJVVSVDhw6Vyy67rNZ6V111lfztb3+r+fqDDz6Q8PDws/7/zuzZrFy5UgBZuHBhrfUeffRReeWVV2q+vu222yQpKanOcQDy5JNP1ny9a9cuiY6OlocffljsdruIiHz00UdijJENGzaIiMh7770nw4cPl7KyMhEROXz4sHTt2lVmzZpVa9vas6nbyzius5nhPOUZY8xwYBIwX0TyrCxOKZ8zLercH+tm/3e9dbPPv+7p/jb03OstefC/6x3a0KTSKysr2bBhQ80dUO+9915sNhtz586tc/077rgDPz8/AgICGDduHN98802tHs/LL7/MT3/605qvJ0+eTEREBAsWLDhvHZdffjmdO3dmzpw5NctEhPfff58pU6Y0amzTpk2jurqap556CudLIRMnTmTw4MH4+Tlexq+66irmzZtHSEgIAK1bt2bSpEm8/vrrjdqnq3l0z0NEPjDGXA38Dsg2xlTjOKJ5DPiTpcUppTzKJ598wnXXXVfzdY8ePRgyZAhvvfUWDz/88Fnrd+vWrebfsbGxABw9epQOHToAEBISwn333cfatWux2+0YYzhx4gT79u07bx3GGG655RZmzpxJfn4+MTExZGVlkZqaSlxcXKPG9sUXX9CrVy/CwsJqLf/qq69q/h0ZGclf/vIX5s2bR0FBAQEBAeTm5pKfn9+ofbqaR4cNgIh8juPUZ6WUu0079/UTtfS/w/FxLkVF//33PV+de73Ttc2o33rn8M4777Bjx45afZr8/Hx+/PFHNm7cSHp6eq31T3/hPnV0cKq/UVJSwvDhw2nfvj3Lly8nJiYGgI4dO9brDLdbb72VGTNmMH/+fO69917eeustpk6d2uix5eXl0a9fv/Ou8/jjj/PSSy+xbNkyBg0aBDiOiKZPn97o/bqSp7+NppRSF3Ts2DHsdjtbt25l48aNNR9btmwhKCjorBMFLmTVqlXs3r2bBx98sCZoGqJLly4MHDiQOXPmUFpaSlZWFmPGjGnwdk5p2bLlBY9Q3n77bUaNGlUTNJ5Gw0Yp5fXmzp3LVVddddbyyMhIhgwZwty5c6mqqqr39k4dvZzqjwDY7XaOHTtW721MnTqV1atX8/zzzzNmzBiCgoJqPR4YGHjqhCdKSkrOe5r2qFGj2Lp1K2VltU/MveGGG8jKyqqp+fR6AXJzc+tdr7tp2CilvN7bb7/N2LFj63xs7NixHDt2jE8//bTe2xs0aBDR0dH8+c9/pry8HIBZs2ZRWlpa723ceOONBAUF8fTTT9f5FlqnTp3Iy8ujoqKC1atX89BDD51zW9OmTcPPz49p06bVBNQ777zDhg0buPTSSwEYM2YMX3zxBZs3bwZg586dzJ8/v971ul1dp6j50oee+uy9fHl8euqza+Tl5UmfPn3E399f0tLSZPv27bUef+mllyQ5OVkAadWqlYwePVoSEhIEkLS0NNm0aZP84Q9/kA4dOgggPXr0kAULFoiI4xTmAQMGSNu2bSUzM1OmT58u7dq1k+joaBkxYoTk5ORIWlqahIeHS3h4uKSlpcmRI0dq7X/ChAlyrtegI0eOSGZmpnTt2lV69eolixcvls8++0zS0tIEkISEBBk3blzN+lu3bpUxY8ZIYmKipKWlybXXXiu7d++uefzEiRMydepUSUhIkIEDB8oNN9wgU6dOrRnrqlWr5Kc//Wmtsb7++utn1eWuU589eiJOV0hJSZEdO3ZYXYZbZGVl1Zzm6Yt8eXxWjm3btm306NHDrfsoKioiMjLSrfuwSkPGNmPGDPz9/XnsscfcXJXrXGh8F/r5McbUORGnx5+NppRS3mrhwoUsWrTI6jI8gvZslFLKhYYNG0ZFRQUrV66kTZs2dU7jfzHSIxullHIhYwzdu3cnPj6+1iwCFzsNG6WUcqFTpyKr2vRtNKWUUm6nYaNUM8s/eZL1R2ys/E/9LxB0NV8/C1W5R1N+bvRtNKXczF5dzbZvP6d4/Qe0OrGOqOrjvFzxV2TDd7x15yUM6xbfrPUEBARgs9kIDAxs1v0q71dVVYW/v3+jnqtho5Sb7M1eS+7Kt+l0eCm9+O/dMCoIpF9YHutKWzHtw+/55IFBhEe0aLa6QkJCKC4ubtScX+riVlhY2OjrpzRslHKho4XlLN54iC1rl/NS0a9qbiWbSzx7215DdNpYOqcN4d7V3/Lxhj08WPACm98eysD/91qz1RgfH8+PP/5IcHAwoaGhZ82npdTpRISqqioKCwvJz88nMTGxUdvRsFGqiUqLC8he/i4Hd23i4WNjsQsYErg3pCOFsX2IvORmug8YRevT3n4I8DM8MDKFpI+O0OnI++xYN4WU/lc0S70hISEkJCSQm5tbr+nyG6O8vLzmJl6+xpfHBnWPz9/fn8jISBITEwkODm7UdjVslGqEapuN7NUfU/79u/Q6+W/6mwoyxPCc3xB6p3RjUt92dE75nuDAc/+KdUkbzDdrfsJlh98haOlDVPb+jqDg5nkRi4qKIioq6sIrNlJWVhYZGU27P42n8uWxgfvGp2GjVAPs3L2b41/MIjn3U3pzwrHQwPaAHhR0m8TSUVcTExNb7+1l3PocB19YRif7ftZ88DwDb37CTZUrZS0NG6Uu4OiJkyzecoIP1x/kWO5Bvg1+jwBjJ8ck8GP78XQYdjvdu6Q2atshYREc7n0f7Tc9QUDuRhdXrpTn0LBRqg6n+jBB2QuILdvPsxV/RPAjKrQlS9v+gm59BpEyYCTt/Jp+qVpQdBsAAqsKm7wtpTyVho1STnX1YQAq8ef2LmVcOnAow7vHExxwpUv3GxzhOAXZVFe6dLtKeRING3XR23a4kOVr1nL9D3fV0YeZTPcRU3kyLsFt+/fvMIDk8jl0jGzBMrftRSlradioi1Leof1sWP05/5fTg22HCzHYmRgMB/1ac6D9uCb1YRqqRXgI1fhTWG5rlv0pZQUNG3XROL0P06tsPcMx/E/Fq0SFtmRsnzbkpXxC7+49aO+CPkxDRIU6po0pLKtq1v0q1Zw0bJRPO18fZmvEZfzfNclc2r8/wQGNm+/JFYID/FgQNINYCigvG0pIaLhltSjlLho2yidtO1zIwg05rNywhU8qf4a/EUcfJrAnBV0n0X3EVDLc2IdpCGMMnf0OE0cBeQXHNWyUT9KwUT7j2KF97F4+G799X3Nj8UMIfkAYn0RcSau2SSRm3k73zr2sLrNOpSacOCmgpOAELVs3bu4ppTyZho3yaiVFJ8le/i4h2QvoWb6BeOO438bw0N20TRvBxIx29E0c7fGTTZb5R4ANyovyLryyUl7IK8LGGDMZ+AUQDsQAJ4CXRERv8H0RqrYL32bvJuhfv6NnwVcMqOnDBPBD+GWYtCn8ddh1zTbPmCtUBESCDSqK8q0uRSm38PiwMcb8ErgVGC8iB40xgcBbwAhAw+Yi8p+d2SzYZVi8MYdjhWWsDl5HmKmo1Yfp6yF9mIaqDGwB5VBZomGjfJNHh40xpiMwExgsIgcBRKTKGPNroK2FpalmcqoP02rPIjpWH2BBxavk04KkuAhWd3yKS/v19dg+TEPYghw3T6suPWlxJUq5h0eHDY4jmpMisvb0hSJyCDhkTUnK3UqKTpK9bC5Rm94jbsWWmj7MSRPB/alVZAwZRN/EaI/vwzTE4ZhLmH2knKjAjlaXopRbeHrYDAL2OXs2DwHxOPo1b4jIPyytTLlUtV1YtSuPpet28viOyQwwZYDzepjwy5A+U+g1bDJ3hYRZXKl7HG5/Dc9t6cQ9IZ2tLkUptzAiYnUN52SM2QJ0BDYDE4GjwGTgPeBJEXnmHM+7G7gbID4+vt/777/fLPU2t+LiYiIiIqwuo0mKcnchB7/liaIJnHTeNHJO4LPEBFSxM2oQ4V2vIDi8hbVFusGZ37sVP1bxVnYlw9oHcEdq4+6E6El84WfzXHx5bND08Q0fPvx7Eel/5nJPP7IJwXEG2m9EJNe5bIExZgrwqDHmRREpPfNJIvIa8BpASkqKZGZmNle9zSorKwtvHNvRnL3sWT6bhL2L6WTfB8BblT04HpfBxIz2JPX5hMRWMeR56fjq48zvXdm6HWzetpQeoW3IzJxiXWEu4q0/m/Xhy2MD943P08OmyPn5zLtKbQAmAT2Bdc1akWqUktJSsv81m5BtC+hVvpFWp/owRLCj5ZXMGJZJz9QMn+rDNES74i18FDyNLUczAO8PG6XO5Olhsx1IB86cGbHa+bl5Z0xUDVJdbWfV7uMs3JDDF1tyWOH3HPGmgEoC2OK8HqbXsOu41Iuuh3GXkBZxAARXF1tciVLu4elh8zGOP/P6AF+ftjwVKAO2WlGUOjex29mzZQ3HVr1N4pEveaD8GQpwvP+7MP4W0hJj6X7FrV57PYy7hEY6bqAWqmGjfJSnh818HGehPW2MGSsixcaYIcB1wAwRKbG2PHXK0Zy97Fn2DxL2LSbZvp9k5/IpLbYQesmtTMpoT2LcGEtr9GRhziObcP2RVj7Ko8NGRKqNMVcDzwFbjTHlQAXwcxF53drqVEmFjc83HaDLlz8ltXx97T5M3EiiLpvK7/oOxzTz/WG8UUS0I2wipRix2/X/TPkcjw4bABE5AfzM6jqUg62qkq3ffMbsnPZ8nn2Usqpq3g0sw+bnz5bwQZi0G7UP0wjBwaGUSjBhpoLS0kLCIqKtLkkpl/L4sFHWO70P0+XIp6Rxkv0V0ymTrgzoGMPxrjMpy+hJ39h4q0v1asUmnDAqKDp5XMNG+RwNG3VO/70eZlGtPsxB04bbMqJ5acRwEuN884p+K/xP+FPsOGHj736x6OkTytdo2KhaSipsfLYll0XrD/B/B29koCkA/ns9TNTAW0jpO5z22lNwuZMRnTl4PJ/CSs+d1UOpxtKwUdiqKsletYSK9fO5P38KR6sc/ZZFgUPp16JA+zDNpEWI49exsKzK4kqUcj0Nm4uU2O3s3vwNeavfpsuRz+iDY2r7K+yJ7O44mYkZ7RnT+0qiwgItrvTicWXFv7g+8N8E7/8Z9LjB6nKUcikNm4vM4ZOl7Pv4OVrvXUgX+366OJcfMG052GE8D1zxU9p17GZpjRerzrbdXOr/HWvyh1tdilIup2FzESguOslnO4pYuDGH1buPMy/wMzr57SefSHa2HFXTh+mgfRhL2YOjAJDyAosrUcr1NGx8lK2qkuyvl1C54T16FXzF25W/Z5MkE+Tvx5oOdxHcMZyeQydpH8aDmFBH2BgNG+WDNGx8yJnXw5zqw2BgUssDTBk8gTF92hAVqn0YT+Qf6ri2xq9Cw0b5HpeGjTHmzyLy/1y5TXVhhwvKWLzxEJd9dRtp1Ztrroc5YNpyMHE8SZl3cHun7pbWqC7MP9wxGWdgVaHFlSjleg0KG2PM1AusMroJtagGKC7Mp3DbF9y1s5ple8sRgRkBrUkM2MfOlqOIHngr3fpmah/GiwQ5wyaoqugCayrlfRp6ZPPmeR7TK9HczFZVydavF1O1fh69Cr9ivKlkZVUlgX5XMLJnKzr0fIbwnknah/FSQdFt+aa6J4f9OtHL6mKUcrGGhs02zj56CQe6Az8BXnVFUeq/RITdm1bXXA+TdlofZrNfCmMGpvH4yJHah/EBwW17clPV43TwC2WS1cUo5WINDZsHRWR/HcuzjTGfAvOAFU0vSx0uKGPRhkMs3HCQZ/N/zUC/nYCjD5OTOJ7EzDs4vj/Xp++FfrE59QdDQanOIKB8T4PCRkSWneexMmOMdqGboLgwn23L5xK6bQG/K7qeLfZOAHwUejXVcanEXDaVrhnDavowO/fnWlmucrGI4ADCKCe8ogR7tR0/f+23Kd/R0BMEhta1GIgBJgDlrijqYvLfPsx79CpcyQBTCcDkgLYkplzGpIz2DEu5hkB94fF5Af5+fBd8PxGmjMKiMbRw3lBNKV/Q0LfRsqj7RAADHARuaWpBFwMRYeuhQvI+mUHq4Q9q9WGyA1MpTpnMpBFTuSOmpbWFqmZXbMKJoIySgjwNG+VTGho2u4G7zlhWDRwFdotItUuq8lFHDu5iyY5yFmw6zs4jxTwTsJfMgJO1rofpqdfDXNTK/CLAnkdpwXGrS1HKpRoaNq+KyL/dUomPKi7MJ3vZO4Rt/4Ce5T+wreoedtqHEhMWyPGUu9jZ9UG6pg/V62EUAGUBEVAJ5UX5VpeilEs19ASBP7qrEF9yZh/mEmcfpoJAhrWpYvSI/gxLidc+jDpLRUALqITKEg0b5Vt0bjQXOdWH+Wh9Dlesv5/Bst7xgIHswN4UdZ9M9ytu5Vrtw6jzsAVGAlClYaN8jIZNE+Ue2MXe5bN542h3lh2PBaDKvw9JQbmOPszwO+nZMcXiKpW3sAe3cHwu1bBRvkXDphGKCk6wbflcwrZ9QM+KH2hthE22MWwIv4NxfdowMf33tO/wovZhVIP9p821vH4wkYEthjLQ6mKUciENm3qyVdvZ8vUSbOvePqsPszViEAP6XcedQ0doH0Y1SUV8Kl/a/Wgv8VaXopRLadich9jtbMkpYOHGwyz54RC/Lp/NlIAsRx8mqDfFKZNJueJW+mofRrlIixDHr2RhuU5Zo3yLhsTcQZsAABfnSURBVE0dTvVh2uxfzBvl41hsHwzAv2PGktSqB0nD79A+jHKLhOpcHvD/iLijSUC61eUo5TJeFzbGmJXAYKCTiOxz1XYdfZh3CN/2AT0qNtHaOCZKGBu0nuiMm5nYtz1p7aMwxrhql0qdJa76GL8K/IDsglTgCavLUcplvCpsjDGTcQSNS9iq7az8Tx5lXzzN8Lz3zurD+KffRObQSYwKCnbVLpU6r5BIxxmNodXFFleilGt5TdgYY4KAPwBLacIdQcVuZ9emVXyyu4q526rJK67kDn87owMryQ7qTUnKZFJG3ErfaO3DqOYX1sJxt85Qu4aN8i1eEzbA/cA6YCeNCBtHH+YftN2/hK72A4TYxpFnu4nO8eG06X0nh7o9oH0YZbnwKMcfOZGiYaN8i1eEjTEmFvgNMAi4vSHPrS4vZOuzQ+hRsbmmD5NPC1I6tGbx6Mvpo30Y5UEiIqOpFkO4KcdWVUlAYJDVJSlVb6WVtnM+5hVhg6NT+o6I7GtoMERXHaVXZSkVBLIl4nL8M26i15CJXKF9GOWB/Pz9KTBhRFFC0cnjxMS3sbokpc6rutpO9qollH//Lkvz259zPSNS1+1pPIcxpgvwFdBDRAqMMdOAJznP2WjGmLuBuwF6tAnrN/OxBwjuPJjg0Ihmqrp5FBcXExHhW2M6nS+P73xjS8h6iDB7CZsyniEyplUzV+YaF+v3zhfUd3wHiuysPmTjm0NVzJNH6O53gGx7Er2e2vy9iPQ/c31vOLJ5HpgpIgX1fYKIvAa8BpCSkiLj75/prtoslZWVRWZmptVluI0vj+98Yxuz6WW2HipkyYDL6dM+unkLc5GL9XvnC843vrxD+9m1fDbxexfzZOnP2S+tAVjYYgLDE8poP+x2eKpPnc/16LAxxgwBUoEbra5FqeYSFRoIQGHZud//Vqq5lBYXkL38XYKyF9CrbD0Dnb3vG4K/JSftASZltKNf0ugL9r49OmyAUYA/sPa0gbR2fl5qjKkEHhWRpVYUp5Q7tAhxho1OWaMsUm0XVu/Ow++fD5OR/zn9TQUAlfizKfwy6HMDdw27nuCQsHpv06PDRkSe4IzLqE/r2Yx25QwCSnmKKQWv81zwYnbsfAR6P2R1Oeoisnfrt3y4vYRHVi/jSGEFLwXmEeZfwfaAHhR0m0z3EVPJiEto1LY9OmyUuhiF+EOUKcVeesLqUtRF4JizD5OwZyGd7fuorPwlR+wDSIwN40T3hznY8wW6d0lt8n68JmyMMaOBZznjbTQR0dkKlU+RkCjH5/J6nxOjVIOUFDn6MMHbHH2Yy5x9mALCGRhbwt3XX0bfxBiXXoPoNWHj7Mtob0b5PL9QxxlofhWFFleifMmpPszC9Tncnn0nA8xuwNmHCRsIaVPoOew6Oq35jn5JsS7fv9eEjVIXC/8wR9j4V2rYqKbbs+Vbjn79FjOPX87GIsfPVoL/AIJCgyjoNomUK6aS0bL1BbbSdBo2SnmYgHDHZJxBlfo2mmqcY4f2sXv5bFrtWURn+z46A0OqKsiPu4UJ6e2YmP6/dIyPbNaaNGyU8jAhEY63MIJsOhmnqr/SShtb/vUmIVvepVfZeuJP68NsjxvFNYNu5+G+gy2bC1LDRikPExjfmf+tup7K0EQetboY5dGqbTZHH2ZjLp9tzeV5WcBY/++d18MMRPrcRK9hk7m0AdfDuIuGjVIeJjyuPa9UT6SVLVjDRtXpVB+mc+6nvFl5B8vs/QBY3eY64tuMbNL1MO6iYaOUh2kR6vi11BkE1Onq6sMAjA3LpvfAKUzMaEdSXLilNZ6Pho1SHiY00J+rAtYTaS+iomI4wcEhVpekLFJSYePzrblELfstmcVLz+rDRA28lQn9rsD4+Vlc6YVp2CjlYYwx/CHgNWIpJC//AYJbd7C6JNWMqm02slct4YMDLXh/h42yqmp+7h/GkAA/NoVfhvSZ4jF9mIbQsFHKA5WYCGKlkNLC46Bhc1HYs+Vbjq6cTfKRz+hNPp9W3UBZ9QT6J8WQmPoLyno+43F9mIbQsFHKA5X5R4ANygqPW12KcqNjOXvZvXw2CXsX1+rDHDSt6dsjma+uGU5inHcdwZyLho1SHqg8IBJsUFGUb3UpysVO9WEWbsjhJ/se5xr/7wA4SQQ74kYSNfBWUvpdQXsv6MM0hIaNUh6oKsBxdXdlqYaNLzjVh6n4/l3+cvISllX2AiA0IJOE8GBM2hR6DbuOS334ZBANG6U8kC2oBQDVJRo23mz35jUc+/rNmj4MwLVSTEHSYCb2bceY3qOIDguyuMrmoWGjlAeyBztvM1B20uJKVEPlFpSz7dO/0nHnbJLt+0h2Lj9oWnOgw3j6Zd7O+M69LK3RCho2Snmgzcl3c+fe4dwW35OBVhejLqik6CRfbs1lwZYCVu3O40G/TQwP3OfzfZiG0LBRygOFRURRRggF5TarS1HnUG2zsfXrxVSuf5eeBSvZZpvI19XjCfL342iX69nQ7gqf78M0hIaNUh6oRWggoFPWeBoRoTB3F2v+Mo8uRz6lD863OQ1cEnmcDsNTGdu7LVFhgdYW6oE0bJTyQG1LdzA/aAbFOV2AuVaXc9HLLShn8cYcWq2axsTKJTXLD5o2HOgwnsTMO7iicw8LK/R8GjZKeaCIQDvd/bazs8JudSkXrZKik2Qvm8vHhyKZcyAOEbjKrzPDAyPYEX8lUQNvIaXv8Iu6D9MQGjZKeaCwFnEAhFTrDdSak62qkuxVH9f0YQaYCg5WX848vwcY0aMVk9J+xvdHLmfEyCutLtXraNgo5YHCWjju1hkhGjbuJiLs3rqWvJX/OKsPsy2wJwm9r2btNSNr+jBZx3dZWK330rBRygNFRDmObCKkBLHbvWIKeW+Te7KMRT8cYuH6HMYc/zsPBiwCavdhemgfxmU0bJTyQCGh4ZRLICGmirKyEkLDI60uySeUFOaTvXwuods+4KOSVP5huwaA4NAruDTOT/swbqRho5SHKjbhhHCSopPHNWyaoK4+DIDdr5DDqXcwMaMdmSmtCAq4zeJKfZuGjVIe6t9BQ6koLebSSqGV1cV4GREh+3Ahu//1GoP2vXJGH6YXhd0m033EVP4SG29toRcRDRulPNTc6HtZX3iSD0yU1aV4jSMHd/PllhzmbIftuUWM8StkfNBJDpo2HOwwnsThd9Cjk/ZhrKBho5SH0lkE6qe4MJ9ty+cSum0BPct/IMQ+mO1V9xETGkBCnwns6JBJN+3DWE7DRikP1T6wiF5mLxX57QDvvR2wO9iqKsn+egmVG96jV8FXDDCVAFQSQJvoMF67uh+Z3eIJCvS3uFJ1iseHjTEmHbgf6Iuj3kDgS+ApETlmZW1KudO4/Dk8HbyQb/f/Di5Lt7ocy4kIWw8VsnBDDi3Xv8R99nmOB87owwzSPoxH8viwAeYBW4GhIlJijGkHLAOuNsakiUiZteUp5R72EEevxn6R39PmyMHd7F3+JlmHA/lrfj8Akk0/xoZ8xcHE8SRl3kGPTt0trlJdiDeEDcAjIlICICI5xpgXgDeA0cCHllamlJsYZ9iY8gKLK2l+Z/ZhEowQYe/I/LCBjEtry8SMQbRvfxcdtA/jNbwhbPqISOUZyw45P8c0dzFKNRe/0GgA/CsKLa6kediq7fzw7XJkzV9q9WEqCGRrxGX4pd3Et5kjtA/jpYyIWF1DgxljfgG8CHQXkZ11PH43cDdAfHx8v/fff7+ZK2wexcXFREREWF2G2/jy+OoztoKdX3HtoVmsCbiE8sGPNVNlrlHf753Y7eQUlrPysD/fHK4m07aKPwW9AsBmv+7sjRtGcOchBId6zkWtvvxzCU0f3/Dhw78Xkf5nLveGI5tajDH+wJ3A3+sKGgAReQ14DSAlJUUyMzObr8BmlJWVha+ODXx7fPUZ22a/AjgEEf6VDPSy/4cLje/Iwd3sWf4mbfYt4mhlZz633Q3AzrihrG5ZTcfM2+ndqTu9m6nehvDln0tw3/i8LmyA3wM24JdWF6KUO4VEOt4lDrEVWVyJazj6MO8Quu2Dmj4MgH9ABbf178CEvh1I7xCNMVdZXKlyB68KG2PMHcANQKaIzr2ufFtwmx5MrJhOcFQr5lldTCPZqu2s3JXH/qy3ufHQc7Wuh9kSMQi/tJvoOXQS04NDLK5UuZvXhI0x5lbgV8AVInLU6nqUcrfIyCg2SFeiKrzrfvZit1N4+D+8+e4hXtnTirziCrqZCG4PrmRbYC+KUq4j5Ypb6avXw1xUvCJsjDG3AI8AI0Uk17lsLNDW2Z9RyudEhjh+PYvKq7DbBT8/Y3FF53eqD9N63yLG239ko70zeZVP07llOOMyRnCoyzp6JHW1ukxlEY8PG2PMzcDrOHo1I42p+YUbAhy2qi6l3C3A348ngt8j2p5PSfEgIlt43pn+xYX5ZC97h7DttfswJySSilYZLBp3CWmJLTnt91ZdpDw+bICXgRDghToem97MtSjVrMaY1ST4Hyf35DGPCZtTfZiF63MI3raAF/xeBU67Hib9JxynJSNGjrK4UuVJPD5sRCTW6hqUskqpXyTYj1NaeMLSOsRuZ/fmb8hb/TbZRyuZUXYdAKH045YWfanoNo6UEVPpG9MScJw+q9TpPD5slLqYlflHgB3KLAqb3AO72Lt8Nm32L6aL/QBdgJ4Sxvy4mxjbtxMTMtrRIXayJbUp76Jho5QHqwyMhCqoLG6+sCkqr+LblV/QZu1MelRsorWzD5NPC3a2HEX0oNv4LH0IRuclUw2gYaOUB6sKdEzTYivNd+t+bFWVfLd5G/N22PlXdi6dbXtZGvxDrT5Mr6GTuDQo2K11KN+lYaOUB6sOcsz8XF3q+tsMOPowq8lb9TZdjn5OiL0lSypnABDRMYNVrV8gddikmj6MUk2hYaOUByuN7MiGI104Ka6biLKuPgxASUAE/zO4HaP7d6NDbBgwyGX7VErDRikPtqfTT3g6O4NbwhO5ugnbKa6w8enmw+xevYjfHv997T5M/JXEXDaVrulDuEf7MMpNNGyU8mC92jreRvtofQ73DkumfUxYvZ9rq6pk68pFbN65i6dzMiivshNOG+4JDmdvRD/8Mn5CryETtQ+jmoWGjVIe7LLkOEb3bs2JrcvZ8fd3aPfw2+c9C0zsdnZtWsXx1W/T9ejnpFFAkoQzo+rPXNIpgUkZ7fDrsYO+kb57PxblmTRslPJw06/uSOB//kh0cTHrPvkb/cffd9Y6hw7uZ/+yv9F2/xK62g9wagayH/3akZM4nuVXX0771jrxpbKOho1SHi4+Lo7v0n7HJT88TvL6ZzgxcDwxLdtw8PBhsg7YWLIxh9L96/lnsGPamBO04D+n9WEStQ+jPICGjVJeYMC197N5+4f0rtjA/r9ehc1eQoG9Bb+vfBaAkMDOLI+5ntjUkdqHUR5Jw0YpL2D8/Ii98S+UvjWMJPsBAPz9qhnTI4YRvRO5sldrIoJHW1ylUuemYaOUl2jXuQfbx71H4f5NJPQcSmJKBq/qW2TKS2jYKOVFuvcfAf1HWF2GUg2mfxYppZRyOw0bpZRSbqdho5RSyu00bJRSSrmdho1SSim307BRSinldho2Siml3E7DRimllNtp2CillHI7DRullFJup2GjlFLK7TRslFJKuZ2GjVJKKbfz+LAxxrQyxsw1xuxwfnxgjGlvdV1KKaXqz6PDxhgTBHwBBAG9gJ5ACbDCGBNhZW1KKaXqz6PDBrgN6AM8IiI2EakGHgE6A/dZWplSSql68/SwmQz8KCJ7Ti0QkVwg2/mYUkopL+DpYdMH2FvH8r1A72auRSmlVCN5+m2hWwLf17G8EAgzxoSKSNmZDxpj7gbudn5ZYYzZ4sYardQSyLO6CDfy5fH58tjAt8fny2ODpo8vqa6Fnh4252LO96CIvAa8BmCMWSci/Zulqmbmy2MD3x6fL48NfHt8vjw2cN/4PP1ttDwgso7lkUBpXUc1SimlPI+nh80moGMdyzsBm5u3FKWUUo3l6WHzEZBkjOl4aoExJgHoAXxYz2285vqyPIYvjw18e3y+PDbw7fH58tjATeMzIuKO7bqE86LOdcA24GbADvwdGAxkiEixheUppZSqJ48+shGRSmAUUI3j2pptQAvgCg0apZTyHh59ZKOUUso3ePSRjaqbMaaNMeYzY4z+paA8kjFmpTFGTu+3qoubT4aNL88UbYyZCHwDJFtdi6sZY9KNMa8bY743xvxgjMk2xvzJGBNvdW1NZYxJNsb8r3Ns3xtjdjpfkMdYXZurGWMm4+ir+gxjTEdjTLExZmMdH9FW1+cKxpjJxpivnD+fe4wx64wxt7pq+z4XNhfBTNG/w9HHWmV1IW4wD4gFhopIGo5xXgmsMsaEWlpZ010DTAFuFJF+QHccfzQsMcYMs7QyF3L+/v0BWGp1LW6wTkTS6/g4aXVhTWWM+SXwGPAT589nCrATGOGqffhc2OD7M0VfLiL/sboIN3pEREoARCQHeAHoCoy2tKqmywGmicguABGxA8/i+B281srCXOx+HGeQrrW6EFU/zrc6ZwL3iMhBABGpAn4NvOKq/XjrdDXnU+dM0caYUzNFv2BZZS4gIjara3CjPs4zEE93yPk5prmLcSURWVjH4hbOz8easxZ3McbEAr8BBgG3W1uNaoBbgZMiUusPBBE5xH9//5rMF49sdKZoL1VH0AB0AwT4qpnLcStjTDvgVWC987MveAJ4R0T2WV2ImyQYY94xxmxw9tzeNcb4wmvKIGCfs2ez0hiz3Riz2hhzpyt34oth0xIoqmN5zUzRzVyPaiRjjD9wJ/B3EdlpdT2u4DxRYBdwEPAHJohIocVlNZkxpgtwA/CM1bW4STVgA14G+gH9gSrgW2PMACsLc4EOOPrbvwaux9HnfhF4zRjzmKt24othcy7nnSlaeaTf4/gF/6XVhbiKiOwWkS5AFI4G7A/GGF84c+t5YKaIFFhdiDuIyAER6S0i34qI3fkHwr04Tj561uLymioECAd+IyK5zvEtABYDjxpjwlyxE18MG50p2gcYY+7A8ZfyNb44W4TzxeqXwBHgzxaX0yTGmCFAKvAXq2tpTs7Xks3AQKtraaJT7wRtPGP5BiAMx5FOk/li2OhM0V7OeW7/r3BMS3TU6npcwRgTaoypdXQtjuk7NgOpxphgaypziVE43hJce+raExx/9QMsdS7z6rMJjTFRztO6z1SNY+zebLvz85l5UH2O5Y3ii2HjipmilUWMMbfgOFV9pIjkOpeNdd591Zt9St1/AXfE0U+s6+QIryAiT4hI8unXngB/dT482rnM26+7eQnH2aw1nOHTG8dJHt7sY+fnPmcsTwXKgK2u2Ikvhs2bOP5afM4YE2CM8cNxDvleLrLDfG9jjLkZeB3H93CkMeYWZ/iMA9paWZuLTDfGxAEYhweAAcCfRCcp9Aa/Mca0gZqTV14A4oHpllbVdPNxXBf19KkL351vjV4HPHPqurem8smJOJ1HMi/iOGNEgC3AQyJywNLCXMAY8wKOty0ScVx78oPzoUvOceqw1zDGnODc19NMF5FpzViOSxljLgfuwhEuNhxN2eM4+jXv+krYON8uexZoDSTgmKm90nm047WcpzjfAwxxLmqJY2zPiMgKywpzEec1Us/hmLGjHKgAXhaR1122Dx/5GVdKKeXBfPFtNKWUUh5Gw0YppZTbadgopZRyOw0bpZRSbqdho5RSyu00bJRSSrmdho1SSim307BRSinldho2Siml3E7DRimllNtp2CillHI7DRulPIAxZq4xptAYYzfGfOlc9qoxJt8Ys9cYc5fVNSrVFDoRp1IewhhzPfA+8DMRecMYk4TjXiODfPFuperiomGjlAcxxnwEjATSgX8Az4rIv6ytSqmm07BRyoMYY1oD2ThuyfuJiNxhcUlKuYT2bJTyIM5bYU/HcXMur78pl1Kn6JGNUh7EeRvzLCAUx91Ye4lInqVFKeUCemSjlGf5BfAtMAHHraNfsrYcpVxDj2yU8hDGmGTgAxxnn5UZY+4B/gqME5FPrK1OqabRIxulPIAx5hnga6A1cKdz8f9zfp5rjPnAksKUchE9slFKKeV2emSjlFLK7TRslFJKuZ2GjVJKKbfTsFFKKeV2GjZKKaXcTsNGKaWU22nYKKWUcjsNG6WUUm6nYaOUUsrt/j8y1HMjIuIu7QAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig = pyplot.figure(figsize=(6.0, 4.0))\n", "pyplot.xlabel('x')\n", "pyplot.ylabel('u')\n", "pyplot.grid()\n", "u0_analytical = numpy.array([u_lamb(0.0, xi, nu) for xi in x])\n", "line1 = pyplot.plot(x, u0, label='Numerical',\n", " color='C0', linestyle='-', linewidth=2)[0]\n", "line2 = pyplot.plot(x, u0_analytical, label='Analytical',\n", " color='C1', linestyle='--', linewidth=2)[0]\n", "pyplot.legend()\n", "pyplot.xlim(0.0, L)\n", "pyplot.ylim(0.0, 10.0)\n", "fig.tight_layout()" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "def update_plot(n, u_hist, u_analytical):\n", " \"\"\"\n", " Update the lines y-data of the Matplotlib figure.\n", " \n", " Parameters\n", " ----------\n", " n : integer\n", " The time-step index.\n", " u_hist : list of numpy.ndarray objects\n", " The history of the numerical solution.\n", " u_analytical : list of numpy.ndarray objects\n", " The history of the analytical solution.\n", " \"\"\"\n", " fig.suptitle('Time step {:0>2}'.format(n))\n", " line1.set_ydata(u_hist[n])\n", " line2.set_ydata(u_analytical[n])" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [], "source": [ "# Create an animation.\n", "anim = animation.FuncAnimation(fig, update_plot,\n", " frames=nt, fargs=(u_hist, u_analytical),\n", " interval=100)" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Display the video.\n", "HTML(anim.to_html5_video())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Array Operation Speed Increase" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Coding up discretization schemes using array operations can be a bit of a pain. It requires much more mental effort on the front-end than using two nested `for` loops. So why do we do it? Because it's fast. Very, very fast.\n", "\n", "Here's what the Burgers code looks like using two nested `for` loops. It's easier to write out, plus we only have to add one \"special\" condition to implement the periodic boundaries. \n", "\n", "At the top of the cell, you'll see the decorator `%%timeit`.\n", "This is called a \"cell magic\". It runs the cell several times and returns the average execution time for the contained code. \n", "\n", "Let's see how long the nested `for` loops take to finish." ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "23.1 ms ± 344 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "%%timeit\n", "# Set initial conditions.\n", "u = numpy.array([u_lamb(t, x0, nu) for x0 in x])\n", "# Integrate in time using a nested for loop.\n", "for n in range(nt):\n", " un = u.copy()\n", " # Update all interior points and the left boundary point.\n", " for i in range(nx - 1):\n", " u[i] = (un[i] -\n", " un[i] * dt / dx *(un[i] - un[i - 1]) +\n", " nu * dt / dx**2 * (un[i + 1] - 2 * un[i] + un[i - 1]))\n", " # Update the right boundary.\n", " u[-1] = (un[-1] -\n", " un[-1] * dt / dx * (un[-1] - un[-2]) +\n", " nu * dt / dx**2 * (un[0]- 2 * un[-1] + un[-2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Less than 50 milliseconds. Not bad, really. \n", "\n", "Now let's look at the array operations code cell. Notice that we haven't changed anything, except we've added the `%%timeit` magic and we're also resetting the array `u` to its initial conditions. \n", "\n", "This takes longer to code and we have to add two special conditions to take care of the periodic boundaries. Was it worth it?" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2.52 ms ± 64.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n" ] } ], "source": [ "%%timeit\n", "# Set initial conditions.\n", "u = numpy.array([u_lamb(t, xi, nu) for xi in x])\n", "# Integrate in time using array operations.\n", "for n in range(nt):\n", " un = u.copy()\n", " # Update all interior points.\n", " u[1:-1] = (un[1:-1] -\n", " un[1:-1] * dt / dx * (un[1:-1] - un[:-2]) +\n", " nu * dt / dx**2 * (un[2:] - 2 * un[1:-1] + un[:-2]))\n", " # Update boundary points.\n", " u[0] = (un[0] -\n", " un[0] * dt / dx * (un[0] - un[-1]) +\n", " nu * dt / dx**2 * (un[1] - 2 * un[0] + un[-1]))\n", " u[-1] = (un[-1] -\n", " un[-1] * dt / dx * (un[-1] - un[-2]) +\n", " nu * dt / dx**2 * (un[0] - 2 * un[-1] + un[-2]))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Yes, it is absolutely worth it. That's a nine-fold speed increase. For this exercise, you probably won't miss the extra 40 milliseconds if you use the nested `for` loops, but what about a simulation that has to run through millions and millions of iterations? Then that little extra effort at the beginning will definitely pay off. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "###### The cell below loads the style of the notebook." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 29, "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 (MAE6286)", "language": "python", "name": "py36-mae6286" }, "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.9" } }, "nbformat": 4, "nbformat_minor": 1 }