{ "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", "\$$\n", "\\frac{\\partial u}{\\partial t} + u \\frac{\\partial u}{\\partial x} = \\nu \\frac{\\partial ^2u}{\\partial x^2}\n", "\$$\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", "\$$\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", "$$\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", "\$$\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", "$$" ] }, { "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", "\$$\n", "u(0) = u(2\\pi)\n", "\$$\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": [ "