{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Step 4: Burger's Equation\n", "\n", "Next, we combine what we have learned about convection and diffusion and apply it to the Burger's Equation. This equation looks like —and is— the direct combination of both of the PDE's we had been working on earlier.\n", "\n", "$$ \\frac{\\partial u}{\\partial t} + \\frac{\\partial u}{\\partial x} = \\nu \\frac{\\partial^2 u}{\\partial x^2} $$\n", "\n", "We can discretize it using the methods we have developed previously in steps 1-3. It will take forward difference for the time component, backward difference for space and our 2nd order combination method for hte second derivatives. This yields:\n", "\n", "$$ \\frac{u^{n+1}_i - u^n_i}{\\Delta t} + u_i^n \\frac{u^{n}_i - u^n_{i-1}}{\\Delta x} =\n", "\\nu \\frac{u^{n}_{i+1} -2u^n_i + u^n_{i-1}}{\\Delta x^2}\n", "$$\n", "\n", "Given that we have full initial conditions as before we can solve for our only unknown $u^{n+1}_i$ and iterate through the equation that follows:\n", "\n", "$$ u^{n+1}_i = u^n_i - u^n_i \\frac{\\Delta t}{\\Delta x} (u^n_i - u^n_{i-1}) + \\frac{\\nu \\Delta t}{\\Delta x^2}(u^{n}_{i+1} - 2u^n_i + u^n_{i-1})\n", "$$\n", "\n", "The above equation will now allow us to write a program to advance our solution in time and perform our simulation. As before, we need initial conditions, and we shall continue to use the one we obtained in the previous two steps.\n", "\n", "## Initial and Boundary Conditions\n", "\n", "The Burger's equation is way more interesting than the previous ones. To have a better feel for its properties it is helpful to use different initial and boundary conditions than what we have been using for the previous steps.\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", "This has an analytical solution, given by:\n", "\\begin{eqnarray}\n", "u &=& -\\frac{2 \\nu}{\\phi} \\frac{\\partial \\phi}{\\partial x} + 4 \\\\\\\n", "\\phi &=& \\exp \\bigg(\\frac{-(x-4t)^2}{4 \\nu (t+1)} \\bigg) + \\exp \\bigg(\\frac{-(x-4t -2 \\pi)^2}{4 \\nu(t+1)} \\bigg)\n", "\\end{eqnarray}\n", "\n", "Our boundary conditions will be:\n", "\n", "$$ u(0) = u(2 \\pi)\n", "$$\n", "\n", "This is a periodic boundary condition which we must be careful with.\n", "\n", "## Obtaining $ \\frac{\\partial \\phi}{\\partial x} $\n", "\n", "Evaluating this initial condition by hand would be relatively painful, to avoid this we can calculate the derivative using sympy. This is basically mathematica but can be used to output the results back into Python calculations. \n", "\n", "We shall start by loading all of the python libraries that we will need for hte project along with a fix to make sure sympy prints our functions in latex.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Adding inline command to make plots appear under comments\n", "import numpy as np\n", "import sympy\n", "import matplotlib.pyplot as plt\n", "import time, sys\n", "%matplotlib inline \n", "\n", "sympy.init_printing(use_latex =True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We shall start by defining the symbolic variables in our initial conditions and then typing out the full equation." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR0AAAAcCAYAAABcQWtqAAAABHNCSVQICAgIfAhkiAAAB3ZJREFUeJzt3WusHVUVwPEf9pYWWm0vFbFWQYgERYlF04daoInYiI/WYGp8Jif4SjTGajQRPxkqUWPVxpho8cEjJsZHVGwQIkVaDGJEa2urpgixBdqC4LNaCA/rhzWTM2fuzNxzDrf3nnPv/icnd87M2nuvM3uvNWvvmTWXRCKRmKbMxzac0EOZk3EAm0r7z8LaCdKrX87Erfgj9mBehczzsD2T+T3WF459BH/Ijn1ZnJeF+A12YS/eW5D/Ef6BH1S0U1UX7M/a3ZXpapw23oB9+DPeU9h/Tiaffx7Bmxr0apL/WKbrXryz9Duq9K0rM4obCjKX4Ov4Ll5tejLINjSQ/bEB7+uxzJVC6fIJ+yA+UVOmlX2ONztwQbZ9CkYqZBZjabb9bBwUzulU3IO5mIXb8Yps++RMfh7+gkXZ99V4o7FOp64uwojnl+Tr2hjBXViSldlXaLvIfDys7WTr9KqSPw87M11Pwq+EE8yp0repzFWF35ozii01ugw7g25DXfXH0/qouF/ejut7kD8bL8SNpf0XYSPejd+pjjDG43asyLa/KSKFXngxHscvsu9/xxMVcofFVRseEMZ3SvZ9RBjS7OzzVzyJo9nxOeKKll/VtuNIjT5VddVR18ZyEU0cxH/EeV9TUX4tbsF/u9CrLP8i3IFHRfSzG69tKGucMj/B20ryl+Nr49Q5rAy6DXXVH5PldObgNDzYQ5lNQuEyO0QI/hqcrz34e2Gj8PIfxf/wpR7Lny0Mc6u4Cn+yizIvF1HGfXhI/L57cUiEzPdkcguFYd2PzwtH1URTXcfE+boT7yiUqWrjOcLh5BwUUU+Zt4grZ7cU5feKyGihuAKuLrVRpW9TmZ14ZaH8RtwsDGm6MQw21FV/VE0JjgeLxLy/yK6a9tdgmQj179L5I3JOF6F4zon4dbadRxIbsr/L8Vip/E0i7Hy9sVfaJr0OZdsjYmq1VEQVNwlDubmiXK7TddrrJ6Ni/eT54up9Iy7EbfgnXioG2A/FtKVpoDXVtUo4j8XCGe0Rg62qjW54huiPt/Ypn685/Rz/ElOlJwvydfrWlXkok4XLxDTvVLHedlWXOg4Lw2BDXfXHZDmdR0X4X2RplWDGSjFQ14s5/mz8G1fgudrGn/NYob5W9veahvqXiRN7QEyTutUr56BYjL0v+/7TrFyV05mDH+Oz+GW272LcLaZlxALcSuEoch4U0cgFmp1CU1155HI40/Flwoir2jikM+pYoj0Ic9bhZ6I/u6FKfov2HP8bYtE6p07fujJzhaOFb2Wf6cow2NDA9ccB/Tm5ls5FsFfh++PItxqOLxGGdqYIB1/Sh04jImQcFVPUrSLaKHMCvoNPlfavzMrni783CAM9DU/PZBaIqcV5hXKrjXVAdXXNK9Q1H78VA6WujRFhzE0LyVvF1atMlV518s/K/p4jHEo+Jur0bSpzftbGTGHQbWjg+uM68WN7paXzhC0QA3IPzq2Rb9XUdZKINi7Mvq/X2/pEkUsyHfbii6Vj+eLxKjHfLd4+zp3IlfiTWLzNb3Mvz2R2C+N6f6HObSJ8PSrWYop3CarqOiurZ3em44cz2aY21opw/G5j75IsEJHRiaX9dXrVyd8hpll3inWunDp9m8pswAfMHAbdhgauP1aavncVElPDLSLanCkMug111R+zJkGRnPtFmLxrPMFEogtGRSQ1k8bTINvQTOyPRCIxE/i0eLai6bN6qpRLJPqgZXLH7YyzobqV8DeL3Jtl4lmLA7gan9P5XMVmfHucNu4VJy6RmCp6yVWaCLq1H5INmSVu8R4Tt0+/Kk7KvmzftVOn2lBRl2RXRTmJ76km4lUlYJaT8RL1tPQfXUyF/Qz9WPuKODmf0RkFzRa5FsdU32JLdFKXZFdFOYmvKRGP8Z+hWK06AbMqGS8xlpb+nc5U2M/QjbVi7tUKcY/9epGvUUxgfFzbS6+QaKIuyW6ReP6FGJC7xBWnmMQ3EYl421UnYFYl4yUmjqmwn6Eca0Vv/KFMsaPGPkFL+6nDycxMH0Y24ePG5rv8TTun5QmRs3W6ziS+PBGvpTMvZiLYqbpfZzL7cUbNsVsr9l2r/so/FfYzlGOt6HTy1xiMdzU88BQVms6s05xkd0Q8DXpEDALGJvGVE/HoPRmvimIyXiLYrPN9PkT+0TrhYPaXjjU9gzLZ9jO0Yy13OnNFNuhtIuxK9EdTkh2RZLcY7xILekd0JvFVJeLRezJeFcVkvESwuWJfSxj0NWL60A1TYT9DO9Zyp5OvZj+zx8oTnVyu/f6SlgiprygcPyzeL3OxeB3AUZHLMiLC4DMymePBC0R+VmLimQr7Gdqxls8vHxHzu3NxaY3sKpObNjEdOSye4bhU++19O7QXF/eK25h1iXjdsE1kEL9OZwLmRcYuOCYmhkG0n6EYa2vEKvsx8V6YL4i3gX1PvInu3j4VSzQzWUl8My05sl9a+rtlPgz2M5BjbZm45/6AOIEPC0+4xfR9w/4gcJnj+9TsqPZ/Y0g009L/czrDYD9prCUSiUQikUgkEolEIpFIJHrm/8ctltqIXYL3AAAAAElFTkSuQmCC\n", "text/latex": [ "$$e^{- \\frac{\\left(- 4 t + x - 6.28318530717959\\right)^{2}}{4 \\nu \\left(t + 1\\right)}} + e^{- \\frac{\\left(- 4 t + x\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}$$" ], "text/plain": [ " 2 2 \n", " -(-4⋅t + x - 6.28318530717959) -(-4⋅t + x) \n", " ──────────────────────────────── ─────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) \n", "ℯ + ℯ " ] }, "execution_count": 2, "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 * np.pi)**2 / (4 * nu * (t + 1))))\n", "phi" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwcAAAAoCAYAAACrSoYxAAAABHNCSVQICAgIfAhkiAAAF8xJREFUeJztnXncXUV5x79kJQiEEDWQEBZZIppAEIGwhcuWylK0GqniwhVcaAENFFuhVd4iEBDKIrQQWvCFurSihYJRCWBClYZNSCTIEpY3QhYghkggYQl5/eM38znnnnv2uyZ5vp/P+dx7z5mZM3PmuWdmnnnmGTAMwzAMwzAMwwAGdjoDLWAE8EYD8T8OfN59f67x7BiGYRiGYRiG0QwuA4amXN8UuB/4SOjcDTHhJgDbZ9zro8DR7vto4Es582gYhmEYhmEY6xtHAf8O/DdweIfzkouRwKUZYf4R+AbB4OBY4DfAVyPhqkAlcm4bYAawM3AtmkW5HhgMfBsYUi7bhmEYhmEYTWFz4C5gkwJxNgMWUduHeh9wXBPzVZadgNnA74FHgXdFro8F5rjrvwM+Gbl+BvCYu/5d9Fy2Ah4C5gELqFXu3gK8AvwkRzqePnfveS6vZNzjWOBJYCHwRXdunAvrjzXAxzLylBbnLJffBcBnM/KaFn4EMJN6RqA+cddzLHASsCfws8gxGJgEnIg6/n5wMAn4ciiN3VDHfzZwm/s+MnR9OvB9NAMBqrDpwPnA/k0uj2EYhmEYRhGmUduvycMFSBMcHhycipSpSVTd0WruAQ5237cGBkWubwtMdN+3ARYTDCDeAzyD+mwDgXtRX20gGhDhwj5H0NerAH9JbUc8KR1PHxqUhUm6xyDgKWCMi/Mktf1M3PnloXLE5SlKOM4E4GGX32HAfWiwkpTXtPAA11Hfx/0OsJf/Ea2UbmJLYCUwHw0UohwGvBf4ALAK+CWwBxpBeZ4CTkEC34dGo54t0Ah1LcEahbfQbMS6ppTAMAzDMAyjPCegjmRedgXeD9wOjHfnDkEWEX8E/ho4CHi9YD7uBc5EptzXI4305QXT+CDwNvBr93tFTJil7gBYhjrIW4fyO4hAoTsYeAl4B1jtzg1FswB+JmAO9ZYjSemkkXSPfZGGfrG79gtgCvCjUNzjgLtDZUjKU5hwnN2BuQR91flIKf5fCXGzwt8GfNqFAcnGncAjPoEBGZnrJAuBHVOuX4hG1D9E9lIgQfoKMhVKYxBwJXA2GjTsHbpmAwPDMAzDMDrNUGAU8GKBOJeivk2Ye5Di9EikHS46MAB1IL+BBgjrKD4wAA1cXkMDl4eBczLC74009s+73y+j8v0BWILMrZ5x17ZCneAXgEvQoCKJtHQA+tEzexD4TOh83D1GEwwMcN/HRO53PJrJKUI4zgI0mNgKmf9UQveIy2taeNCzP8B9PwkNPj9B8RmqjjAAjU7bxVhkTmQYhmEYhtFpRiNb8jDe5j16jEaOVS5x4arUmhU9G5P+EAL79j+4w/+OW3f5W6TNHhxzLS1fnqlo9mIsGvjMRgOWOLZGGvkDQudGAHe4a8OQBn5yJN4oNMsxKnSuQq0JT1Y6viO9rcvDHin3mApcHbr2dWTv79kSzUpsSi3RPIWJi/MV1KmfDdyIlONpeU0KD6q/8ICmjk6aFc1LuP8UNJJbB/wnEqA325CfUWg2wTC6icnoRbM3esl+AejtZIYMwzCMtvAG9Z3KiXEBHZOAT6FFvJujTuCryIvjkpjwb4XSq7rP3oS090Gd6UXINChKWr48i9GiXj8T8HMX785IuKHArcBFwP+Hzh8BPE1gjjQTlfn/QmFeRNr9g0nufGel4zvOS10eP0StyXr4Hkuo1cqPAR4I/f4oMItiLvbj4swgWDD8H8i6Ji2vSeFBMrWmQH4Mw+gyjkYmdFORvWO1o7kxDMMw2skiyilyqwQzBwcCN+cIX024NgZ1hndC2ujxCeGyGITs2kcg65DbqV9Tugmy1++JiT/JxfcLiWeijvQotI4UYDiasZgQilehdqCQlA5oAbBPa3M0W7JPyj0GoY530oLk24lfMxLNU5i4OO91n+NQ539QSl6Twnv2cvfoCGOAm9AU0krgp9RO82TR32WHYXSa17DBgWEYxsbETahzX5QqweBgOOo4PoqcuCSFr8acH4a0997s5pMUt58Pc5TLxwK0l5XHm08dhCxHwi49wx39C4DHkQmNd0G6rws3H3WEvxIKfxdaY7AarRXYPyUdkMvX+e5YAHzNnU+7x3HIAc7T1NrtD0ezDFETraQ8pcWZi9yuPkiwTjYpr0nhPdOAv6UD7IQKNx2tmp6IFkzcUjK9uM3OwoQ3S8u7CdqOLqxfvR3eBK0sNyI7sajfXqP57I0GbV/MCrgB0W2DA5N3wzA2ZLqhnZmE3LAbRrO4G83etJ07kClEmCOQ7VsZwpudHY9Gru9Bq+fDm6UV2QTN4wcHfhO0suyDRrtnNpBGFscgO7QXkL3Ys2iqMGlPhjPQi+2EFuYpiZHohXoLGk2vAf6E6udkmuMp6xZkZxf18buh0k2Dg26Ud+iczLdD3tvBVOAq5GrwVfQsv58Svtnl7iN59nZZRtzDXT6WoXVqS1BblKT0KRL+YtSgPo/KuAKZJZxLvU9z0P80azb6nYR8bUdgI/4meiZXkN6YFylL0TqO47ME5YjrOBd9XlCu7jdBm1Hdj96PryOb9lOIl70y8toN7cxJFNsEzTCSGEGwsVpb2QH9mVejP6s/1iATo6JENzs7EE2JDEcvTL9Zmg9bZBM0qPUTW3ZmA9SJWYmm4FrBxei5LkeLSy5C9mpvoU7aZ2PifN/FGdeiPKVxirv3EuAHaBbpBvSM+lHeG33Z7evSynKHtqHQTYODbpR36JzMt0Pe28E8lN9VaMo9q+PY7HL3ubg9McdZsTHEd9z9nkcb/FyIXFw/7K41Gv4ttJHQDUgWr0LT9f1oQeDYSPiJCWXoQZ3mfrShZ5Sd0ax7P8GCzF+5308Q37EuWpaidRxlLKqjVSQPDoo+LyhX9z9wab6IynwlMqXoR+Y4UcrIa7PbmfPJHjhWmnQvw+gajkMj8V1ijqhpTx7OQYOAWWjdwo5opuC7qANwAvBxF/bLaIAQpUr2zAHIjq6Mhm831GG5rkTcPGyDtEzLCBaZeA5FL5M4N2WPoxd42U5JlfIvqsPQgpro89wGuUvrR351G+VxtGBrfdHMNkK3DA66Vd6hMZmv0v3y3moORb7IN0HPIavj2Oxy97mjCF9y9+kl3v1i1O1i0fBQ7zHGc4FL699y5NMz18U5LubaHe7a6ZHzl7nzUfOSMmUpWsdhNiHwCX8JyYODMs+rj2J1/1cE74J3h84PQYst+wn6Bp6y8trMdubdaJOytGOzxNiGsZ5yFNp1uNlTcFU0czAEaUgq7vw+BKYNxwHfo34TtCr1Df5W6EX7NIFWouwim4vQS+XwkvGz2M+l/78J119FHaJofuKOzxW4b5XWaDHOceleFTk/i/iX8yao8etHZQtzrjv/F03OYzfSLYODbpP3cJ4akfkq7ZV3KCfz7aJCOZMTT1q5k+ijWAdxKFr3soj4znGj4bPYE5Ux6ooxiQku/AvIlDXMzu7ac9R3QrcgMJvxa3yaUZYKxer4a0gxMBlp9Iva46c9rz6K1f1NLq1TY65NdNd+VSC9NHnN0858Au2OuxzNnCx0aUbruQhZMwx22NGSoxX7HNwHvIL2KDgPNeTvQzZOp1N+B+Le0PfwlORv0VQhyHTotoy4npWheD7NhTHh8nAE0nTeVzJ+FgvRy2ZfpHUI7/w3GTUct4bOPYwWi56IvAyEX8RzWpTHIngfyWsj57+O8v5tVB5vk3spKst1aO1JmHvd55FI67ahsTmadQN1GLZHDd8KpOnqBN0m79DdMp8k71BO5tcX0sqdxlBkNrY96gz/Dvkfj7PRPxKtP7sCtS3HIDeLbyBf43MbDJ+Fdzf4u9RQAd7s9Xrqy3Oo+5xFfTu5Cr3rpqDZ8btpflmy2B0NVK9E9XFYiTSynleRut/GfcbNIvpzB6OB01s58pYmr2ntzEA0uPoUUjbejNZ9HIVmSsah/3IZ1gdTxGayGZqluZl0M0JQ23grqpN+1M8cT3wfMItb0ED5brQmxzMC1e0xJdI0YvgwGrF7u8T51G/n3UwqBN6KyvJh9KItyrvQy+TRBu+fxTTUALyEOgzTgR+jhmAW9eYXX0Z/mEa2w666NCoNpBFlEHpW/cRrYXrdtar77bU5SSZfw931B2KubQhUiB/Z93YoP90q79C4zFdpv7xDcZlvFxWXjzIzB3nKHUcf8fL+LHBITPh/dtenh+4XPu6h9r1eNHyUs5DG/HK0oLcftW952o5hSHG2lnibe2+m83cJ8a921/+mSWWB/HU8CC30fZJgnVGPi5s2c1DkefXFlCGt7n/orse5ZPQzB/3ITCeLLHlNa2d8vUyn1pf8YDSo6CfZfahRywXovXdpVkDURoTf9aeSrEipkj7rXkED17h9B64j3QmGYcSyG/rzz2rDvT6GNMbhF+dC4j2zXOuuf7iB+1VpfmfpUpfmzITrY9Hi9eeA01zYX5I+bb6GbE8mRnPoVnmHxmW+SvvlHcrJfDuoUH5wkKfccZyLNNKjkBZxPKrXdcjJxZ6R8Ne4+6xFWuaDkEZxAoH9/pwGwkdZRq08/oL8+/ec6OLELUQGdULSOtveXt8r2hotC+Sv4/OQ9j7cSerJyC8Ue15F6/4zLs2n0c69nsHIJNHfM0/HLo+8xrUz+7n8RWcyPV5h8YUcedjY2RWtK60SDA5GovoFDbzmEcymPEAgS4egtmIh8ooVda9dJdskt0L84OBYtMZ1o6aPYnZJZW1RNyT2J9DyJdFH48/171EjcBmaPtsMbZPtG4GoV4oH0FRq3hmVonnszZlumK+6uI9T+zKPMj10n3vJXpy1mHymC300Xg8bO90q71BM5ovmsTdHmlHyyjsUl3lovTxXSsYrUu68+M5b1KPcDHf+DeSsIsxmaH1auINYNHwSo9CC2CeR55sP5SiD1yLH7ZYKxQcHzShLhew63g/9F6P/uZ6M/IYp87w8SXU/EA2i+1GnfQYyeXoMdRQXuWv7ZaSfV17j2hnvHe2HxHtZ+om7fnJGHgwN6HajdnAAqkvPLKRMGYrqN8wc6v8HnirlBwej0azZRkV0zcEz6EWTlyVNzMv6yhr3meSdARp/rhXk2vEWav3KP4xeuE+hqehr0RTsIKQ9+j2yfczDFWiRdpiJaHO4G6lfKDaPYpxG4GbucGr/8FFeDn0/GWmN0hhGUA9pmHw3TjfKOxSX+W6Sdygu89Cd8ly03Hm5FtX55Mj5le7zEerrbDUaTJ6M1q/MLRE+iReRfD6M5PEmpOlO4oPAAWgh8s8TwvzJfQ5PuO7Pr4x8NlqWNAahsj0FfLNkGlD8eYVJqvt30EDrTLRO4UT0f5iDFgf7jt5LKWkXkde4dmaK+/x0WgGo78gatXwUycVT6H8SZhWS/VVoLcpKNKPwSiTc9tT+D4YQmIH5Qd8097kv+dahgN7N2+YMa0Qoor3aUA7PaPf7Nw08vyy85iTq3s7zP+6693iyh/sdt1t0EaounUqD6Uxz6TxKvK14mBPQNO1SF+eajPADXPhnGsxjUTotf52S+W6Ud2iOzFdpv7xDcZlvFxWUn7wzB0XLXQRv8x0dDJ3kzv8iIZ634/9GyfB5eMTFeXdKmCtdmJ6UMF90YWYkXPezZt5LWDPKUiG9jrci//vhipT7hMnzvMIk1X0amyIlwcspYYrIa1w7s6mLf0+BfBnxTEczXX3IAcWfgG+5a3PRupHT0PoOUGf/iVD87Uhvk6qUnznYgsC0aaOhWd6KNrYV9WGWohdQKzdd8mYSSYvL/Hk/Ep7oPh9pWY7y8w/Iw8U85FVgeUrYo5H5xgLUAP4aNZhXoOnoOMYh+Suq2W2UjVXmu1HeoXtkvoi8QzmZ70aKlrsofv+aqGcav6HYBwg6cGG8dvq5kuHzMNp9Ju14vClyp/sO8lKUxGz3OSUmb1ugDUBXE3gJa0VZorxJcp4/BOyFOmVPkn92Iut5RUmq+zQ+hTTHP0q4XlRe49oZ3wbkHeQYyZxNYC5XRbJ7nvu9FDgeecn7iDu3As3kDEKmXju4cK1gF2RyZnSYEWgKMrpXQae4mWTvER5vV7hLRriyHO/SXwaMiVw7CjUKawh2zzzDhT+xwftWaUyT+k0X/yGybY8PQg3fswRTeFNd/KTFXqCFXv1Iq7A+YvJeT1F5h+bIfJX2yTuUl/l2USHfzEHRcoPk/f3UbtK1O/ULCUF2xAvdPeJ2qfWLT8+InJ+CZOUVak11iobfjXhTnwEE6wDujbnu+ZwLc3tKGE/RTdCKliVKhWKzQ2F6XNzomoMyz6ts3W8Zc24iUmCsIBiIhCkjr0ntzHx3PrrZmucgGtvnYGOkSu2ag39Fzzk6CLsJDZhB8vZbNAsU5xmqSvrMwV1IZlYj07/wOp1pxHvEMkpwNvrDXJ0VMIZL0EZm0XOt8GE/GfnEXYzyW40JMwG9ZNJesJ928eM2Y2kGA5Df9n60b8SNyCb7NtQA9KPNaTyT3bkX0LPrAT5Z4r5VyneWTnRx1yI3dj0xR9WFnYjsCJdS30l+0KVzcMJ9fkSya8B2Ulbm4+Tdn2+2zG+o8g7Nkfkq7ZF3aEzmW8nH0ExGL8Fiz2dC56IuBouW29Pn4u0YOteD7Ipnoh10L0YD0TUu7EziPThtR7C77V2o/n/i8vQ29RvMFQ0/zeXhTgK3ujeg59KP6jDNVaV34Zm0EDnMzkhZ4AeI05Fb8H6knR8ZCV+0LFC8jpPoIX5wUOZ59VCu7u9Hawyudve5FZX7VeLdn5aV16R2Zoq7X78r77+4dH/sytupvWg2BiZRP1huBXcjJZ7RIJPQNOZ8ineUNkOajoMi5+9BmwTlpZd0207P0cCFSFu3muSR5UOkd4SGoBf6/XkzWILB6KV7H3rxrUULrX5GsCgqzGmoMXkDvbguKHHPKuU7Sz1k26jOQdrnZaje94hJ5wgXNm7DreGo8ei0lrWszCfJOxST+V5M3qFxma/SenmHxmS+1fSQXoa+guHD5Q7TR/3g4BDUCXsCDZzeRtq8O4HPk27K9x60s+0iZG62HC2A3bcJ4cej//U8F24tsol+EJU/Tfu8Oyrn8+TXII9FCoOlLm+LkJlZUgelaNl7KFbHSfh0ooODMs+rbN1/HWmMVyITqGeRpnm7jDwXkdesdmYfNJBZ5vK9HGmwZ9C6XeQNcRKtNfEdgQbTRoMMR6PlQwlG82Fuo9YF4j6ok+KnE6ciraWvbL+zYfiP+/sc+eglX2cpzGskd5a+RfYCTK853qvgfY3ynI6eeVznul2kyXxReYdyMt+LybthGEYr6IZ2xjDaSrN337wOjaBnJ1xfTK0N8YNoi/Qj3O+DkRag3/1eS2D7tR+yyz2Q9vMA0sIMSwlzOZpCPC8ljNE8hqEO6k9preecLNJkvqi8Q3fIvMm7YRhG97QzhtFWmjk4+BKaLv+nlDCLqZ/ue43A+8gO1PriXoc6R6tQx8pPxbebJcjMIW5xk+cNtPDsIeIXVhnNZUfUMT+rg3nIkvmi8g7dIfMm74ZhGN3RzhjGess46t0bzqHerOgLyC7Qm1HsijpDPt4d1Pt4/iZa0JXGOajT5Y+3kWlG+FzW4r40M4tdkXY376YtxoZPHpkvI++QLfMm74ZhGIZhtIRm7XOwP3Iz9Vjo3EDkQeQUpFl8E2lShyDN6ctokdXtBP68l1O/6Goi2b7Lr0XeATwXu3t9N3Rucb6ixOIXUKVtqGJsXOSR+TLyDtkyb/JuGIZhGEZLaNbg4FZkXhDme8g/8YUEmxX5Dst2yBvKB9EiTc8j1Gsz9yR5B0jPCmq3PV/lfjdrV7vxKO8vNik9Y/0nj8yXkXfIlnmTd8MwDMMwWkKzBgcr3RHmddRhWRA65ztL5yN/6odQq528A2lBRwJ/DOXx/cj+eXXMfRphc4KNnAYA2yOt7Qpq/RMfTGv2WTDWX/LIfBl5h9bJvMm7YRiGYRgdYw7xPt9fR9rVnRLizaXWx/pn0MZG64Brct67l3yuHSvE+znuDYXZFPlonoRhpDOHepkvKu9QXOZ7MXk3DMMwDGMD5SPIJrtbthw/FZjV6UwYGywm74ZhGIZhdA3d0iEJ8zQyeViKNJidZi+kvf1jVkDDKIHJu2EYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYhmEYXc2fAa7lURfgXv5aAAAAAElFTkSuQmCC\n", "text/latex": [ "$$- \\frac{e^{- \\frac{\\left(- 4 t + x\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}}{4 \\nu \\left(t + 1\\right)} \\left(- 8 t + 2 x\\right) - \\frac{1}{4 \\nu \\left(t + 1\\right)} \\left(- 8 t + 2 x - 12.5663706143592\\right) e^{- \\frac{\\left(- 4 t + x - 6.28318530717959\\right)^{2}}{4 \\nu \\left(t + 1\\right)}}$$" ], "text/plain": [ " 2 \n", " -(-4⋅t + x) -(-4⋅t + x - \n", " ───────────── ─────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν\n", " (-8⋅t + 2⋅x)⋅ℯ (-8⋅t + 2⋅x - 12.5663706143592)⋅ℯ \n", "- ─────────────────────────── - ──────────────────────────────────────────────\n", " 4⋅ν⋅(t + 1) 4⋅ν⋅(t + 1) \n", "\n", " 2 \n", "6.28318530717959) \n", "───────────────────\n", "⋅(t + 1) \n", " \n", "───────────────────\n", " " ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "phiprime = phi.diff(x)\n", "phiprime" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In python code:" ] }, { "cell_type": "code", "execution_count": 4, "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 - 12.5663706143592)*exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))\n" ] } ], "source": [ "print(phiprime)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lambdifying\n", "\n", "Now that we have the expression for $ \\frac{\\partial \\phi}{\\partial x} $ we can finish writing the full initial condition equation and then translating it into a usable python expression. To do this we use the lambdify function which takes a sympy simbolic equation and turns it into a callable function." ] }, { "cell_type": "code", "execution_count": 5, "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 - 12.5663706143592)*exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)))/(exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1))) + exp(-(-4*t + x)**2/(4*nu*(t + 1)))) + 4\n" ] } ], "source": [ "u = -2 * nu * (phiprime / phi) + 4\n", "print(u)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "3.4917066420644494\n" ] } ], "source": [ "ufunc = sympy.utilities.lambdify((t,x,nu), u)\n", "print(ufunc(1,4,3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pretty neat right?!\n", "\n", "## Solving the Burger's Equation\n", "\n", "Now that we can set up the initial conditions we can finish up the problem. We can generate the plot of intiial conditions using the lambifyied function." ] }, { "cell_type": "code", "execution_count": 7, "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": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "#New initial conditions\n", "grid_length = 2\n", "grid_points = 101\n", "nt = 150\n", "dx = grid_length * np.pi / (grid_points - 1) \n", "nu = .07 \n", "dt = dx * nu #Dynamically scaling dt based on grid size to ensure convergence\n", "\n", "#Initiallizing the array containing the shape of our initial conditions\n", "x = np.linspace(0,2 * np.pi, grid_points)\n", "un = np.empty(grid_points)\n", "t = 0\n", "\n", "u = np.asarray([ufunc(t,x0,nu) for x0 in x])\n", "u" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5oAAAJmCAYAAAA90v9dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3X145XddJ/z3N5mHTNtkoNWSQktbOlSM+ECllQIKjB3ErZV191Yr1gVx1VZHl66uyl2loqss6671kbK4uwW3stzi5YqwFqUgD5bq4A7IFhRKO+WhzEDpwGTaTuYpv/uPk9OeZM7J/Gbyy5yTc16v68o1zck533yTGa7wyeepVFUVAAAAaMpYvy8AAADAcBFoAgAA0CiBJgAAAI0SaAIAANAogSYAAACNEmgCAADQKIEmAAAAjRJoAgAA0CiBJgAAAI0SaALAECmlVKWUX+73PQAYbQJNgBFUSnnZQkDS+fbFUspfl1K+s9/366cu35fOt9f3+35JUkr5Z2s1mCylPLuU8sullMfVfP5LSimvaPgOY6WUnyul7CqlzJVSPlpK+YEmPwfAqFvX7wsA0FevSrIrSUnyhCQvS/IXpZSrqqp6Rz8v1mfvSvKHXR7/5Km+SA//LMlPJvnlLh/blOTIKb3NiXl2khuTvDHJV2o8/yVJnp7ktxq8w68l+YUkf5DkQ0lenOTNpZSqqqq3NPh5AEaWQBNgtN1WVdXft98ppfy3JF9I8gNJGgk0SymnV1X1cBNn1fhcY0k2VFU1t8KjPllV1a1N3OlUa+BrH2qllCcl+Zkkv19V1faFx/5rkvcl+Y1SylurqjrazzsCDAOlswB0+kqSA+nIiJVSnr9QNvr8zieWUi5YePxlHY+9sZTyUCnlolLKX5RS9if5o46P/2Qp5d5SyoFSyo5SyreWUt5bSnnvkrM3llJeXUr5VCnlYCnls6WU/1hK2bjkeVUp5fdKKT9YSvlYkoNJXrTwsatLKf+nlLK/lDJbSvm/pZR/09Q3auFz/Fgp5Z7lvp6OMuULlrz2mO/rwuvfWkr5TMfXfVMpZVPHc96YVjZzUZnvku/JLy/5XM8opdy28H14qJTy7lLKs5Y8p33P55RSfrOU8kAp5eFSyv8qpXx1je/FNyz8/d+7UI66p5Ty30spZ3U855eT/MbCu7s67n9BlyOz8H28Msn5Hc+973h3OY4XJ1mf5HXtB6qqqpLcnOTcJJev8HwAIqMJMOo2l1K+Kq3S2bOT/FSSM5KsJJu3LslfJvmbJD+b5JEkKaVcl+T3knwgyU1JLkjyZ0m+nORz7RcvZCX/PMlzk7whyT8m+fok1ye5OMk/X/L5tib5voWzv5TkvlLKtiT/M8m7k/z8wvO+Nslzkvx2ja9hYuH7stRsVVWHFu75I0n+S5IPplXW+ZSFe+9N8tkan6Ob701yWlpBz4NJLkvr7+TchY9l4XM+Mcm2JD90vANLKV+X1vd8Nsl/THI4yY8neW8p5XlVVf3dkpf8blp/J69O6+/oFWl9b7//OJ9qW1rfg1uS7EnydUl+LMnXlVKetRDM/Wlaf4c/kNbf55cWXvtAjzN/LcnmtL7+6xcee6jja+v2d9TN/qqqDi789zOSPJzWv6tOOzo+/jc1zwWgB4EmwGi7fcn7B5O8vKqqd63gzI1J3lpV1SvbD5RSNiT51bT64bZWVXVk4fGPptWr97mO178kyRVJnldV1d90nHFXkteXUp5dVdUHO57/NUm+vqqqj3c897fSCqy+4yTLIH9k4W2pH0jyllLK+iS/nuQjSV7QEXx+PK3g+GQDzZ+vqupAx/tvKKV8Ksmvl1KeXFXVZ6qqurOU8skk22qW9/77tDJ4z62q6t6Fe/5hkk+kFXg+b8nzH0zywoXAsB34/3QpZXNVVfuW+Tyvq6rqP3c+UEr527QC/ucm+UBVVR8tpexM6/v4Z1VV3bfcxauqelcp5f4kj+/xtfYKUJf64bT+nSXJOUm+0P76Ouxe+POJNc8EYBkCTYDR9pN5bMDNE5Jck+S/llL2V1X1pys49+Yl7z8zyVlJXtkOMhf8UVrZzU7fm1a26Z+WZKzes/DnC9LKIra9rzPIXPCVJKenlWV754lfP29LK4u31P9d+POZaWWAX9UOMhe8MY+Vhp6wziCzlHJ6WoN9PphWxvkZST5zIueVUsaTvDCtoO7ejs+zu5Ty5iQ/WkqZqqpqtuNlb1gShH0grWzi+Uk+WvPuE2llxv924aFLFs5p2raaz/tYx39vSusXKkvNdXwcgBUSaAKMth1LhgH9zyQfTvJ7pZR3LAmi6jqSxRnKpBWkJMmnOh+squpIl567p6ZV5torW3X2kvd3dXnO69Iqp71tISP2V0n+uKqqukHn56qqWprt7dT+eu7ufLCqqsOllHu7PL+WUsqTk/xKku9O8vglH958Ekd+dVqluJ/o8rF/TGtWw3lZHIgtDWa/vPDn0vssUko5M61pslfn2L+jk7n7cR3n76iXA2ll3Zea6Pg4ACsk0ATgUVVVzZdS/jrJv0kr4PtYkqUlhm3jPR4/WFXV/AquMZZW5vDf9vj40rLUYwKDqqq+WEr5piTfkeQ7F95+uJTyh1VVvXQFdzsZtb5/C9nHdyU5M8lrk/xTWr2ET0orU3qqBvj1KjUux3ndH6e1uuQ30iopfiitO78zq3T3Usp0zafu68i47k7yglJau0w6nnPOwp+fb+yCACNMoAnAUu2fDWcs/NnOaD1uyfPOT32fXvhzS5K/bj9YSlmX1sCZzpLMe5J8Y5J3d+mjq20hG/v2JG9f6DN8XZIfL6X8alVVn1r+1cfV/nqemsdKerPQu3lhkn/oeG7d79/XpzUo56VVVT26w3NhsNFSdb8vD6Q1jOlrunzsaUnmc/L9pI8qpTw+ybcnubGqql/pePypXZ5+on+nyz1/9zIf69TZo/mRJP86rax5Z8n1t3R8HIAVst4EgEctBEovTHIoj03l/HRaWa5vW/L0nziBo/8+rSEzP7oQXLb9YI4tyfzjtLJ4P9rlfpsWeheX1blSI2llavNYMNutbPJE/X1aQdy1C4OO2l6WYwPKexb+fPT7t5C9/LElz2tnEkvH80pa2eWlHl74+NLPtcjCIKS/SvLizhUipZQnpDV06W+W9GeerGPuvuAVXZ7b3qm67N2XPL9X6e22mm9/2fGat6U1effRf78L3+drk9yfxf2/AJwkGU2A0fadpZSnLfz32WkFH09N8h/aAUhVVftKKW9N8lML+xrvSfJdObYPr6eqqg4t7FD83STvKaX8cVqZzJctnNeZtfofafVXvr6U8oIkd6RVZvq0hce/I61Abzn/daFn8D1p9Yuen9aakI/k2LUW3VxcSrmmy+NfqKrqXQu9mL+Y1qqR95RS/r+0Mpk/nGRRj2ZVVR9bmL76moU77U2rj3Hpz+B/Sut78Z9KKU9Ka2ruv0z33sj/s/Dn75RS/jLJ0aqq3tLja/nFtIKtvymlvC6tHtofTyvg/rne34L6qqqaLaW8P8nPLfyy4v60fmFx4TJ3/7VSylvSCvreXlXVw12e237+95dSfjOtqcUPVVX19oXPe8I9mlVVfW5hKvG/W7jrh9JamfOtSX7wJKcUA7BUVVXevHnz5m3E3tIK8KolbwfSGgR0bZKy5PlfleRP0sou7U3y+rT2JFZJXtbxvDemFQj0+rw/leS+tCZ8/l1aPX1/n+S2Jc9bn1YQdNfCc/cuPO9VSaY6nlcl+b0un+dfppXF+kJaE0Y/vXDn6Rrfm6Xfl8639y557nVpBZZzaQUs35rkvV2e95S0+i/n0tox+WtprXCpkjy/43lfu/C8/WllTN+Q5Bu6fJ/Hk/xOki+mVf5aLbn/Ly/5/M9Iq1dy/8Lf4XuSXN7j38Qzlzz+/KX37PF9e1JaezK/nNbU3z9Oq++x231+Ma1fABxd+PgFy5x7elrTib+88Nz7Gvj3P5bklQv/Fg8u/Dv7wX7/79KbN2/ehumtVNVJt78AwIos9E4+kORPq6o6plR2LSqlvDdJqqp6fn9vAgD9o0cTgFOilDKx0AvX6V+lNWX1vaf+RgDAatGjCcCp8qwkNy30ez6Y5JIkP5JW2eJb+3kxAKBZAk0ATpX70lql8dNpZTH3JvnDJL9QtVaRAABDoq89mqWUb0vy75J8c1oDA76nqqo/6/h4SfLqtEbcPy6tyYPXVVV1dx+uCwAAQA397tE8Pa2l1j/Z4+M/l9Zvvq9Na5Hyw0n+spQycWquBwAAwIkamKmzC7vZHs1oLmQzP5/kP1dV9Z8WHtuc1qj6l1W994UBAADQR4Pco3lhkukkjy5jrlpLw/8uyeVJugaapZSNaS2h7tTuBQIAAODETCb5fHUCWcpBDjSnF/78wpLHv9DxsW5emeTGVbkRAADAaDo3yf11nzzIgebJek2S3+x4fzLJ5z772c9mamqqT1cCAABYe2ZnZ3Peeeclyf4Ted0gB5p7Fv58QpLdHY8/IclHer2oqqqDSQ6232/vBp+amhJoAgAAnAL9njq7nF1pBZvf3n6glDKV1vTZO/t1KQAAAJbX14xmKeWMJFs6HrqwlPJNSfZWVfWZUspvJfnFUsrdaQWev5rWJNo/O/Y0AAAABkG/S2efmeSvO95v91a+KcnLkvzHtHZtviHJ45L8TZIXVVU1dwrvCAAAwAkYmD2aq2Wh3Hbfvn379GgCAACcgNnZ2WzevDlJNldVNVv3dYPcowkAAMAaJNAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABol0AQAAKBRAk0AAAAaJdAEAACgUQJNAAAAGiXQBAAAoFECTQAAABo10IFmKWW8lPKrpZRdpZQDpZR7Sim/VEop/b4bAAAA3a3r9wWO4+eTXJfkpUk+luSZSW5Jsi/J7/TxXgAAAPQw6IHms5O8raqq/73w/n2llB9Iclkf7wQAAMAyBrp0NskHk3x7KeXiJCmlfGOS5ya5rdcLSikbSylT7bckk6fmqgAAACSDn9H8D0mmkvxTKeVokvEkN1RV9UfLvOaVSW48FZcDAADgWIOe0fy+JD+Y5CVJLkmrV/NnSykvXeY1r0myuePt3NW+JAAAAI8Z9IzmbyT5D1VVvWXh/f9bSjk/razlm7q9oKqqg0kOtt83oBYAAODUGvSM5mlJ5pc8djSDf28AAICRNegZzbcnuaGU8pm01ps8I8m/TfLf+3orAAAAehr0QPOnkvxqktclOTvJ55P8lyS/0s9LAQAA0NtAB5pVVe1P8oqFNwAAANYAvY4AAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRKoAkAAECjBJoAAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRKoAkAAECjBJoAAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRKoAkAAECjBJoAAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRKoAkAAECjBJoAAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRKoAkAAECjBJoAAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRKoAkAAECjBJoAAAA0SqAJAABAowSaAAAANEqgCQAAQKMEmgAAADRqXb8vAABr0dH5Kjt27c0X98/l7MmJXHbhmRkfK/2+FgAMBIEmAJygd961O69++8eze9/co4+ds3kiN141kxc9/Zw+3gwABoPSWQDo4uh8lTvveTBv+8j9ufOeB3N0vkrSCjKvu3XnoiAzSfbsm8t1t+7MO+/a3Y/rAsBAkdEEgCV6ZSx/6cqvza/+739M1eU1VZKS5NVv/3i2zUwrowVgpMloAkCH5TKWP/HmDx/zeKcqye59c9mxa+8q3xIABptAEwAWHJ2v8uq3f7xnxrKuL+7vHYwCwChQOgvAyOk1MXbHrr3LZizrOntyooFbAsDaJdAEYKR067+c3jyRa77lyfmbu7+04vPP2dwKXAFglAk0ARgZ7f7LpWWwe/bN5T/91SdP6KyS7uW0r/quGYOAABh5ejQBGDrdVpMs13/ZabkQsaSVsXzdSy7J9Obu5bEHDh892WsDwNCQ0QRgqPRaTfLCr3tCrf7L7S+4KL/31/ckWZyxbAegN141kxc9/Zx8x9OnH+3z/PxXDuS17/xEkuTX/+KfcsXMEzI1sb6pLwkA1hwZTQCGRq/VJLv3zeVNH/x0rTO2PGEyN19zbMZyevNEbr7mkrzo6eckScbHSi6/6Ky8+JuelOuevyXf8XVPSJJ86aGDueldJ1aGCwDDRkYTgDWl18TYuqWxx3P25EQuv+isbJuZ7vp5evml75rJ+z75QOYOz+dNH7wv3/vN52XmiVMrvA0ArE0CTQDWjF5lsTdeNZN1Y2MrWk1S0spatifGtjOWdZ37+NPyU1ufmt/4y09kvkpe9ba78tZrL08pBgMBMHpKVa30d7+DrZQylWTfvn37MjXlN8sAa1WvibFtYyWZr/kjbenE2HYo2FkaezIOHjmaF/3WB7LrSw8nSX7j//mGnPv402pnRQFg0MzOzmbz5s1Jsrmqqtm6rxNoAjDwjs5Xee5r37OijGXb9VdcnLd86DNds6IrCTLb3vfJB/LS/74jybHBb5OfBwBOhZMNNJXOAjAwevVf7ti1t1aQuXHdWA4eme/6sXZp7PatW7J965YT6r88Ec+7+KvzjPMelw9/9ivHZFj37JvLdbfuXHHmFAAGnUATgIHQq//y2uddlA/c/UCtM37wW56cW+64L0nv1STtgPJE+i9PxNH5Kp/7yoGuH6sW7vLqt38822amldECMLSsNwGg75ZbS3Ljn38st//jF2uds21mutZqktW0Y9fePLD/YM+PV2l9XTt27V31uwBAv8hoAnDKdCuNTbLitSSdE2PHx8oJryZp0hf31+sjrfs8AFiLBJoAnBK9SmP/+TOeVKv/8se/7cK84f27khy/LPZEV5M06ezJieM/6QSeBwBrkdJZAFbdcqWxN7/3nlpnzDxxc9/LYuu47MIzc87mifTKn5a0Aux2NhcAhpGMJgCN6DUx9uh8teLS2KSVAbz8orP6WhZbx/hYyY1XzeS6W3ce87Fu2VcAGEYCTQBWrFdZ7I1XzeSMjetWtP+ys/8y6W9ZbF0vevo5ufmaS/Izb/2HPHzw6KOPT9ujCcCIEGgCsCLtstilGcvd++Zy7a07s2l9/S6NkuP3X64VL3r6Odn1pYfz2nd+Iknyiiuemp/a+tQ193UAwMnQowlALUfnq9x5z4N520fuz533PJij81WtstgDh+drnX/9FRcPfP/lidq8acOj/33O5glBJgAjQ0YTgOPqVRp79aXn1SqLXT9ecvho93C0XRq7feuWbN+6ZaD7L0/U5MRjP2b3zx3p400A4NQSaAKwrF6lsXv2zeWm2++udcYPPev83HLHfUmOXxo76P2XJ6Iz0Jw9cLiPNwGAU0vpLAA9LVcaeyJTZLfNTK+J1SRNm5xY/+h/z8poAjBCZDQB6LmaZMeuvY1NjB0fKwO/mqRpU0pnARhRAk2AEder//KXrpzJR+//Su1z6kyMXQurSZrUmdHcP6d0FoDRIdAEGGHLrSb5iTfvrH3O9VdcnLd86DOLglU7I5OpTTKaAIwmgSbACOhWGpvkuKtJjmeYJ8Y2YdP68YyPlRydr7L/oIwmAKNDoAkw5Fa6muS7v+GcvP2ju5OM1sTYJpRSMjmxLl955LCMJgAjxdRZgCHWLo1dGlDuPoHVJN8+84SRnBjblPaKE+tNABglMpoAa1yvibHLrSY5EWdPTuTyi84auYmxTZncuD7JgeyfO5KqqlKK7xkAw0+gCbCG9SqLvfGqmWzetKGx1STJ6E2MbUo7o3lkvsrc4fls2jDe5xsBwOoTaAKsUb0mxu7ZN5drb92ZC77qtNpn1VlNwsmZ2rR4xYlAE4BRoEcTYA1ariy2/dh9X3qk1lnXX3Gx/stV1M5oJsmsgUAAjAgZTYAB160Hc8euvbXKYpdmKpd+zGqS1Tc1sTijCQCjQKAJMMB69WDW7ZV86eXn5013fjqJ1ST90pnRtOIEgFGhdBZgQC23muRPd95f64zvePo5VpP02eLSWRlNAEaDjCZAH63WapLOibHjY8Vqkj6aXFQ6K6MJwGgQaAL0yXKrScbHSu3VJHUmxlpN0j+LS2dlNAEYDQJNgD7otZpk98JqkrrJxpc/54LcdteeRUHp9EKwqix2MEzJaAIwggSaAKuoW2lskuOWxc7XrJndNjOdG66cURY7wAwDAmAUCTQBVkmv0tirLz2vVlnsxLqxzB2Z7/qxpT2YymIHV2ePpmFAAIwKU2cBVsFyE2Nvuv3uWme85FuenJLHei7buvVgMrimOqfOHpDRBGA0CDQBGrbSibFt22amrSYZAounzspoAjAalM4CnKReq0l27Npbe2JsN1aTDJeJ9WNZN1ZyZL7SownAyBBoApyEXv2XP/vCr8md9z5Y+xyrSYZfKSWTE+vy5UcOZ/9BGU0ARoNAE+AELbea5Gfe+g+1z7n+iovzlg99xmqSETC1aX0r0JTRBGBECDQBejjZ1STH0y6N3b51S7Zv3aIsdgS0V5zsnzuSqqpSir9jAIbbwAeapZQnJXltku9MclqSTyX54aqq/r6vFwOGWq/S2O975rm1+i9/6FlPzq1/+5kkxy+NVRY7/CY3tgYCHZ2vcuDw0Zy2YeB//ALAigz01NlSyuOT3JHkcFqB5kySn0ny5X7eCxhuy60m+e13f6rWGc+84EwTY3nUpBUnAIyYQf+V6s8n+WxVVT/c8diufl0GGB69JsY2tZrk7MmJXH7RWSbGkuTYFSdLfwEBAMNm0APN707yl6WUtyZ5XpL7k7yuqqo/6O+1gLWsV1nsjVfNZPOmDY2tJklMjKVlUUbTQCAARsCgB5pPSXJdkt9M8utJLk3yO6WUQ1VVvanbC0opG5Ns7HhoctVvCawZvSbG7tk3l2tv3ZmzTt9Q+6w6q0kgaU2dbds/Z8UJAMNvoHs007rfzqqq/t+qqj5cVdUbkvxBkmuXec0rk+zrePvc6l8TWAuWK4ttP/bgw4dqnXX9FRfrv6S2qY6MphUnAIyCQc9o7k7y8SWP/WOSf7nMa16TVga0bTKCTRg53Xowd+zaW6sstt2r2Y3VJJyMSYEmACNm0APNO5J8zZLHLk7y6V4vqKrqYJKD7fftKoPR060Hc3pqIl9/7lSt17/08vNzyx33JbGahGYsHQYEAMNu0Etnb0ryrFLK/1tK2VJKeUmSH0vy+32+FzCgeq0m2TM7l3d9/Iu1ztg2M201CY1aPAxIoAnA8BvojGZVVR8qpXxPWuWwr0prtckrqqr6o/7eDOin1VpN0jkxdnysWE1CYxZnNJXOAjD8BjrQTJKqqt6R5B39vgcwGJZbTbLvwOHaq0nqTIy1moSm6NEEYNQMfKAJ0NZrNcnuhdUkdb38ORfktrv2LO7hXAhWlcWyGqb0aAIwYgSawMDpVhqbZEVlsZ22zUznhitnlMVyyizu0ZTRBGD4CTSBgdKrNPbqS8+rVRZ72obxHDh0tGtAurQHU1ksp8rE+vFsGB/LoaPzSmcBGAmDPnUWGCG9Jsbu3jeXm26/u9YZV196XpLHei7buvVgwqnUzmoqnQVgFAg0gYGw0omxbVaTMKjagebsAYEmAMNP6SxwSvVaTbJj197aE2O7sZqEQddecfLQwSOpqiql+PcIwPASaAKnTK/+yx/91qfkA3c/UPscq0lYi6Y2tX7kzlfJw4eO5oyNfgQDMLyUzgKnxHL9l7/yjo/nrz9RL9C8/oqLlcWyJk1utOIEgNHh16lAo1ZrNUm7NHb71i3ZvnWLsljWnM4VJ/vnjuSczX28DACsMoEm0JhepbEv/sYn1uq/vO55T8nr33dvkuOXxiqLZa1p92gmMpoADD+ls0AjliuNff377611xtPOmTIxlqHVmdGcPWCXJgDDTUYTqK3XxNimVpOcPTmRyy86y8RYhtKiQFNGE4AhJ9AEaulVFnvjVTM5Y+O6xlaTJCbGMpymFpXOymgCMNwEmsBxtctil2Ysd++by7W37szE+vpV+HVWk8Awaq83SQSaAAw/PZrAsuqUxc4dnq91ltUkjDLDgAAYJTKawKO69WDu2LW3VlnshvGSQ0e7h6NWk8Cx600AYJgJNIEk3XswnzC1MReedXqt11/zrPNzyx33JbGaBLqR0QRglAg0gZ49mF+YPZgvzB6sdca2melcduGZxwSr0wsDg5TGMuoWT52V0QRguAk0YUSs1mqSzomx42PFahLoYXHprIwmAMNNoAkjoNdqkhuu/NrseuDh2qtJ6kyMtZoEutu4bjwb1o3l0JF5PZoADD3gntApAAAgAElEQVSBJgy55VaTbH/zh2uf8/LnXJDb7tqjLBZWYGpifb700EGBJgBDT6AJQ2C1ymI7bZuZzg1XziiLhRWYmliXLz10MLNKZwEYcgJNWON6lcXeeNVMNm/aUKssdnLjujx08EjXgHRpD6ayWDh57T7Nhw4eyfx8lTG/qAFgSAk0YQ3rVRa7Z99crr11Z848fX3X1y31vc88N7fccV+tHkzg5LVXnFRV8vChI4tWngDAMBnr9wWAk7NcWWz7sb0P1yvP2zYznZuvuSTTmycWPT69eSI3X3OJHkxoiBUnAIwKGU1YA7r1YO7YtbdWWWy7V7Mbq0ng1Dp2xcmm/l0GAFaRQBMGXLcezOmpiXz9k6Zqvf6ll5+fW+64L4nVJNBvUx2lsibPAjDMBJowwHr2YM7OZc9svd2X22amc9mFZx4brFpNAqfc5KJA0+RZAIaXQBP6bLVWkyiLhcGzuHRWRhOA4SXQhD5abjXJVx45XKsHM0mtabHKYqH/DAMCYFQINKFPepXF7l5YTVLXy59zQW67a4+yWFgDOktnZw8onQVgeAk0YZV1K41NsqKy2E7bZqZzw5UzymJhDZhSOgvAiBBowirqVRp79aXn1SqLPW3DeA4cOto1IF3ag6ksFgafYUAAjIqxfl8AhlW7NHZpQLl731xuuv3uWmdcfel5SR7ruWzr1oMJDL6pTTKaAIwGgSasgpVOjG3bNjOdm6+5JNObJxY9Pr15Ijdfc4keTFhjZDQBGBVKZ2EFeq0m2bFrb+2Jsd1YTQLDyXoTAEaFQBNOUq/+y3/9rU/J+z/5QO1zrCaB0bF+fCwT68cyd3heoAnAUFM6Cydhuf7LX33Hx/O+moHm9VdcrCwWRky7fHZW6SwAQ0xGE5axWqtJ2qWx27duyfatW5TFwgiZnFiXB/YflNEEYKgJNKGHXqWxL/7GJ9bqv7zueU/J6993b5Ljl8Yqi4XR0c5oPnTwSI7OV36xBMBQUjoLXSxXGvv6999b64ynnTNlYixwjKmOgUAPHZTVBGA4yWgysnpNjG1qNcnZkxO5/KKzTIwFFplasuJk86b1yzwbANYmgSYjqVdZ7I1XzeT0DesaW02SmBgLLGbFCQCjQKDJyGmXxS7NWO7eN5drb92ZjevqV5TXWU0C0EmgCcAo0KPJSKlTFnvwyHyts6wmAU7GZEfp7OwBK04AGE4nldEspbxquY9XVfUrJ3cdaE63Hswdu/bWKovdMF5y6Gj3cNRqEmAlFmU0Dwo0ARhOJ1s6+z1L3l+f5MIkR5Lck0SgSV9168F8wuTGXPhVp9d6/TXPOj+33HFfEqtJgGYtHgakdBaA4XRSgWZVVc9Y+lgpZSrJG5P8rxXeCVakVw/mF/YfzBf2H6x1xraZ6Vx24ZnHBKvTCwODlMYCJ0uPJgCjoLFhQFVVzZZSbkzy9iT/o6lzoZvVWk3SOTF2fKxYTQI0blGP5pzSWQCGU9NTZzcvvMGq6bWa5IYrvzb3fPHh2qtJ6kyMtZoEaJqMJgCj4GSHAf300oeSnJPkh5LcttJLQS/LrSbZ/uYP1z7n5c+5ILfdtUdZLHDKTZk6C8AIONmM5vVL3p9P8kCSNyV5zYpuBOleGptkRWWxnbbNTOeGK2eUxQKnnIwmAKPgZIcBXdj0RaCtV2ns1ZeeV6ssdnJiXR6aO9I1IF3ag6ksFjjVzlgUaMpoAjCcxvp9AejULo1dGlDu3jeXm26/u9YZ3/vN5yZ5rOeyrVsPJsCptn58LKdtGE8iownA8BJocsodna9y5z0P5m0fuT933vNgjs5Xjz7eRGnstpnp3HzNJZnePLHo8enNE7n5mkv0YAJ91y6fFWgCMKyanjoLy+pVFnvjVTPZvGlD7Ymx3VhNAqwVkxPr84XZg0pnARhaAk1OmV4TY/fsm8u1t+7MeY/fVPssq0mAtayd0Xz40NEcOTqfdeMKjAAYLn6ycUosVxbbfuyzXz5Q66zrr7hYWSywpk12rDh56KDyWQCGj4wmjeu2mmTHrr21ymJLSaoeTZrt0tjtW7dk+9YtymKBNWvpipPHnbahj7cBgOYJNGlUrx7M9h7M43nZsy/IG++4L8nxS2OVxQJr1VRHRnNWnyYAQ0jpLI1ZbjXJ2z7y+VpnvNDEWGAETC3JaALAsJHR5IR0K4sdHysrXk1iYiwwSpaWzgLAsBFoUttyq0mqpPZqEhNjgVHXOQzIihMAhpFAk1p6rSbZvbCapG6u8eXPuSC33bVnUVA6vRCsKosFRoWMJgDDTqDJIt1KY5Mctyy2bsnstpnp3HDljLJYYKR1ZjRnD8hoAjB8BJo8qldp7NWXnlerLHZi/VjmDs93/djSHkxlscAoW5TRtEcTgCFk6ixJlp8Ye9Ptd9c64yWXPTklOaaMtlsPJsAom9KjCcCQE2iy4omxbdusJgGopTOjOatHE4AhpHR2hPRaTbJj197aE2O7sZoE4MQszmgKNAEYPgLNEdGr//L6bRfnjk99qfY5VpMArNwZi6bOKp0FYPgINEfAcqtJfu5PPlr7nOuvuDhv+dBnrCYBWKHxsZLTN4zn4UNHTZ0FYCgJNIfIya4mOZ52aez2rVuyfesWZbEADZicWJ+HDx1VOgvAUBJoDolepbHf+81PqtV/+bJnn583ffDTSY5fGqssFmDlJifWZc+sHk0AhpOps0NgudUkv/Oee2qd8YwnP97EWIBTaGpTayDQgcNHc/ho9x3EALBWyWiuEb0mxja1muTsyYlcftFZJsYCnCKdK04emjuSx5++oY+3AYBmCTTXgF5lsTdeNZPNm9Y3tpokMTEW4FSZXLLiRKAJwDARaA64XhNj9+yby7W37szjNq3v+rpu6qwmAeDU6MxozlpxAsCQ0aM5wJYri20/9pWaY/Gvv+Ji/ZcAA0SgCcAwk9EcEN16MHfs2lurLHbdWMmR+e5dmlaTAAymqSWlswAwTASaA6BbD+b01ES+7omTtV7/ry4/P7fccV8Sq0kA1oqpjoymQBOAYSPQ7LOePZizc9kzW2/Iz7aZ6Vx24ZnHBqsLA4OUxgIMnsXDgJTOAjBcBJqnwGqtJumcGDs+VqwmAVhDJmU0ARhiAs1V1ms1yau+ayZ7Hz5UezVJnYmxVpMArB0ymgAMM4HmKupVFrt731yu+6Odtc95+XMuyG137VEWCzBEZDQBGGYCzRVarbLYTttmpnPDlTPKYgGGiPUmAAwzgeYK9CqLvfGqmWzetKFWWezpG8bzyKGjXQPSpT2YymIBhsek9SYADDGB5knqOS1231yuvXVnznv8plrnfP+l5+WWO+6r1YMJwPCY3LgupSRVlcwKNAEYMmP9vsBatFxZbPuxz375QK2zts1M5+ZrLsn05olFj09vnsjN11yiBxNgSI2NlZyxofX7XsOAABg2MprH0a0Hc8euvbXKYsdKMt+jSdNqEgAmJ9Zl/8EjSmcBGDoCzWX06sG87IIza73+pc++IG+8474kVpMAcKzJifXJvjkZTQCGjtLZHto9mEszl7v3zeVt//D5Wme8UFksAMtoT56dOzyfQ0fm+3wbAGjOSGc0V2s1ibJYAOpYvEvzcM46Y2MfbwMAzRnZQHO51SRVUqsHM0mtabHKYgHoZumKE4EmAMNiJAPNXqtJdi+sJqmba3z5cy7IbXftWRSUTi8Eq8piATieqU2dGU0DgQAYHiMTaO64d29e8A2TSXLcsti6JbPbZqZzw5UzymIBOCmnb3zsx/Df3ftgZp445WcIAAPj6HyVHffuPanXrqlAs5TyC0lek+S3q6p6xYm89uVv+lCedPY9ufrS82qVxU6sH8vc4e6DGZb2YCqLBeBEvfOu3Xnz337m0ff//V/8Y/7bHbtUxQAwENqthvd/8eQCzTUzdbaUcmmSH0/y0ZM9Y/e+udx0+921nvuSy56ckhxTRtutBxMATkS7hWP/wcXlsnv2zeW6W3fmnXft7tPNAKD3Bo4TsSYCzVLKGUn+KMmPJvnyqfic26wmAWAVLDfZvP3Yq9/+8RydP9nZ5wBw8la6gaNtrZTO/n6S/11V1e2llF9c7omllI1JOsf2TZ7IJ7KaBIDVtGPX3mV/Q9yefL5j116tGQCsqm7rHo/3c6qugQ80SylXJ7kkyaU1X/LKJDfWOjtWkwBwan1xf70f3nWfBwAno9u6x+mpiXzTeY9r5PyBLp0tpZyX5LeT/GBVVXV/4r4myeaOt3O7Pen6Ky5WFgvAKXf25MTxn3QCzwOAE9WrB3PP7Fze+bE9jXyOQc9ofnOSs5PsLOXRctXxJN9WStmeZGNVVUc7X1BV1cEkB9vvd7yu9X5aAeX2rVuyfesWZbEAnFKXXXhmztk8kT375rr2v3S2cADAyepWFjs+VhrrwTyeQQ80353k65c8dkuSf0ry2qVB5vF0K41VFgvAqTQ+VnLjVTO57tadtVo4AOBEdSuLPWfzRG68aiaPHDpauwdzJT+JBjrQrKpqf5K7Oh8rpTyc5MGqqu7q/qrephe+uUpjAeinFz39nNx8zSX55T//ePbMdvTG+DkFwAq1y2KXZix375vLtbfurH3Oy59zQW67a0/u/+IjJ3WPUlVra3x6KeW9ST5SVdUraj5/Ksm+d314V17wDef7DTEAA+PofJWn3/jOHDg8n+mpidzxC1v9nAKglm6lsUny3Ne+p5Gpsf/zR5+Vyy48M3/90U9n2zMuTJLNVVXN1n39QGc0u6mq6vkn87rLnqL/EoDBMj5WcubpG3P/Vw7kyPy8n1MA1NKrNPbqS8+rFWSetmE8jxzq3oW4dN3jZU85uZkBAz11FgCG3eRE63e+s3NH+nwTANaCXhNjd++by023313rjKsvPS8lx/ZgNjkrQKAJAH00NbE+SXLoyHwOHjmhGXcAjJimJsZum5nOzddcsqrrHtdc6SwADJN2RjNJ9s8dycYzxvt4GwAGQa/VJDt27V1R/+XSsthtM9Ortu5RoAkAfbQ00PyqMzb28TYA9Fuv/svrnndRPnD3A7XPqbNCa3ysrNq6R4EmAPTR5ELpbJLMHjjcx5sA0G/LrSZ51Z9/rPY5119xcd7yoc8sClZP9QotgSYA9NHSjCYAw6/XapKV9l+2S2O3b92S7Vu3rFpZbB0CTQDoo6lNj2U098/JaAIMu16lsd/zjCfV6r/88W+7MG94/64kxy+NXa2y2DpMnQWAPpLRBBgdy60med1776l1xswTN6/6xNgmyGgCQB8t6tGU0QRY83pNjG1qNcnZkxO5/KKzVnVibBMEmgDQRzKaAMOjV1nsjVfNZHJifWOrSZLVnRjbBIEmAPTRlEATYCgsNzH22lt35rT19fck11lNMuj0aAJAHymdBVj76pTFPnL4aK2zrr/i4oHvv6xDRhMA+mhqwtRZgLWkWw/mjl17a5XFrh8vOXy0ezg6SKtJmiDQBIA+0qMJsHZ068F8wtTGbPnqM2q9/oeedX5uueO+JIO9mqQJAk0A6KPTNow/Oo1QoAkwuHr1YH5h9mC+MHuw1hnbZqZz2YVnHhOsTi8MDFpLpbHHI9AEgD4qpeSMjeuy78BhpbMAfbZaq0k6J8aOj5WBX03SBIEmAPTZ5EQ70JTRBOiXXqtJfunKmXx678O1V5PUmRg76KtJmiDQBIA+a02ePZDZucOpqiqlDNdvtQEG3XKrSX7izTtrn/Py51yQ2+7aM/RlsXUINAGgz9oDgQ4frXLwyHwmTmDXGgD1dSuNTbKisthO22amc8OVM0NfFluHQBMA+mxqyS5NgSZA83qVxl596Xm1ymLP2LguDx880jUgXdqDOexlsXWM9fsCADDqpqw4AVhV7dLYpQHl7n1zuen2u2ud8X3PPDfJYz2Xbd16MBFoAkDf2aUJsHpWOjG2bdvMdG6+5pJMb55Y9Pj05oncfM0lI9eDeTxKZwGgzyY7SmetOAE4Ob1Wk+zYtbf2xNhuRnE1SRMEmgDQZzKaACvTrf9yemoi33fpuXn/J79U+xyrSZqjdBYA+qwzozl7QEYT4ET06r/cMzuX33n3p/KRz36l1jnXX3GxstgGyWgCQJ9NbZLRBDielawmKSWpejypXRq7feuWbN+6RVlsQwSaANBnejQBltdrNckLZ55Qq//yp16wJb/7nk8lOX5prLLYZiidBYA+6+zRnJXRBFhkudUkb7rz07XOuOjsM0yMPcVkNAGgz+zRBEZdr4mxTa0mOXtyIpdfdJaJsaeQQBMA+kzpLDDKepXFtstZm1pNkpgYeyoJNAGgz6w3AUZVuyx2acZy9765XHvrzpxIsrHOahJOHT2aANBnm9aPZ93C/wmaldEEhtDR+Sp33vNg3vaR+3PnPQ/m6HxVqyx2vmbNrNUkg0dGEwD6rJSSyYl1+fIjh2U0gaHTqzT26kvPq1UWO7FuLHNH5rt+zGqSwSXQBIABMDmxfiHQlNEEhsdypbE33X53rTNe8i1Pzi133JfEapK1ROksAAyAdp/m/rkjqXptFQdYQ5qaGLttZtpqkjVIRhMABkA70DwyX2Xu8Hw2bRjv840A6um1mmTHrr2NTYwdHytWk6wxAk0AGABLV5wINIG1oFf/5c++8Gty570P1j6nzsRYq0nWFqWzADAAOlecmDwLrAXt/sulWcvd++byM2/9h/zJ//lcrXNMjB1OMpoAMACmOjKasybPAgOkW2lskhX3X5oYO9wEmgAwAKY6MppWnACDoldp7Pc989xa/Zc/9Kwn59a//UwSE2NHjdJZABgAS3s0AfptudLY3373p2qd8cwLzjQxdkTJaALAAJiU0QT6oNfE2KZWk5w9OZHLLzrLxNgRJNAEgAEgowmcar3KYm+8aiabN21obDVJYmLsKBJoAsAAkNEETqV2WezSjOWefXO59tad+aozNtQ+q85qEkaPHk0AGABTmzqmzh6Q0QRWz3Jlse3HvvTQoVpnWU1CLzKaADAAZDSBpvXqv9yxa2+tsth2r2Y3VpNwPAJNABgAnYGmPZrASvXqv/ylK2fyD5/7Sq0zXnr5+bnljvuSWE3CiRNoAsAAmDIMCGhIr/7L3fvm8hNv3ln7nG0z07nswjOPCVinFwYGKY1lOQJNABgAG9eNZf14yeGjldJZoJZupbFJVryWpHNi7PhYsZqEkyLQBIABUErJ5MT67H34UPYflNEElterNPbqS8+r1X/53d9wTt7+0d1Jjl8WazUJJ8PUWQAYEO0+TRlNYDnt0tilAeXufXO56fa7a53x7TNPyM3XXGJiLKtGRhMABkS7T3P/3JFUVZVSlKbBqOo1MXa51SQn4uzJiVx+0VnKYlk1Ak0AGBDtjObR+SqPHDqa0zf6MQ2jqFdZ7I1XzWTzpg21SmN76ey/TJTFsnr8BAOAAbF0l6ZAE0ZPr4mxe/bN5dpbd+aCs06rfVbJ8fsvYbXo0QSAATFpxQmMtOXKYtuP3ffgI7XOuv6Ki/Vf0ld+VQoAA6IzozlrIBAMtW49mDt27a1VFrs0U7n0Y9ObJ7J965Zs37pF/yV9I9AEgAEhowmjoVcPZt1eyZc++/y86YOfTnL80lj9l/SL0lkAGBBTMpow9JZbTfKnO++vdcZ3fN05VpMw8GQ0AWBATMlowlBYrdUknRNjx8eK1SQMNIEmAAyIpVNngbVnudUk46XUXk1SZ2Ks1SQMMoEmAAwIPZqwtvVaTbJ7YTVJ3WTjy59zQW67a8+ioHR6IVhVFstaIdAEgAEhowlrQ7fS2CTHLYudr1kzu21mOjdcOaMsljVNoAkAA0KgCYOvV2ns1ZeeV6ssdmLdWOaOzHf92NIeTGWxrGWmzgLAgJjapHQWBtlyE2Nvuv3uWme85FuenJLHei7buvVgwlom0ASAAdGZ0Zw9IKMJg2SlE2Pbts1MW03CSFA6CwADYuO68WxYN5ZDR+YzK6MJfdFrNcmOXXtrT4ztxmoSRo1AEwAGyNTEunzpoUN6NKEPevVf/uwLvyZ33vNg7XOsJgGlswAwUNorTvRowqm1XP/lz7z1H/InOz9X65zrr7hYWSxERhMABkq7T/Ohg0dSVVVKUUoHTTrZ1STH0y6N3b51S7Zv3aIslpEn0ASAAdIONOer5OFDR3PGRj+qoSm9SmO/75nn1uq//KFnPTm3/u1nkhy/NFZZLKNO6SwADJCpCStOYDUsVxr72+/+VK0znnnBmSbGQk1+TQoAA2TpipNzNvfxMrDG9JoY29RqkrMnJ3L5RWeZGAs1CDQBYIBMymjCSelVFnvjVTPZvGlDY6tJEhNjoQ6BJgAMkM6MphUnUE+7LHZpxnLPvrlce+vOnHn6+q6v66bOahLg+PRoAsAA6cxozspownEtVxbbfmzvw/X+t2Q1CTRHRhMABoiMJvTWrQdzx669tcpi272a3VhNAs0TaALAAFk8dVagCW3dejCnpyby9U+aqvX6l15+fm65474kVpPAqSDQBIABMtU5dVbpLCRZpgdzdi57ZusN+dk2M53LLjzz2GB1YWCQ0lholkATAAaIqbOMqtVaTdI5MXZ8rFhNAqeIQBMABogeTUbRcqtJ9h04XHs1SZ2JsVaTwKkh0ASAASLQZNT0KovdvbCapK6XP+eC3HbXHmWxMCAEmgAwQJTOMqy6lcYmWVFZbKdtM9O54coZZbEwIASaADBANqwby8Z1Yzl4ZF5Gk6HRqzT26kvPq1UWe9qG8Rw4dLRrQLq0B1NZLAyGsX5fAABYbGpTK6sp0GQYtEtjlwaUu/fN5abb7651xtWXnpfksZ7Ltm49mMBgEGgCwIBp92nOHlA6y9q20omxbdtmpnPzNZdkevPEosenN0/k5msu0YMJA0jpLAAMmHaf5kOHjmR+vsqYTA0Drtdqkh279taeGNuN1SSwdgk0AWDATC1kNKuqFWxOdQwIgkHTq//yR7/1KXn/3Q/UPsdqEhguSmcBYMBYccJasVz/5a+84+N57yfqBZrXX3GxslgYMjKaADBgJjcuXXGyqX+XgazeapJ2aez2rVuyfesWZbEwRASaADBgpjbJaDI4epXGvvgbn1ir//K65z0lr3/fvUmOXxqrLBaGh9JZABgwkxNLM5rQH8uVxr7+/ffWOuNp50yZGAsjSEYTAAZMZ4/m7AEZTVZXr4mxTa0mOXtyIpdfdJaJsTBiBJoAMGBkNDlVepXF3njVTE7fsK6x1SSJibEwagSaADBgFmU09WiyStplsUszlrv3zeXaW3dm47r6HVZ1VpMAo0WPJgAMGOtNWG11ymIPHpmvdZbVJEA3A53RLKW8Msm/SPK0JAeSfDDJz1dV9Ym+XgwAVtGU0lka1K0Hc8euvbXKYjeMlxw62j0ctZoEWM5AB5pJnpfk95N8KK27/nqSvyqlzFRV9XBfbwYAq2RxoCmjycnr1oP5hKmNufCs02u9/ppnnZ9b7rgvidUkwIkZ6ECzqqoXdb5fSnlZki8m+eYk7+/HnQBgtS3u0ZTR5OT06sH8wuzBfGH2YK0zts1M57ILzzwmWJ1eGBikNBboZaADzS42L/y5t9cTSikbk2zseGhyVW8EAA07Q48mNa3WapLOibHjY8VqEuCErZlAs5QyluS3ktxRVdVdyzz1lUluPDW3AoDmrR8fy6b14zlw+KgeTXrqtZrkhiu/Nvc+8HDt1SR1JsZaTQKcqDUTaKbVq/n0JM89zvNek+Q3O96fTPK51boUAKyGyYl1C4GmjCbHWm41yfY3f7j2OS9/zgW57a49ymKBxq2JQLOU8ntJvivJt1VVtWzQWFXVwSSPNh6UoqwDgLVncmJdvrj/oEBzxHUrjU2yorLYTttmpnPDlTPKYoHGDXSgWVpR4u8m+Z4kz6+qalefrwQAp8TUptbk2YcOHsnR+cr/8R9BvUpjr770vFplsZMb1+Whg0e6BqRLezCVxQJNG+v3BY7j95Nck+QlSfaXUqYX3jb1+V4AsKomO1acPHRQVnPUtEtjlwaUu/fN5abb7651xvc+89wkj/VctnXrwQRo2qAHmtelNWn2vUl2d7x9fx/vBACrbtGKkwMGAo2SlU6Mbds2M52br7kk05snFj0+vXkiN19ziR5MYFUNdOlsVVV+zQbASJqy4mTo9VpNsmPX3toTY7uxmgQYBAMdaALAqOosnbXiZPh067+cnprI9zzjSfnA3Q/UPsdqEmBQDXrpLACMpMmNMprDqlf/5Z7Zudz8vnty1+dna51z/RUXK4sFBpaMJgAMoPbU2STZf1BGcy1ayWqSsZLM93hSuzR2+9Yt2b51i7JYYCAJNAFgAE3q0VzTeq0m2fq0r67Vf/nTW5+a3353a7rs8UpjlcUCg0jpLAAMoM4eTVNn15blVpP80d99ttYZF3716SbGAmuajCYADCAZzcHWa2JsU6tJzp6cyOUXnWViLLBmCTQBYAAt2qMp0Bwovcpib7xqJkkaW02SmBgLrF0CTQAYQFPWmwykdlns0ozl7n1zufbWnTmRXGOd1SQAa5UeTQAYQEpnB0+dsti6JbNWkwDDTkYTAAbQGYv2aMponkq9+i937Npbqyx2Yv1Y5g7Pd/2Y1STAqBBoAsAAWjc+ltM3jOfhQ0dlNE+hXv2XP/KtF+b9n3ig1hkvuezJueWO+5JYTQKMLoEmAAyoyYn1efjQ0czKaJ4Sy/Vf/vt3/GPtc7bNTOeyC888JmCdXhgYpDQWGAUCTQAYUJMT67JnVo9m07qVxiZZ8VqSzomx42PFahJgpAk0AWBAtQcCPXLoaI4cnc+6cTP8VqpXaeyLv/GJtfovr3veRXn9++5JcvyyWKtJgFHmJxYADKjJjhUnDx2U1Vypdmns0oBy9765vP7999Y642nnTObmay4xMRbgOGQ0AWBATW3q3KV5JI87bUMfb7M29JoYW2c1SR1nT07k8ovOUhYLcBwCTQAYUJ27NA0EOr5eZbE3XjWT0zesq1Ua20tn/2WiLBbgeASaADCgOgNNAzmr+rIAAA5RSURBVIGWt9zE2Gtv3ZmN6+p3C5Ucv/8SgOXp0QSAATXV0aM5e0BGM2mVxt55z4N520fuz533PJij81WtstiDR+ZrnX/9FRfrvwRogIwmAAwoGc3FepXGXn3pebXKYjeMj+XQ0e4BZ7s0dvvWLdm+dYv+S4AVEmgCwIBaHGiOdkazV2nsnn1zuen2u2udcc2znpxb7rgvyfFLY/VfAqyM0lkAGFCdpbOjnNFcrjT2RKbIbpuZtpoE4BSR0QSAAdW5R3P/COzR7LWaZMeuvY1NjB0fK1aTAJwCAk0AGFCjVDrbq//yhiu/Nh/+9Jdrn1NnYqzVJACrT6AJAANq0R7NA8Ob0VxuNcn2N3+49jnXX3Fx3vKhzywKVqcX9mgqiwU4tQSaADCgOktnZ4cgo9mtNDbJcVeTHI+JsQCDR6AJAAPqjI2P/Zj+zIOP5M57HlyzgdNKV5P8i2c8Kf/rw/cnMTEWYC0wdRYABtS7Pr7n0UDq03sfyQ/8wd/mua99T9551+6+3utEtUtjlwaUu09gNcnzvuarTYwFWENkNAFgAC23N/K6W3cOXHDVa2LscqtJTsTZkxO5/KKzTIwFWCMEmgAwYI63N7Kk1de4bWZ6IIKsXmWxN141k82bNjS2miQxMRZgrRBoAvD/t3evMZqWdxnAr/8uhy10d0sLy0EWQSjqCpZwWFuwghzarT2SNIJoIzZRQYjSqkmhKi2JoUmVtgqpCTEUA5rWL2BjCJVatVIEFgoUi1grFCincujMtrAcdm8/vDO4TGdmB/aB53lnf79ksrzPvu/M9eEOm2vuEwOztXsjW0bLTm+654neS9d8M69nXHFrVu/2mgV/r4VcTQLAeLBHEwAG5tENC5sBXOj7Xilbm3lNkvuffHpB3+tDJx5s/yXAImJGEwAGZtXyZVt/00t4Xxdm24O5tZnXaVVJm2OTpqtJABYnRRMABmbtAa/P3iuX5eGJjbPOFs7ct/hKm2sP5tr9F/bzTz96/3zu+nuTuJoEYHth6SwADMzSJZXz370myf+XsS21vHr7Fue7muTq2x9c0Pd425q9XE0CsJ0xowkAA7TukL3z2V87/EdmEpPkx9+wS962Zq/OftYrdTXJljOvS5eUq0kAtiOKJgAM1LpD9n6hnD008XT+/Et357vf35jvPP5UvrD+/py6dr9t/hnzXU0yfbrtQizkxFhXkwBsP6rNtTt/kaiqFUkmJiYmsmLFir7jAMDL9rVvP5bTLr0xSbLbLjvmK39wXF63y04v+/vNdTXJtJnlcS4fPGb/XHPnw7OWVctiAcbb5ORkVq5cmSQrW2uTC/2cGU0AGBNHH7h73vOmffIPtz+YJ596Lp+89u786cmHbvVzsy2NTbLVZbEL/VX0SWv2ykffucayWABeoGgCwBj56Dt/Ol++65H88NlN+dub7sspR63Oz+77ujnfP9fS2FOPWr2gZbHLdlySjc9tnvXvZu7BtCwWgGlOnQWAMbLnimU558SDk4zupvzjq+7M5s2zzz3Od2Lsp6771oJ+3mlr90vlR0+/nW0PJgBMUzQBYMycfsz+eeOq1yZJbn9gIhdec1euvu27ueHbj2fTVOnc1hNjp53kahIAXgZLZwFgzOy4dEkueO8h+ZVL/yNJculX73nh76YP4Vm249IFnxg7G1eTALAtFE0AGEMTTz876/OHJjbmjCtuTb2EDuhqEgC6ZuksAIyZ6WWx81no7WUfOvFgy2IB6JwZTQAYMzfd88SClsXuvMOSPPP8/CfGnn38QTn7+IMsiwWgU4omAIyZRzcsbO/lr/7cfrns+nuTbH1prGWxAHTJ0lkAGDOrli/b+pvixFgA+mNGEwDGzNoDXp+9Vy7LwxMbZ72+xImxAPRN0QSAMbN0SeX8d6/JmVfc6sRYAAbJ0lkAGEPrDtnbslgABqvaQs8/H1NVtSLJxMTERFasWNF3HADo1KbNzbJYAF4xk5OTWblyZZKsbK1NLvRzls4CwBizLBaAIbJ0FgAAgE4pmgAAAHRK0QQAAKBTiiYAAACdUjQBAADolKIJAABApxRNAAAAOqVoAgAA0ClFEwAAgE4pmgAAAHRK0QQAAKBTiiYAAACdUjQBAADolKIJAABApxRNAAAAOqVoAgAA0ClFEwAAgE4pmgAAAHRK0QQAAKBTiiYAAACdUjQBAADolKIJAABApxRNAAAAOqVoAgAA0ClFEwAAgE4pmgAAAHRK0QQAAKBTiiYAAACdUjQBAADolKIJAABApxRNAAAAOqVoAgAA0ClFEwAAgE4pmgAAAHRK0QQAAKBTiiYAAACdUjQBAADolKIJAABApxRNAAAAOqVoAgAA0ClFEwAAgE6NRdGsqrOq6t6q2lhVN1bV2r4zAQAAMLvBF82qOiXJRUk+nuTwJLcnubaqVvUaDAAAgFkNvmgm+XCSS1trl7XWvpnkjCRPJflgv7EAAACYzaCLZlXtlOSIJNdNP2utbZ56/Za+cgEAADC3HfoOsBW7J1ma5JEZzx9J8lOzfaCqdk6y8xaPlifJ5OTkK5EPAABg0Xq5PWroRfPlODfJ+TMfrl69uocoAAAAi8LyJAtunUMvmo8l2ZRkzxnP90zy8ByfuTCjw4OmLU/yQJJ9k2zoOiBjz/hgLsYG8zE+mI/xwVyMDeYz5PGxPMmDL+UDgy6arbVnq+qWJCckuSpJqmrJ1OuL5/jMM0memX5dVdP/uaG1Zv0sL2J8MBdjg/kYH8zH+GAuxgbzGfj4eMl5Bl00p1yU5PKqWp/kpiTnJNk1yWW9pgIAAGBWgy+arbXPV9UeSS5IsleS25Ksa63NPCAIAACAARh80UyS1trFmWOp7AI8k+Tj2WI5LWzB+GAuxgbzMT6Yj/HBXIwN5rOoxke11vrOAAAAwCKypO8AAAAALC6KJgAAAJ1SNAEAAOiUogkAAECnFnXRrKqzqureqtpYVTdW1dq+MzEMVfULVfXFqnqwqlpVva/vTAxDVZ1bVTdX1YaqerSqrqqqn+w7F8NQVWdW1R1VNTn1dUNVvaPvXAxPVX1k6t+XT/edhf5V1cemxsOWX//Vdy6Go6p+rKquqKrHq+rpqvpGVR3Zd65tsWiLZlWdkuSijI4IPjzJ7UmurapVvQZjKHbNaEyc1XcQBufYJJckeXOSk5LsmORLVbVrr6kYigeSfCTJEUmOTPLPSa6uqp/pNRWDUlVHJfntJHf0nYVB+c8ke2/x9fP9xmEoqmq3JNcneS7JO5KsSfL7SZ7sM9e2WrTXm1TVjUlubq2dPfV6SZL7k/xla+0TvYZjUKqqJTm5tXZV31kYnqraI8mjSY5trf1b33kYnqp6Iskfttb+uu8s9K+qXpvk1iS/k+SPktzWWjun31T0rao+luR9rbXD+s7C8FTVJ5Ic01p7a99ZurQoZzSraqeMftt83fSz1trmqddv6SsXMJZWTv35RK8pGJyqWlpVp2a0QuKGvvMwGJck+cfW2nVbfSfbmzdObdn536q6sqr26zsQg/GeJOur6u+ntu18vap+s+9Q22pRFs0kuydZmuSRGc8fSbLXqx8HGEdTKyE+neT61tqdfedhGKrq0Kr6QZJnkvxVRisivtlzLAZg6hcPhyc5t+8sDM6NSU5Psi7JmUkOSPLVqlreZygG4ycyGhffSvL2JJ9N8hdV9eu9ptpGO/QdAGDALklySOyj4cXuTnJYRrPd709yeVUdq2xu36pqdZLPJDmptbax7zwMS2vtmi1e3jG1xes7SX45iWX3LEmyvrV23tTrr1fVIUnOSHJ5f7G2zWKd0XwsyaYke854vmeSh1/9OMC4qaqLk7wryS+21h7oOw/D0Vp7trX2P621W1pr52Z0sNjv9Z2L3h2RZFWSW6vq+ap6PqPDxX536vXSfuMxJK217yf57yQH9Z2FQXgoycxfVt6VZKyXVy/KotlaezbJLUlOmH42tQTuhNhHA8yjRi5OcnKS41tr9/SdicFbkmTnvkPQuy8nOTSj2e7pr/VJrkxyWGttU4/ZGJipQ6MOzKhgwPVJZl6ldnBGs95jazEvnb0oo+VM65PclOScjA5suKzXVAzC1P/gt/wt4gFVdViSJ1pr9/UUi2G4JMlpSd6bZENVTe/rnmitPd1fLIagqi5Mck2S+5Isz2isHJfRnhq2Y621DUletJe7qn6Y5HF7vKmqP0vyxYyKwz4ZXb+3Kcnf9ZmLwfhUkq9V1XlJvpBkbZLfmvoaW4u2aLbWPj91LcEFGR0AdFuSda21mQcEsX06MslXtnh90dSfl2e0WZ/t15lTf/7LjOe/keRzr2oShmhVkr/J6A68iYzuSXx7a+2fek0FDN2+GZXKNyT5XpJ/T/Lm1tr3ek3FILTWbq6qk5NcmORPktyT5JzW2pX9Jts2i/YeTQAAAPqxKPdoAgAA0B9FEwAAgE4pmgAAAHRK0QQAAKBTiiYAAACdUjQBAADolKIJAABApxRNAAAAOqVoAgAA0ClFEwAAgE4pmgDwKquqParq4ao6b4tnR1fVs1V1Qp/ZAKAL1VrrOwMAbHeq6peSXJXk6CR3J7ktydWttQ/3GgwAOqBoAkBPquqSJCcmWZ/k0CRHtdae6TcVAGw7RRMAelJVr0lyZ5LVSY5orX2j50gA0Al7NAGgPwcm2Sejf4/37zcKAHTHjCYA9KCqdkpyU0Z7M+9Ock6SQ1trj/YaDAA6oGgCQA+q6pNJ3p/kTUl+kORfk0y01t7VazAA6IClswDwKquq4zKawfxAa22ytbY5yQeSvLWqzuw1HAB0wIwmAAAAnTKjCQAAQKcUTQAAADqlaAIAANApRRMAAIBOKZoAAAB0StEEAACgU4omAAAAnVI0AQAA6JSiCQAAQKcUTQAAADqlaAIAANApRRMAAIBO/R/qfXl+7cvT5wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(11, 7), dpi= 100)\n", "plt.plot(x, u, marker='o', lw=2)\n", "plt.xlim([0, 2 * np.pi])\n", "plt.ylim([0, 10]);\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.title('Burgers Equation at t=0');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This new function is known as a `sawtooth` function. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Periodic boundary conditions\n", "\n", "The biggest difference between this step and the previous ones is the use of periodic boundary conditions. If you have experimented with steps 1-2 you would have seen that eventually the wave moves out of the picture to the right and does not show up in the plot. \n", "\n", "With periodic BC, what happens now is that when the wave hits the end of the frame it wraps around and starts from the beginning again.\n", "\n", "Now we will apply the discretization as outlined above and check out the final results." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "for n in range(nt): #Runs however many timesteps you set earlier\n", " un = u.copy() #copy the u array to not overwrite values\n", " for i in range(1,grid_points-1):\n", " u[i] = un[i] - un[i] * dt/dx * (un[i]-un[i-1]) + nu * (dt/dx**2) * (un[i+1]- 2*un[i] + un[i-1]) \n", " \n", " u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu*(dt / dx**2) *(un[1] - 2* un[0] + un[-2])\n", " u[-1] = u[0]\n", "\n", "u_anal = np.asarray([ufunc(nt* dt , xi, nu) for xi in x])" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5oAAAJmCAYAAAA90v9dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3Xl8nWWd9/HPddI2bWmTtKUbpSttU0oBwbbaugAKFoGKDjIMwgg6OiMqPo7jzOioA4iMyuLgBj76FBiWAcWFTRaxCsjajVagbVq6A+lCl3SjW871/HGfpNmXcqcny+f9ep1XOOfefuck1Xxz3b/rCjFGJEmSJElKSybfBUiSJEmSOheDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEkdVAghhhCuzHcdkiTVZdCUpE4uhHBpLpDUfGwMIfw5hPDhfNeXTw18LjUfP8t3fQAhhLM6apgMIUwPIVwZQihp4f6fCCF8OeUavhFCeCCEsKG5YB5CGBZC+FUIYVsIYXsI4f4Qwpg065GkrqJbvguQJB02/wmsAgIwGLgUeDiEMDPG+FA+C8uzx4HbG3h92eEupBFnAV8ArmxgWy/gwGGtpnWmA1cAtwHbWrD/J4BJwI0p1vAdYD3wIjCjsZ1CCH2APwPFwH8B+4F/Bp4MIbwjxrg5xZokqdMzaEpS1/FIjHFe1ZMQwixgA3AhkErQDCEcEWPclca5WnCtDNAjxrjnbZ5qWYzxzjRqOtxSeO9dwegY4+oQwpHApib2+zwwDpgaY5wLEEJ4BHgZ+BfgP9q8UknqRLx1VpK6rm3AW9QYEQshnJq7vfDUmjuGEEblXr+0xmu3hRB2hhCOCSE8HELYAdxVY/sXQggrQwhvhRDmhBDeF0J4IoTwRJ1zF4YQrgohvBpC2BtCWBdCuDaEUFhnvxhC+EkI4aIQwivAXuDM3La/CyHMDyHsyN3y+FII4f+k9UHlrvGPIYQVTb2fGrcpj6pzbL3PNXf8vSGEtTXe93+HEHrV2Oc2ktHMWrf51vlMrqxzrZNCCI/kPoedIYTZIYR319mnqs73hBB+EELYFELYFUL4XQhhYAs+ixNy3/+VIYQ9IYT1IYRbQggDauxzJXBd7umqGvWPauCU5D7Hs4GRNfZd3VwtzYkxtvQcHwfmVoXM3LFLgdnA377dOiSpq3FEU5K6juLcqE4ABgGXA32AtzOa1w14DHga+CqwGyCEcBnwE+AvwH8Do4D7gK3Aa1UH50YlHwDeC/wcWAIcT3LL4njgo3Wu9wGSX/p/ArwJrA4hnAHcTRII/j2337HAe4AftuA99Mx9LnVtjzHuy9X5D8D/BZ4lua1zTK7uLcC6FlyjIecDvYGbgc3AVJLvydG5beSueRRwBvD3zZ0whHAcyWe+HbiW5PbPfwKeCCGcEmN8oc4hPyb5nlxF8j36Mslne0EzlzqD5DO4leS21OOAfwSOCyG8O8YYgd+SfA8vJPl+vpk7trFRxWtIbls9Orc/wM4a762h71FDdsQY97Zw36pzZ4ATgFsa2DwH+FAIoW+McUdrzitJXZlBU5K6jj/Web4X+HSM8fG3cc5C4N4Y49erXggh9ACuBuYCH4gxHsi9/leSXr3Xahz/CeB04JQY49M1zvEy8LMQwvQY47M19i8Fjo8xLq6x740kwWpGjLHyEN7DP+QedV0I3BNC6E7Ss7cQOK1G+FxMEo4PNWj+e4zxrRrPfx5CeBX4rxDCiBjj2hjjcyGEZcAZLby99ztAd+C9McaVuTpvB8pIgucpdfbfDHwoFwyrAteXQgjFMcaKJq5zU4zxhpovhBCeJwn87wX+EmP8awhhAcnneF9zI4sxxsdDCK8D/Rp5r03d9lrTp0h+zlqjP8nPcnkD26peO4rkc5QktYBBU5K6ji9wcIKbwcDFwP8LIeyIMf72bZz35jrPJwMDgK9Xhcycu0hGN2s6n2QUc2mdEas/5b6eRjKKWOXJmiEzZxtwBMko26OtL5/7SUbx6nop93UyyQjwf1aFzJzbOHhraKvVDJkhhCNIJvZ5lmTE+SRgbWvOF0IoAD5EEupW1rhOeQjhf4HPhhCKYozbaxz286qQmfMXktHEkcBfW1h7T5KR8edzL52cO0/azmjhfq8cwrmrblduaCR0T519JEktYNCUpK5jTp3JgO4mmYnzJyGEh+qEqJY6QO0RSkhCCsCrNV+MMR5ooOduHMltro2NVg2q83xVA/vcRHI77SO5EbE/AL+KMbY0dL4WY6w72ltT1ftZXvPFGOP+EMLKBvZvkRDCCODbwEeAfnU2Fx/CKQeS3Irb0KjbEpJ5GYZTO4jVDbNbc1/r1lNLCKE/yWyyf0f979Gh1N6sZr5Hb1dVcC5sYFvPOvtIklrAoClJXVSMMRtC+DPwf0gC3ytAbGT3gkZe3xtjzL6NMjIkI4dfaWR73dtS6/2yH2PcGEJ4B8nSFR/OPT4VQrg9xnjJ26jtULTo88uNPj5Ocsvm94GlwC5gGMlI6eGarK+xW41DM8f9imTpkutIbineSVLzo7RR7SGEIS3ctaLOLcktsYVkNHNoA9uqXnujleeUpC7NoClJXVvV/w/0yX2tGtEqqbPfSFpuTe7rWJJ1CQEIIXQjmXCm5i2ZK4ATgdl1buFsldxo7IPAg7k+w5uAfwohXB1jfLXpo5tV9X7GcfCWXnK9m6OBRTX2benndzzJRDmXxBir1/DMTWxUV0s/l00kkzGVNrBtApDl0PtJq4UQ+gEfBK6IMX67xuvjGti9td/TpvZvqH+yIa3u0cz90eUlktuk63oXsNKJgCSpdVzeRJK6qFxQ+hCwj+TWSkhCVSXw/jq7f74Vp55HMsnMZ3PhsspF1L8l81cko3ifbaC+XrnexSbVXFIDktDAwTDb0K2QrTWPJMR9LjfRUZVLqR8oV+S+Vn9+udHLf6yzX9VIYqixXyAZXa5rV2573WvVkpsI6Q/AuTWXEAkhDCaZdOnpOv2Zh6pe7TlfbmDfqjVVm6y9zv6N3Xp7Rgsfj7XwWnX9GpgSQqgOmyGEUpKZju89xHNKUpfliKYkdR0fDiFMyP33IJLwMQ74XlUAiTFWhBDuBS7Prde4AjiH+n14jYox7sutofhj4E8hhF+RjGRemjtfzVGrO0j6K38WQjgNeIbkNtMJuddnkAS9pvy/XM/gn0j6RUeSLBOykIMBuinjQwgXN/D6hhjj47lezG+SLDXypxDCL0lGMj8F1OrRjDG+kpt99bu5mraQ9DHW/f/bpSSfxfUhhGEks+aeR8O9kfNzX38UQngMqIwx3tPIe/kmSdh6OoRwE0kP7T+RBO5/a/wjaLkY4/YQwlPAv+X+WPE6yR8sRjdR+zUhhHtIllt5MMa4q4F9q/a/IITwA5JZi3fGGB/MXfeQejRDCH9P8jPRO/fS+3PfT4A7YoxVI9Y3kfzB4/chhOtztX4F2ADUmmFXktQCMUYfPnz48NGJHyQBL9Z5vEUyEdDngFBn/yNJRnd2kQSln5GskxiBS2vsdxtJEGjsupcDq0lm7XyBpKdvHvBInf26k4Sgl3P7bsnt959AUY39IvCTBq5zHsko1gaSPrs1uZqHtOCzqfu51Hw8UWffy0iC5R6SEPQ+4IkG9htD0n+5h2SNyWtIlnCJwKk19js2t98OkhHTn5Os5Vj3cy4AfgRsJLn9Ndap/8o61z+JpFdyR+57+CdgWiM/E5PrvH5q3Tob+dyGkayTuZVk1t9fkfQyNlTPN0n+AFCZ2z6qifMeQTI78dbcvqtT+Pl/oonv8al19j2aZPSyIvf5PQiMzfe/YR8+fPjoiI8Q4yG3xEiS1GK53slNwG9jjPVule2IQghPAMQYT81vJZIktS/2aEqSUhdC6JnrOazpkySzrD5x+CuSJEmHkz2akqS28G7gv3P9npuBk4F/ILk91olVJEnq5AyakqS2sJpkKY0vkYxibgFuB74Wk6VIJElSJ5bXHs0QwvuBfwXeSTKJwMdijPfV2B6Aq0hmgSshmY3wshjj8jyUK0mSJElqgXz3aB5BstD1FxrZ/m8kfw3/HMmCybuAx0IIPQ9PeZIkSZKk1mo3s87m1murHtHMjWa+AdwQY7w+91oxyfT1l8bG1xCTJEmSJOVRe+7RHA0MAaoXaI7JQuIvANOABoNmCKGQZGHqmqr6gyRJkiRJrdMXeCO2YpSyPQfNIbmvG+q8vqHGtoZ8HbiiTSqSJEmSpK7paOD1lu7cnoPmofou8IMaz/sCr61bt46ioqI8lSRJkiRJHc/27dsZPnw4wI7WHNeeg+b63NfBQHmN1wcDCxs7KMa4F9hb9bxqvfCioiKDpiRJkiQdBvmedbYpq0jC5gerXgghFJHMPvtcvoqSJEmSJDUtryOaIYQ+wNgaL40OIbwD2BJjXBtCuBH4ZghhOUnwvJpkJtr76p9NkiRJktQe5PvW2cnAn2s8r+qt/B/gUuBakrU2fw6UAE8DZ8YY9xzGGiVJkiRJrdBu1tFsK7nbbSsqKirs0ZQkSVKHV1lZyf79+/NdhjqRHj16kMk03FW5fft2iouLAYpjjNtbes58j2hKkiRJaoEYI+vXr2fbtm35LkWdTCaTYfTo0fTo0SO1cxo0JUmSpA6gKmQOGjSI3r17V6+uIL0d2WyWN954g/LyckaMGJHaz5VBU5IkSWrnKisrq0PmgAED8l2OOpmBAwfyxhtvcODAAbp3757KOdvz8iaSJEmSoLons3fv3nmuRJ1R1S2zlZWVqZ3ToClJkiR1EN4uq7bQFj9XBk1JkiRJUqoMmpIkSZJ0mJx66ql8+ctfPuzXDSFw3333HbbrGTQlSZKkLqIyG3luxWbuX/g6z63YTGU2Hpbrrl+/nssvv5wxY8ZQWFjI8OHDmTlzJrNnzz4s1387brvtNkpKSlp93BNPPEEIod5yNL/97W+5+uqr0yqv3XLWWUmSJKkLePTlcq56cDHlFXuqXxta3JMrZk7kzElD2+y6q1ev5j3veQ8lJSVcd911HH/88ezfv5/HHnuML3zhCyxdurTNrt0e9e/fP98lHBaOaEqSJEmd3KMvl3PZnQtqhUyA9RV7uOzOBTz6cnmbXfvzn/88IQTmzJnDeeedx/jx4znuuOP4yle+wvPPPw/A2rVrOffcc+nTpw9FRUX87d/+LRs2bKg+x5VXXsk73vEObrnlFkaMGEGfPn34/Oc/T2VlJddeey1Dhgxh0KBBXHPNNbWuHULg5ptv5sMf/jC9evVizJgx/PrXv67e3tCo48KFCwkhsHr1ap544gk+9alPUVFRQQiBEAJXXnklAHfccQeTJ0+mb9++DBkyhE984hNs3LgRSML1aaedBkC/fv0IIXDppZcC9W+d3bp1K5/85Cfp168fvXv35sMf/jDLly+v3l41ovrYY49x7LHH0qdPH84880zKyw9+z+bOncsZZ5zBkUceSXFxMaeccgoLFix4O9+2t82gKUmSJHVildnIVQ8upqGbZKteu+rBxW1yG+2WLVt49NFH+cIXvsARRxxRb3tJSQnZbJZzzz2XLVu28OSTT/L444+zcuVKLrjgglr7rlixgkceeYRHH32Uu+++m1mzZnH22Wfz2muv8eSTT/L973+fb37zm7zwwgu1jvvWt77Feeedx6JFi7jooov4u7/7O5YsWdKi+qdPn86NN95IUVER5eXllJeX89WvfhVIlpy5+uqrWbRoEffddx+rV6+uDpPDhw/nN7/5DQBlZWWUl5fzwx/+sMFrXHrppcybN48HHniA5557jhgjZ511VvWSNgC7d+/m+uuv54477uCpp55i7dq11XUA7Nixg0suuYSnn36a559/nnHjxnHWWWexY8eOFr3PtuCts5IkSVIHNPPHT7Npx95m99t7oJKtu/c3uj0C5RV7mPydxynsVtDs+Qb2LeTBy9/bohpfffVVYoxMmDCh0X1mz57NSy+9xKpVqxg+fDgAt99+O8cddxxz585lypQpAGSzWW655Rb69u3LxIkTOe200ygrK+Phhx8mk8lQWlrK97//ff785z/zrne9q/r8559/Pp/5zGcAuPrqq3n88cf58Y9/zE033dRs/T169KC4uJgQAkOGDKm17dOf/nT1f48ZM4Yf/ehHTJkyhZ07d9KnT5/qW2QHDRrUaI/n8uXLeeCBB3jmmWeYPn06AHfddRfDhw/nvvvu4/zzzweSUPuzn/2MY445BoAvfvGLfPvb364+zwc+8IFa5/35z39OSUkJTz75JOecc06z77MtGDQlSZKkDmjTjr2s376n+R1bKAmjjQfSQxFj86OkS5YsYfjw4dUhE2DixImUlJSwZMmS6qA5atQo+vbtW73P4MGDKSgoIJPJ1Hqt6vbVKtOmTav3fOHChYf0fmqaP38+V155JYsWLWLr1q1ks1kguQ144sSJLTrHkiVL6NatW61gPGDAAEpLS2uNuvbu3bs6ZAIMHTq01vvcsGED3/zmN3niiSfYuHEjlZWV7N69m7Vr177dt3nIDJqSJElSBzSwb2GL9mtuRLNKv97dWzyi2VLjxo0jhJDKhD/du3ev9TyE0OBrVYGvJapCas1AXPOW1cbs2rWLGTNmMGPGDO666y4GDhzI2rVrmTFjBvv27Wvx9VuqofdZs+ZLLrmEzZs388Mf/pCRI0dSWFjItGnT2qSWljJoSpIkSR1QS29frcxG3vv9P7G+Yk+DfZoBGFLck6f//QMUZEKqNfbv358ZM2bw05/+lC996Uv1+jS3bdvGsccey7p161i3bl31qObixYvZtm1bi0cGm/L888/zyU9+stbzk046CYCBAwcCUF5eTr9+/QDqjXb26NGDysrKWq8tXbqUzZs3873vfa+65nnz5tU7Dqh3bE3HHnssBw4c4IUXXqi+dXbz5s2UlZW16r0/88wz3HTTTZx11lkArFu3jjfffLPFx7cFJwOSJEmSOrGCTOCKmUloqRsjq55fMXNi6iGzyk9/+lMqKyuZOnUqv/nNb1i+fDlLlizhRz/6EdOmTeP000/n+OOP56KLLmLBggXMmTOHT37yk5xyyilMnjz5bV//3nvv5ZZbbmHZsmVcccUVzJkzhy9+8YsAjB07luHDh3PllVeyfPlyfv/733PDDTfUOn7UqFHs3LmT2bNn8+abb7J7925GjBhBjx49+PGPf8zKlSt54IEH6q2NOXLkSEIIPPTQQ2zatImdO3fWq23cuHGce+65fPazn+Xpp59m0aJFXHzxxQwbNoxzzz23xe9x3Lhx3HHHHSxZsoQXXniBiy66iF69eh3Cp5Ueg6YkSZLUyZ05aSg3X3wyQ4p71np9SHFPbr745DZdR3PMmDEsWLCA0047jX/5l39h0qRJnHHGGcyePZubb76ZEAL3338//fr14/3vfz+nn346Y8aM4Ze//GUq17/qqqu45557OOGEE7j99tu5++67q0cLu3fvzt13383SpUs54YQT+P73v893vvOdWsdPnz6dz33uc1xwwQUMHDiQa6+9loEDB3Lbbbdx7733MnHiRL73ve9x/fXX1zpu2LBhXHXVVXzta19j8ODB1eG2rltvvZV3vvOdnHPOOUybNo0YIw8//HC922WbMmvWLLZu3crJJ5/M3//93/OlL32JQYMGtfKTSldoSYNuRxZCKAIqKioqKCoqync5kiRJUqvt2bOHVatWMXr0aHr27Nn8AY2ozEbmrNrCxh17GNS3J1NH92+zkcz2IITA7373Oz760Y/mu5R2ramfr+3bt1NcXAxQHGPc3tJz2qMpSZIkdREFmcC0Ywbkuwx1Ad46K0mSJElKlSOakiRJkjqlzt4m2J45oilJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJ6tBGjRrFjTfe+LbO8cQTTxBCYNu2banUtHr1akIILFy4MJXzdTQGTUmSJElt6rnnnqOgoICzzz4736UAcOqpp/LlL3+51mvTp0+nvLyc4uLiPFXVuRg0JUmSJLWpWbNmcfnll/PUU0/xxhtv5LucBvXo0YMhQ4YQQsh3KZ2CQVOSJElSm9m5cye//OUvueyyyzj77LO57bbbqrdV3a46e/ZsJk+eTO/evZk+fTplZWXV+6xYsYJzzz2XwYMH06dPH6ZMmcIf//jHRq/36U9/mnPOOafWa/v372fQoEHMmjWLSy+9lCeffJIf/vCHhBAIIbB69eoGb5195plnOPXUU+nduzf9+vVjxowZbN26FYBHH32U9773vZSUlDBgwADOOeccVqxYkdKn1vEZNCVJkqSOJkbYtys/jxhbVeqvfvUrJkyYQGlpKRdffDG33HILsc45vvGNb3DDDTcwb948unXrxqc//enqbTt37uSss85i9uzZvPjii5x55pnMnDmTtWvXNni9z3zmMzz66KOUl5dXv/bQQw+xe/duLrjgAn74wx8ybdo0PvvZz1JeXk55eTnDhw+vd56FCxfywQ9+kIkTJ/Lcc8/x9NNPM3PmTCorKwHYtWsXX/nKV5g3bx6zZ88mk8nwsY99jGw226rPp7Pqlu8CJEmSJLXS/t3wX0fl59r/8Qb0OKLFu8+aNYuLL74YgDPPPJOKigqefPJJTj311Op9rrnmGk455RQAvva1r3H22WezZ88eevbsyYknnsiJJ55Yve/VV1/N7373Ox544AG++MUv1rve9OnTKS0t5Y477uDf/u3fALj11ls5//zz6dOnD5DcJtu7d2+GDBnSaN3XXnstkydP5qabbqp+7bjjjqv+7/POO6/W/rfccgsDBw5k8eLFTJo0qaUfT6fliKYkSZKkNlFWVsacOXO48MILAejWrRsXXHABs2bNqrXfCSecUP3fQ4cOBWDjxo1AMqL51a9+lWOPPZaSkhL69OnDkiVLGh3RhGRU89ZbbwVgw4YNPPLII7VGSVuiakSzMcuXL+fCCy9kzJgxFBUVMWrUKIAm6+pKHNGUJEmSOpruvZORxXxdu4VmzZrFgQMHOOqog6OvMUYKCwv5yU9+cvCU3btX/3fVZDxVt6B+9atf5fHHH+f6669n7Nix9OrVi49//OPs27ev0et+8pOf5Gtf+xrPPfcczz77LKNHj+Z973tfi+sG6NWrV5PbZ86cyciRI/nFL37BUUcdRTabZdKkSU3W1ZUYNCVJkqSOJoRW3b6aDwcOHOD222/nhhtu4EMf+lCtbR/96Ee5++67mTBhQrPneeaZZ7j00kv52Mc+BiQjnKtXr27ymAEDBvDRj36UW2+9leeee45PfepTtbb36NGjuteyMSeccAKzZ8/mqquuqrdt8+bNlJWV8Ytf/KI6wD799NPNvpeuxKApSZIkKXUPPfQQW7du5R/+4R/qrU153nnnMWvWLK677rpmzzNu3Dh++9vfMnPmTEIIfOtb32rRhDuf+cxnOOecc6isrOSSSy6ptW3UqFG88MILrF69mj59+tC/f/96x3/961/n+OOP5/Of/zyf+9zn6NGjB3/+8585//zz6d+/PwMGDODnP/85Q4cOZe3atXzta19rtqauxB5NSZIkSambNWsWp59+er2QCUnQnDdvHn/961+bPc8PfvAD+vXrx/Tp05k5cyYzZszg5JNPbva4008/naFDhzJjxoxat+5CcjtuQUEBEydOZODAgQ32VY4fP54//OEPLFq0iKlTpzJt2jTuv/9+unXrRiaT4Z577mH+/PlMmjSJf/7nf25RaO5KQt2phTubEEIRUFFRUUFRUVG+y5EkSZJabc+ePaxatYrRo0fTs2fPfJfTIezcuZNhw4Zx66238jd/8zf5Lqdda+rna/v27VV/LCiOMW5v6Tm9dVaSJElSp5HNZnnzzTe54YYbKCkp4SMf+Ui+S+qSDJqSJEmSOo21a9cyevRojj76aG677Ta6dTPy5IOfuiRJkqROY9SoUXT29sCOwMmAJEmSJEmpMmhKkiRJHYQjdWoLbfFzZdCUJEmS2rnu3bsDsHv37jxXos5o3759ABQUFKR2Tns0JUmSpHauoKCAkpISNm7cCEDv3r0JIeS5KnUG2WyWTZs20bt371QnTjJoSpIkSR3AkCFDAKrDppSWTCbDiBEjUv3jhUFTkiRJ6gBCCAwdOpRBgwaxf//+fJejTqRHjx5kMul2VRo0JUmSpA6koKAg1V46qS04GZAkSZIkKVUGTUmSJElSqgyakiRJkqRUGTQlSZIkSakyaEqSJEmSUmXQlCRJkiSlyqApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJklJl0JQkSZIkpcqgKUmSJElKlUFTkiRJkpQqg6YkSZIkKVUGTUmSJElSqgyakiRJkqRUGTQlSZIkSakyaEqSJEmSUmXQlCRJkiSlyqApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJklJl0JQkSZIkpapdB80QQkEI4eoQwqoQwlshhBUhhG+FEEK+a5MkSZIkNaxbvgtoxr8DlwGXAK8Ak4FbgQrgR3msS5IkSZLUiPYeNKcD98cYf597vjqEcCEwNY81SZIkSZKa0K5vnQWeBT4YQhgPEEI4EXgv8EhjB4QQCkMIRVUPoO/hKVWSJEmSBO1/RPN7QBGwNIRQCRQA34gx3tXEMV8HrjgcxUmSJEmS6mvvI5p/C1wEfAI4maRX86shhEuaOOa7QHGNx9FtXaQkSZIk6aD2PqJ5HfC9GOM9uecvhRBGkoxa/k9DB8QY9wJ7q547Qa0kSZIkHV7tfUSzN5Ct81ol7b9uSZIkSeqy2vuI5oPAN0IIa0mWNzkJ+ApwS16rkiRJkiQ1qr0xJrEWAAAgAElEQVQHzcuBq4GbgEHAG8D/Bb6dz6IkSZIkSY1r10EzxrgD+HLuIUmSJEnqAOx1lCRJkiSlyqApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJklJl0JQkSZIkpcqgKUmSJElKlUFTkiRJkpQqg6YkSZIkKVUGTUmSJElSqgyakiRJkqRUGTQlSZIkSakyaEqSJEmSUmXQlCRJkiSlyqApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJklJl0JQkSZIkpcqgKUmSJElKlUFTkiRJkpQqg6YkSZIkKVUGTUmSJElSqgyakiRJkqRUGTQlSZIkSakyaEqSJEmSUmXQlCRJkiSlyqApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJklJl0JQkSZIkpcqgKUmSJElKlUFTkiRJkpQqg6YkSZIkKVUGTUmSJElSqgyakiRJkqRUGTQlSZIkSakyaEqSJEmSUmXQlCRJkiSlyqApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKqDJqSJEmSpFQZNCVJkiRJqTJoSpIkSZJSZdCUJEmSJKXKoClJkiRJSpVBU5IkSZKUKoOmJEmSJClVBk1JkiRJUqoMmpIkSZKkVBk0JUmSJEmpMmhKkiRJklJl0JQkSZIkpcqgKUmSJElKlUFTkiRJkpQqg6YkSZIkKVUGTUmSJElSqgyakiRJkqRUtfugGUIYFkK4M4SwOYTwVgjhpRDC5HzXJUmSJElqWLd8F9CUEEI/4Bngz8CHgU3AOGBrPuuSJEmSJDWuXQdN4N+BdTHGT9V4bVW+ipEkSZIkNa+93zr7EWBeCOHeEMLGEMKLIYTP5rsoSZIkSVLj2nvQHANcBiwHZgA3Az8KIVzS2AEhhMIQQlHVA+h7eEqVJEmSJEH7v3U2A8yLMf5H7vmLIYRJwOeA/2nkmK8DVxyO4iRJkiRJ9bX3Ec1yYHGd15YAI5o45rtAcY3H0W1TmiRJkiSpIe19RPMZoLTOa+OBNY0dEGPcC+yteh5CaJvKJEmSJEkNau8jmv8NvDuE8B8hhLEhhE8A/wj8NM91SZIkSZIa0a6DZoxxLvAx4ELgZeBbwJdjjHfltTBJkiRJUqPa+62zxBgfAh7Kdx2SJEmSpJZp1yOakiRJkqSOx6ApSZIkSUqVQVOSJEmSlCqDpiRJkiQpVQZNSZIkSVKq2v2ss5IkSZK6qLe2QcVrucc62LUJhk+FMadBpiDf1akJBk1JkiRJ+XNgL6x/CV6fD5vKkkBZFS73bm/4mKJhcOKFcNJF0H/M4a1XLRJijPmuoU2FEIqAioqKCoqKivJdjiRJktR1ZbOwZSW8Pi8Jlq/NS0Jmdn/jx/TqD8VHQ/FwKOwDyx6DPdsObh/5XjjpYpj4EehxRNu/hy5m+/btFBcXAxTHGBtJ/vUZNCVJkiS1jb07k1C59nlY90ISLvdU1N+v9wAY9k4YcjyUjDgYLIuGJeGypv17oOxhWHgXvDobyOWZHn1h0sdgymdh6Alt/ta6CoNmIwyakiRJ0mFS8Tqsex7WvpB8Xf8yxMra+xQUwtAT4ejJSbgc9k7oNwpCOITrvQaL7oYX74Ktq5LXuvWEf/oLDBz/tt+ODJqNMmhKkiRJbSCbhU1LYO1zyYjl2hegYm39/YqHw/B3wYh3J+Fy8CQo6J5uLTHCmmfhD9+AN16EaV+EGdeke40uyqDZCIOmJEmSlIL9b8HrCw4Gy3VzYG+d22BDJrn9dfi7YcS7kq/Fww5fjUsfhnsuhCMGwleWpB9ou6BDDZrOOitJkiSpvl2bk77KqmD5xov1J+3pfgQMnwIjpiWjlkdPqd9TeTiNOyMJmbs2wfLHYcJZ+aulizNoSpIkSV1djLB1de4W2FywfLOs/n59BiehcsS7k8fg46GgHUWKgu5wwgXw3E+SyYIMmnnTjn4qJEmSJB0WlQdgw8u1g+XO9fX3O7L0YKgcMe3QJ+05nE66OAmayx6FnZugz8B8V9QlGTQlSZKkzm7frmRpkapguW4u7NtRe59MdzjqpIOhcvi74IgB+an37Rh0bDKT7evz4a+/hOlfzHdFXZJBU5IkSepsdm7KLTOSC5bliyB7oPY+hUUHZ4Md8e4knHXvlZ960/aOi5KgufAumPaF9j8K2wkZNCVJkqSOLEbYsrLGbbDPweZX6+/X9ygYOe1gj+WgiZApOPz1Hg6TzoPH/gM2Lk4mMRp2cr4r6nIMmpIkSVJHUnkA1v+1dn/lro319xs0MQmUw9+dBMzi4V1nZK9XCUw4B17+dTKqadA87AyakiRJUnu2dye8Pi8JlGuehdfmwf5dtfcp6AFHnVyjv3Iq9O6fn3rbi5MuSoLmS/fCh66B7j3zXVGXYtCUJEmS2pMdG+r0V/4VYmXtfXoWJyOVVf2VR51skKpr9CnJKG7FOlj6EBz/8XxX1KUYNCVJkqTDrDIbmbNqCxu3v8WIWM6JcTGZdS8kwXLLynr7x+LhhOr1K6fBwAmQyeSh8g4kUwAnXghPXZvcPmvQPKwMmpIkSVKKqkPkjj0M6tuTqaP7U5DJ9UYe2Mdzz/yJuU89TOm+V3hvpowBofYyI5HA8jCS5/aPY362lLnZUigcxhWlEzlz0tCWXUeJd3wiCZor/gwVr0Hx0fmuqMsIMcZ819CmQghFQEVFRQVFRUX5LkeSJEkdXFMB79GXy7nqwcWUV+wBoA+7+WDftVx+zCbG7nmZynVzKajcU+t8e2J3FsVjmJct5ciJp3DNX/uynSNq7VMVH2+++GTOnDS03nUAhhb35IqZtcOogFvPhjVPwwe+Ce//13xX0+Fs376d4uJigOIY4/aWHmfQlCRJklqoqYAHcMWds5mcKWNypowpmTKODWsoCLV/394a+zAvW8rc7HjmZUt5OY5mH91bdP2hxT351tkT+cL/LqDub/F1w6hyFv4v3HcZ9BsNX3qx68y8mxKDZiMMmpIkSWqp5kYrL7vzYMALZDkmvMHUTBnvzJTxroJlHE39ZUbWZAcxLya3wM7NlrIyDiVy6P2V3TKBA9mGf4cPwJDinjz97x/wNtoq+3bB9eNh30649GEY9Z58V9ShHGrQtEdTkiRJounRyjMmDuG/HljESWEZU3IjlpMzy+gXdtY6R2UMLIkjmZstzY1alrKRfqnW2VjIBIhAecUe5qzawrRjBqR63Q6rxxFw3MfgxTuSSYEMmoeFQVOSJEldRmMjlnVHKwGK2MWEHS+y4p7bGHPEKh7fu5TCwv21zvdW7MGL2bHMjUmwfDE7lp30PrxvqgEbd+xpfqeu5KSLk6D5yn3w4WuhsE++K+r0DJqSJEnqEhobsfzW2cdy9e+XMITN1aOVUzJllIZ1ZKr6K/cCAd6MRdX9lfOzpbwcR3GgFb9S9z+iB1t37avXX1klEyBGGtwegH5HdGfLrv0NbK1tUN+ezkpb0/B3wYCxsPlVWHxfEjzVpuzRlCRJUqfQ2v7KceH1WsHy6PBmvXOuyg5OgmVuxHJlHMrBaXdarqp3smoiH6gdJqvO+I/vH83Pn1rV6PaffuIkrv79EtZX7Gk0rAJ8cMJAXnljO+u3761+rcvPSvuXG2D2t5N1SD/9aL6r6TCcDKgRBk1JkqTOr7n+ytO+9yiDdixmSmZZ9aywxWF3rXMciBleiaOYly1lXm5G2E2UVG/vU1jArr2VjY42FvfuTsXuZLSxoZDY0qVJWrL9sjvrh9XmdPlZabe/Af99HMQsXL4ABhyT74o6BINmIwyakiRJnUPr+it38s7McqZkynh/4XLGHVhOYThQ63y7YyELsmOrZ4RdmB3LLno1ev1/Pn08N/5xGdB4kARatL5lc7e1Nre9sTA6ZVQ/HlhU3uh76PKz0t75cXj1cXjfv8AH/zPf1XQIzjorSZKkTqvJ/sqHFnMUm6pvgZ2cWcaEzLqDB1cCATbF4urZYOdlx7M4jmxRf2VVOPviB8ZSOqRPvTqG1AmSZ0wc0mxvZEEmNDkrbHPbz5w0tMHrzFm1pcmg2eVnpZ10XhI01zyb70o6PYOmJEmS8q41/ZUZspSGdUzeWcaBe6/jN5kyjuq5pd45V2SHJsEyN2K5Jg6muf7KQMOjlVfMnEhBJjQa8GoGyeZCYloauk5LZ5vtsrPSHjk++bp1TX7r6AIMmpIkScqr5vorv3v/i0wNS6pHLE/OLKMovFXrHPtjAS/H0czNljI/O5552fFsprh6e3Gv7oS39jfaX1k1Uc/Vv296tBIOX5A8FIP69kx1v06nZETydUc5HNgL3QrzW08nZtCUJElSm2tNf2U/tnPcjvmsuWcWS7ov5/HsCnoUVtY6387YkwXZcdUjli9mx7KHxkPDp98zmhv/uKzJEcszJw1lxqTmb3ttz6aO7s/Q4p5Nzko7pDh5X13SEUdC996wfzdUvOaEQG3IoClJkqQ21Vx/5dFhA1NC1TIjyxiXef3gwREIsCGWMDc7gXnZ8czNTmBpHE4lBc1eu7X9le15tLIlCjKBK2ZO5LI7F9QL1VWOHdK3Q4XnVIWQjGpuWgrb1hg025BBU5IkSW9La/orC6hkQljLlJ1lZO/9PvdlyhhcuK3eOZdnhzE3m/RWzo2lvBYHcjj6KzuDMycN5eaLT64Xqqv8uWwTTy3bxPvHD8xDde1AycgkaNqn2aYMmpIkSTpkzfVXfu/++bw7s5jJIemvPCnzKn3r9FfuiwW8FMfU6q/cysFl6Uq6SH9lmhoK1UvKK/j2Q0sA+Oq9i3jsy++n3xE98lxpHlT1aW5bm986OjmDpiRJkprUmv7K/mznhB1zWXfP/2Np9+U8nl1J9x61+yu3x17Mz46vXmpkUTyGvTQeeD7VRfor01Y3VL9rdH+eWPYmTy3bxMYde/mP373ETRedTAhd7DPqNzL5us0RzbZk0JQkSVKjmuqv/PZDixkR1idrV+ZGLI/J1FjDMddfWR7750Jl0l9ZFoeTJdPstbtaf2Vby2QC1338BGbc+BTbdu/nkZfXc+/81xjer3fXCueOaB4WIcbG5qPqHEIIRUBFRUUFRUVFze4vSZLUlbS2v3JiWJMEy9xSIwNDRb1zLs0Oz4XKZMTydY7kUPsrb7745Oog2VStarlHXy7nc3cuAOp/7kMbuN2403ljIfz8FDhiEPzr8nxX0+5t376d4uJigOIY4/aWHueIpiRJUhfVfH/lPKZnXqkesTwp8ypHhL21zrE3dmNRPIZ5uRHL+dnxVNCnerv9le3PmZOGMm3MAJ5bubne92V9xR4uu3NBrYDf6VSNaO7aCPvfgu698ltPJ2XQlCRJ6sRa0185kG2cuGMZr9/zC5Z1X8Yfs6vo1iNb63wVsXcuVJYyNzuel+IY+ys7mMpsZOWbOxvclrvbmaseXMwZE4d0zu9Br35QWAR7tye3zw4szXdFnZJBU5IkqZNqrr9ydHgjuQU2t4bl6MyGgwfnEsdr8cjq2WDnZCewPA4j2l/Zoc1ZtYUN2/c2uj0C5RV7mLNqS+f8nlStpbnhZYNmGzqkoBlC+M+mtscYv31o5UiSJKmlWtNf2Z0DHBdWM3lnGQX3lvH7TBkDCnfUOl82BpbGEczNjq8etSyn+aDh+pUdy8Yd9dfWfDv7dUglI5OguXV1vivptA51RPNjdZ53B0YDB4AVgEFTkiSpDTXXX3nt/XN5X+bl3IjlMt6ReZVeYV+tc+yJ3VkYxzIvFywXZMexnSOqt9tf2TkN6tsz1f06JGeebXOHFDRjjCfVfS03u+ttwO/eZk2SJEmidf2Vg9nCyTuWUX7P/+XVbmU8HtdQ0KN2RNwa+1T3Vs7LlvJSHMP+Jn4dtL+yc5o6uj9Di3uyvmJPk39EmDq6/+Eu7fBxLc02l1qPZoxxewjhCuBB4I60zitJktQVNdVfefWDr3BMeO3gMiOhjBGZTbVPEGBNdhDzYjIb7JzsBFbGofZXioJM4IqZE7nszgX1/ohQpeq2507LEc02l/ZkQMW5hyRJkprQmv7KHuxnUljFlJ1ldL+3jIczy+hXWHvW0MoYWBxH5kYsk/7KjfRrtg77K7umMycN5eaLT673RwSAqz5yXOdd2qRKSW5Ec6sjmm3lUCcD+lLdl4ChwN8Dj7zdoiRJkjqz5vorr7t/DqdU9VdmynhHWEFh2F/rHG/FHryYHcvcWFrdX7mLg+sB2l+p5tT8I8Ldc9bywKI3AFi/vRNPAlSlakTzrS2wdwcU9s1vPZ1QiLGh//lp5qAQVtV5KQtsAv4EfDfGuKP+UfmR6x2tqKiooKioKN/lSJKkLqI1/ZVD2Vx9G+y7uy1jbFxLJtT+He3NWFSrv/KVOIoDTYwZ/PPp47nxj8uAhkcsb774ZM6cNLTJkVV1HZt27GXad2dzIBsZ1LeQZ7/2AboVNH+bdYf2/VHw1la47FkYfFy+q2m3tm/fTnFxMUBxjHF7S4871MmARh/KcZIkSV1Bc/2V48K6g/2VmTKODm/WPkGAVdnByRIjMbkVdmUcysGY2Dj7K3UoBvYt5AMTBvGHxRvYuGMvTy7bxAePHZzvstpWyYgkaG5dY9BsA2n3aEqSJHV6remvLGQfx4eVTNm5jMJ7y3g0U0Zx4e5a5zsQM7wSR1X3V87PjmcTJc3WYX+l0nTBlOH8YfEGAH45d10XCJojoXyREwK1EYOmJElSKzTXX3n9fc9zWubl6hHLE8JKCsOBWufYFQuT/srsBObF8byYHcduDq5ZaH+l8uGU8QMZ1LeQjTv28qelG9m0Yy8D+xbmu6y24xInbcqgKUmSVENrRishMow3mbyjjDfv+Skruy3jj6yDHrXPuSkWV88EOzdbypI4osn+StevVD50K8jw8XcezU1PrOBANvLbBa/xT6cck++y2k7VzLOOaLYJg6YkSVJOc6OV377/JSaENdW9lZMzZRwVttQ7z4rs0OQW2JisX7kmDsb+SnUEfzt5ODc9sQKAX85bxz++fwwhdNI/XrjESZsyaEqSpC6lNbPBFrKPEduX8NLd91BSuJxHK5dSVPhWrfPtjwW8HEfnRizHMz87ns0tWFbc/kq1R6OOPIJ3je7PC6u2sHLTLuav2crkUf3zXVbbqFrixBHNNmHQlCRJXUaTs8H+fgklbGdyZln1iOWksIoeoTLZMQsE2BF7sSA7rnrE8sXsWPZwsI/N/kp1dBdMGc4Lq5KR+l/OXdf5g+beimT22V798ltPJ3NI62h2JK6jKUlS19Ha/srhYSNTQnIL7NRMGWMzb9Q754ZYwtzshOoey6VxOJUUNFqD61eqo3trXyVTr/kjO/YeoHePAuZ843T6FHbS8anrxsKuTfBPT8HQE/NdTbt0WNfRlCRJam+a66+8+v6XmBhW1Vq/cnDYVu88y7LDkvUrs+OZEyfwWhxIVUws6dWd7Fv7G7y+/ZXqLHr1KOAj7ziKu15Yy+59lTy06A3+buqIfJfVNkpGJkFz6xqDZsoMmpIkqcNoTX9lL/YwascrvHL3/9K/cDmPVZbRp3BPrfPtiwW8FMdUj1bOy45nG30bvX5LZoO1v1KdwQVThnPXC0nv4i/nrevEQXMEvD7PPs02YNCUJEkdQnP9lf2pqNVfeVxYTfc6/ZXbYy/mZ8cn/ZXZUhbGY9hbdy2SBrR2tBIcsVTHdvywYiYM6cvS9Tt4ce027p23jh7dMp3vjyaupdlmDJqSJKldaG1/5ciwgak7l7LjVz/izkwZx/Qsr3fO8tifubm1K+dmJ7AsHk2WTJN1OFopQQiBC6YM56oHFwPwr7/+a/W2oQ38YaXDcubZNmPQlCRJedd8f+VfmRRW1uqvHBgq6p1naXY487LjmZcLl69zJDX7K2Mz/ZXOBisd1LeRCYDWV+zhsjsXVE9s1aG5lmabMWhKkqTDojX9lb3Zw5gdL7H0njsZ0P1V/pAt44jCvbXOtzd2Y1E8pjpULsiOo4I+jV6/Jf2VZ04ayoxJjlZKldnIDY8va3BbJPl3c9WDizlj4pCO/e+jKmhuWwsxQujA76WdMWhKkqQ211R/5bd/v4Qj2VY9Ujk5U8bEsIZuIZvsmPuttiL2zk3YU8rc7HheimPapL/S0UoJ5qzaUuvfSV0RKK/Yw5xVWzr2v5eS4cnX/btg92Y44sj81tOJGDQlSdLb1tr+yjGhnCk7y9h9743cE8oY1XNDvXO+Fo+sng12braU5XEY0f5K6bDYuKPxkHko+7Vb3Qqh71DYUZ5MCGTQTI1BU5IkvS3N9Vd+5/6FnBherTViOSDsqHWObAwsjSOYmx3P/FywLOfgKElJr+5gf6V02Azq2zPV/dq1kpFJ0Ny6Boa9M9/VdBoGTUmS1KzW9Ff2YTfjdixi2T23M7D7ch7PLqdX4b5a59sTu7Mwjq0esVyQHccOejd6ffsrpcNr6uj+DC3uyfqKPbX+zVWp+gPP1NH9D3dp6SsZAeued+bZlBk0JUlSk5rrrxzEllqzwU4IaykIuV9Nc/2VW2Of6t7KedlSXopj2N+CX0Psr5TyoyATuGLmRC67c0G9bXVvSe/wXEuzTRg0JUnq4lrTXxnIckx4gyk7y3jr3h9wbyhjeM9N9c65NjuQuXFC9RqWK+NQ+yulDubMSUO5+eKT6/2BZ1BRIVd95LiOv7RJFdfSbBMGTUmSurCW9FeeFJZVj1hOziyjX9hZ6xyVMbAkjmRutpT52fHMzZaygYO309lfKXVcVX/gufx/F/Dwy+sBuOKciZ0nZIJrabYRg6YkSZ1ca/ori9hF6Y6FvHrP/zC4+zJmZ1+lsLB2SHwr9uDF7FjmxqS/8sXsWHbaXyl1WgWZwMcnH10dNJ9buYWzTjgqz1WlqGpEs2Kda2mmyKApSVIn1lx/5RA21+qvLA3ryNTpr3wzFlX3V87PlvJyHMUB+yulLmXq6AF0ywQOZCPPrHgz3+Wkq/hoCBk4sAd2boC+Q/JdUadg0JQkqQNrbX/luPA6U3aWse/e6/l1ZhlH96z/C+Oq7OAkWMZS5mYnsCoO4eD4Y8Psr5Q6tz6F3ThxeAnz12xl5aZdrK/Yw5DiTrC0CUBBdygaloxobltr0EyJQVOSpA6quf7Ka+5/kXeGpUzJLMv1V5ZRHHbXOseBmOGVOIp52VLm5WaE3URJ9faSXt0Jb+1vcnkD+yulrmH6MQOYv2YrAM+ueJO/OfnoPFeUopKRSdDcugaGT813NZ2CQVOSpHasNf2Vxezk2B0LWHXPrbzUbTl/jK9SWHig1vl2xUIWZMdVj1guzI5lN42PSthfKanK9GOO5Md/ehWAZ17d3MmC5ghYg0ucpMigKUlSO9Vkf+WDizmKTUzJ9VZOzpRRmnmt9gkCbIrFzMmWMj9bypxsKUviSCopaPba9ldKquvkkSX07J5hz/4sz654kxgjobNMnONamqkzaEqSlCet6a/MkKU0rGPyzjIO3Hsdv82UcVTPLfXOuSI7lLnZUubFUuZkJ7A2DsL+SklpKOxWwJRR/fnL8jcpr9jD6s27GX3kEfkuKx2upZk6g6YkSXnQXH/lf92/gKlhSfVssCdnllEU3qp1jv2xgJfj6FrrV26hqHq7/ZWS0jbtmAH8ZXkyidgzr77ZiYKma2mmzaApSVIbaU1/ZT+2M2nHPNbcM4tXui3jj3ElPQora51vR+zFguy46hHLhdlj2ENho9e3v1JS2t5zzJFAGZBMCHTxu0fmt6C0VK+l+RpkKyHTfIuBmmbQlCSpDTTXX3l02MDUcHD9yrGZN2qfIMCGWMLc/9/enUfHed33/f98HywDkMRGEiRBcQcJUDQlcwFkS7ItWZYaO7+4rpM2SVu5TdOktmo3dpLT4zinjSwnjaLGR1Ec66c4rqvIlSxabGzqyItsSd5lWwAI7ssM9x3cAGIhVuK5/WMGQwyxDAZ4yGcGeL/OwQExy32+pkfgfObe773+WjX6tWr0a3XQLZMvL+216a8EcLOsv61MJUX56uy9pl8cuSzfd/KmwwdTpYslL1/yB6TOc/GzNTElBE0AACYhk/7KPA1qrZ1UfVdU/tYn9LIX1cLIlRFjHvJvS4bKRler065S9FcCyCZ5numdq+bptf3n1dY9oAMtHXrb4rKwy5o6Ly8eLtuOx5fPEjSnjKAJAECG0vZXbtuud3r7VWfx2cqN3mGV3NBf2e/ytMetSvZXNvk1aqO/EkAOuLc6HjQl6eeHL0+PoCnF+zTbjic2BLo37GpyHkETAIBRZNJfOU/turOzUae2/C/tz4/pDXdMBYWp/ZUdrji5YU+TX6tdrlp9Khzz+vRXAshW966en/zzm0cu6fffsyrEagJUsVw6Jo44CQhBEwCAG6Trr1xu5+JnV1pMdV5U1d651AFMOuvmqsmvVZNfo0Z/raJuKf2VAKaF1QvmqLIkooudfWo41qqBQV8Feel/v2U9jjgJFEETADDjZNpfuc5OqL4rKtv6uF7xoqqMdIwY86C/NBEq4zOWZ1SZtg76KwHkIjPTPdXz9PLOs+ruH9SuU1dUt2Ju2GVNXfmK+HeOOAkEQRMAMKOk6698fFuT7vH2JWYs4/2Vs60vZYw+l69drjo5Y9nk16hDc5L3018JYLq7t3q+Xgb7W70AACAASURBVN4Z3y37zcOXp0nQZEYzSARNAMC0k0l/ZaWuaENnVGe3/IMO5sf0hjuu/EI/ZbwrbnYiUMZ3hN3rVtJfCWBGu2f19Q/A3jxySZ98cE2I1QSkInEmaMdpaXBAyisIt54cR9AEAEwr6forV9pZ1XlR3ZWYsVzhnU8dwKRTfqUaXXw32AZ/rQ67xXL0VwJA0pKKWVo2d5ZOtnZrx8k29fQPqrgwL+yypmb2AikvIg32SR1npIoVYVeU0wiaAICckkl/ZYGu6W12XHVdUeVtjerbXlTzIp0p4/nOdNAtU+OwGcsWpQ9/9FcCmOnuXT1PJxu6NTDo1Hi8Ve+pSd+bntU8L7589vKheJ8mQXNKCJoAgJyRrr/yiW2Nere3V3VeVPUW0wbvsIqtP2WMXlegnW518vzK7X6NOjUreT/9lQAwMfdUz9eLDackxZfP5nzQlK4HTfo0p4ygCQDIGpnMVkrSArVpU2dU57Z8SbH8qF53J5RXmBoR29ycxExl/JiRvW6lBsb554/+SgCYmOEfpP388OUQKwnQUJ8mZ2lOGUETAJAV0s1WPvbyHlXb6fhusF5U9RbVMu9i6iAmnfQr1ejWqinRX3nUVdFfCQA3wfw5Ea1dVKKDLZ3ae7Zd7d0DKpuV4xvosPNsYHIqaJrZn0h6XNLfOuc+FXY9AIDMZLIbbKEGVNURU/OLWzU7ckjfGTyoikhX6njOdMAtT55d2eTX6LzSb7FPfyUABOOe6vk62NIp56QvvHFID65bmNu/L8sTM5qcpTllORM0zaxe0kcl7Q67FgBA5sbdDfbbB1Siq9rkHYrPVnpRvd2OqMgG4g/0JZnU4wq1w1+d3BG22V+jLvorASA0RQXXV4x85c1j+sqbx5KrUYb/Ts0ZQ0GTGc0py4mgaWZzJL0g6fcl/beQywEAjCLT/soqXVZ9Z1SXXnpGz3pR1UZOybPUiHjJlSb7K5v8Wu1zK3SN/koAyAqv7j2nZ350ZMTtLe29euT5Zj3z8KbcC5tDPZqd56RrfVJ+JNx6clhOBE1JT0v6tnPudTMbN2iaWUTS8FdEyU2tDAAwof7KGjuZ7K+s82JaYpdGjHPUXxRfAuviG/ccc4s0FBPLiws02DMw6vXprwSAW2vQd3rslf2jriBxiv9efuyV/Xpo3aLc+iBv1jypYJY00C21n5bmVYddUc7K+qBpZr8taZOk+gk+5TOSHr15FQHAzJRJf2VE/bqt46B2vfiSSgpjetU/qLJId8p415ynfW5F8uzK7X6NLqp8zOtPZLaS/koAuDUajrWmfKB3IyfpXHuvGo615tYHe2bx5bMXD0htxwmaU5DVQdPMlkr6W0kPOefGfiWnelzSk8N+LpF0OujaAGAmSddfWaoubfZiyRnLO+2oInYt/sDER9tXXUTN/prkjOUOf426VZT22pnOVkrMWALAzXahc2JvzSf6uKxSviweNOnTnJKsDpqSNktaIKnZLPlJdJ6k95jZJyRFnHODw5/gnOuT1Df087DnAQDGkFl/pdMSu6S6zqjaXnpaz3lR1RaN/DzvoitTQ2I32Ea/Vgfccg0qb9w6mK0EgNywoCT9B4WZPC6rcJZmILI9aL4h6Y4bbntW0kFJT9wYMgEAmUvfX7lbt9vx5G6wdV5MVdY6YpzD/uL4ElgXP7/ypFug4f2Vfpr+SnaDBYDccdfKuaoqK1JLe++4O33ftTL9kVNZh7M0A5HVQdM51ylp7/DbzOyqpMvOub2jPwsAcKNM+iuL1KdlHQe058UtKiuM6Xt+VKWRnpTxBlye9rqVifMr4zvCtqp0zOuzGywATC95nunRD67TI883j/jdPmRoNUrO4SzNQGR10AQATF26/spydajOiyVnLNfbMRVaYsFIor+y0xWr2V+TnLHc4a9Wr9Jv+c5usAAwfb1/fZWeeXjTiN/tcyJ5+vy/envuHW0yhBnNQJhzo33+MH2YWamk9vb2dpWWjv1pOwDkqkz7K5faBdVbPFTWe1Gt9s6OGPO8K1ejvzYxY1mrA26ZfHkjHjfcWLOVw89RG69WAEBuGvSdXt3boo9/rVmSdNeKCr30sXtCrmoKetqkJ1bE//yn56TCWaGWE7aOjg6VlZVJUplzrmOiz2NGEwByWLr+ys+9vFvr7FhyN9h6L6qFdmXEOIf829SYWALb4NbqtKvU8P5KR38lAGAMeZ7p/7uzSn/5nWKdudKj3WfaNTDoqyBv/A8os1ZRuRQplfo64rOaC9aGXVFOImgCQJbLpL+yWL1a0blP+1/8msoLYvq+i2lOJHVr+X6Xpz1uVXK2cru/Rm30VwIApqhuRYXO7OxR74CvfWc7tGHp2GcjZ7WhszTP7yFoTgFBEwCy2Lj9ld86oLlqT+mvfJsdV4EN25DbpA5XrO1+Tby/0q/VTletPhWmvTb9lQCATNQtr9DLO+PtGE3HW3M3aErxI07O7+GIkykgaAJAiDLtr1xhLarviqrzpS/oBS+q6qJzI8Y86+am9FfG3JJJ91dyfiUAYKI2L79+lMn2E236vXeHWMxUJTcEImhOFkETAEKSvr9yl+6wI8n+yjovqkpL7cH3nSnqlqgpccxIo79WZzU/eT/9lQCAW6V2UYlKIvnq7LumxuNtcs7JLEc/kOSIkykjaALATZRJf+Us9WpV5x4d3PK8KvIP6TUX0+xIX8p4fS5fu1y1mvzaxFLYNerQnDGvT38lAOBWyfNMG5dX6Cexi7rU1aeTrd1aPm922GVNDkecTBlBEwBuknT9lfPVpjovlpyxXGcnlG/+9QFMuuJmqymxG2yTX6M9bhX9lQCArFWXCJqS1HS8LXeDZkViRpOls5NG0ASAScq0v7LazqquK6burU9pi0W1ouj8iDFP+ZVqdLXJGcvDbrEc/ZUAgBxRt7wi+eemE236jc1LQqxmCoZmNHvapN4OqWjs3dkxOoImAExC2v7KbTv1djusOi+qu7yD2uzFNM86U8bwnemgW5Y8v7LRr1WLrs8qlhcXSPRXAgByyIZl5crzTIO+0/YTrWGXM3mREql4rtTTGl8+u2h92BXlHIImAIwhk/7KOerWms5dim35qublH9Ib7pCKI/0p4/W6Au1y1Wrw16rJr1Wzv0admjXm9emvBADkmlmF+VpXVao9Z9oVO9+l9u4Blc0qCLusySlfRtCcAoImAIwiXX9lpdqSvZX1XlS32wnl2bA4aFKbmzNsN9ha7XUr1a/0/9jSXwkAyGWbl1doz5l2SVLzyTa9d+2CkCuapIrl0rmd9GlOEkETwIyUSX+lyVe1nVV9V1Q9W5/USxbVsqKLI8Y84S9Qk4svgW30a3XUVdFfCQCYcepWVOgff35cktR0ojV3gyY7z04JQRPAjJOuv/LPt+3URosNO78ypgrrShlj0JkOuOVq9K9v3HNB1zdAoL8SADBT1S2fm/xz0/G2ECuZIs7SnBKCJoBpKZP+ylJdVW3nTh3e8pzm58f0hjusokhqSOxxhdrhr07uCLvDX60u+isBABhhUVmRllQU63Rbj3aeuqL+a74K88df4ZOVhoImM5qTQtAEMO2k669cpMsp/ZW1dkreDf2Vl1xpSn/lPrdC1ybwK5P+SgAA4secnG7rUd81X/vOtmvjsor0T8o2w8/SdE4yPgjOBEETQM7JtL+yxk6rviuq/q2f1//1YlpSdGnEmEf9RfElsIkZy2Nuka7PP46O/koAAEa3ecVcbdt5VpK0/URbbgbNsqXx730d8fM0Z80d//FIQdAEkFPS9Vf+xbYdqrODyRnLzV5MZdadMsY152mfW5Hsr2zya3VJZcn7y4sLZD0DKSFyCP2VAACkV7f8erBsOt6m33t3iMVMVuEsafYC6eqF+PJZgmZGCJoAsk4m/ZVl6tLtnc06tuVZ7cqP6Q13RJHItZTxrrqImv01yRnLnf5qdatozOvTXwkAwNTULCxRSSRfnX3X1HSiTc45WS4uPa1YngiaJ6TFG8KuJqcQNAFklXH7K1/Zp9vsouosmpyxrPVOpw5g0kVXpga/Vtv9WjX4tTrglmtQeWmvTX8lAADByPNMG5dX6Cexi7rU1aeTrd1aPm922GVlrnyZdLqRDYEmgaAJ4JbKpL/Sk6+1dlJ1XVENbv2f+qYXU1WkdcSYh/3F148ZcbU66RaI/koAAMJVlwiaUnz5bG4GTY44mSyCJoBbZrz+ygdvX6i/2Nasd3j7kzOWm7xDKrGelDEGXJ72upWJYFmjJr9WrSpN3k9/JQAA2SGlT/NEm35j85IQq5mk8mXx78xoZoygCSBQmfRXVqhD6zubdGLLV7Q3L6Yf6KgKCwdTxut0xWr218SDpavVTr9avYqMeX36KwEAyA4blpUrzzMN+k7bT4xckZQThh9xgowQNAEEJl1/5RI7r3obOr8ypjXemRFjnHflavTXJs6vXKuDbin9lQAA5KBZhfl62+JS7T7drtj5Ll3p7lf5rMKwy8rM0NLZKyc5SzNDBE0AE5ZJf2WeBrXWTqq+Kyp/6xN62YtqYeTKiDFj/m3x3spEf+VpVyn6KwEAmB42L6/Q7tPtkqTmk216YO3CkCvKUNkSSSYNdEtXL0lzKsOuKGcQNAFMSLr+yv/xze16p7c/OWO5yTukOdabMka/y9NuV504u7JGTX6NrqgkeT/9lQAATC91y+fq2TePS4pvCJRzQTM/IpVUSZ1n47OaBM0JI2gCSMqkv3Ke2nVnZ6NObflf2p8X1Q90XAU39Fd2uFnanuivbPTXardbpT6NvWSG/koAAKaXuhWpGwLlpIrliaB5XFqyOexqcgZBE4Ck8fsrH3tln5bbufjZlRZTnRdVtXduxBhn3dyU/sqYWyJfXtpr018JAMD0tLC0SEsqinW6rUc7TrTpG9tPq6q8OLc+JC5fJp38BTvPZoigCcwQmfRX5uua1tkJ1XdFZVuj+pYXVWWkI2U835mibkmyv7LJr9VZzU9bB/2VAADMLIvLinS6rUcDvtMfbd0l6Xr7zfC2l6zFWZqTQtAEZoB0/ZWPf7NR93r7EzOWUW30DmuW9aWM0efytctVJ0Pldn+NOjQneT/9lQAA4Eav7j2nhuMjl8y2tPfqkeeb9czDm7I/bHKW5qQQNIFpIJPZSkmqVJs2dMZ0dss/6GBeTG/ouPIL/ZQxr7jZiQ174jOWe91K+isBAMCEDfpOj72yf9T7nOLvER57Zb8eWrcou98LcJbmpBA0gRw33mzlQ+sW6dGX92qVnVGdF0vOWK7wzo8Y55RfqUZXq+1+jRr8tTrsFsvRXwkAACap4VhrynuCGzlJ59p71XCsNbvfGyRnNE9Jvi956d8fgaAJ5IRMdoMt0DUt7Nirphf/SUUFMX3HHdS8SGfKeL4zHXTL1DhsxrJF6X/B018JAAAm6kLn2CFzMo8LTekSyfKkwT6p67xUmuVLfbMEQRPIcuPtBvu5bx3QbHVrk3dIdV5U9RbTBu+wiq3/+gAm9boC7XSr1ejHZyy3+zXq1KzkQ+ivBAAAQVtQUhTo40KTly+V3ia1n4wvnyVoTghBEwhZpv2VC9WqTZ0xnX/pS/qKd1BrIyeVZ6kRsc3NScxUxmcs97hVGhjnP3f6KwEAQNDuWjlXVWVFamnvHffD7LtWzr3VpWWuYnkiaJ6Ulr0z7GpyAkETCFG6/srPvrxH1XY63lvpRVVvUS3zLo4Y54S/QE2uVk2J/sqjrirZX1leXKBrPQOjXp/+SgAAcLPkeaZHP7hOjzzfnLb9JuuVL5f0U444yQBBE7jJMumvLNSAqjpi2vHiVhUXxPRdF1VFpCt1PGfa75annF95QRVjXn8is5X0VwIAgJvh/eur9MzDm0Z8mF0+q0CP//od2X+0yZDkhkAEzYkiaAI3Ubr+yhJd1aah3WC9mDbYEUVs2OyjSd0uoh3+6uSMZbO/RldVnPbamc5WSsxYAgCA4A19mP0PPzmiJ16NSpL+2dsW5k7IlDjiZBIImsAUZNpfWaXLqu+M6vJL/7+e9aKqjZyWd0N/5SVXmtJfuc+t0LU0/6kyWwkAALJZnmf6nXtW6m9eP6T+a75+Ersk55zMcuS9SHJG82S4deQQgiYwSWnPr9y2RzV2MtlfWefFtMQujRjnqL9ITX6tmlyNGv21OuYWaSgqlhcXaDBNfyW7wQIAgFxQXJind6ycq58euqRz7b06dKFLNQtLwi5rYsoTM5rtpyV/UPLywq0nBxA0gXFk0l8ZUb+WdBzU7he/rtkFMX3fRVUW6U4Z75rztNetSPZXbvdrdUllY16f3WABAMB0cn/tAv30UPyD9x9FL+RO0CxZJHkFkj8gdZyVypeGXVHWI2gCY0jXX1mmTm32Yqr3YqrzorrDjipi164PYNJVF1GzvyYeLF2tdvqr1a30Z0WxGywAAJiO7qup1J8n/vzj2EX9p/dUh1rPhHl58XDZejTep0nQTIugiRkrs/5KpyV2SfWdB3XlpS/qq15UNUVnRox50ZWpIbETbKNfqwNuuQY1/tIK+isBAMBMUV05W0sqinW6rUeNx9p0te+aZkdyJJKUL0sETfo0JyJH/l8FgpX2/Mptu3W7HY+fXZnor6yy1hHjHPGr4ktgXfz8yhNuoYb3V/r0VwIAACSZme6rqdQLb51U/6CvXxy5rAfXLQy7rIkZ6tPkLM0JIWhi2sq0v3JZxwHteXGL5iT6K0sjPSnjDbg87XUrE2dXxneEbVXpmNenvxIAAGCk+2sX6IW34rOCP45dzKGgyc6zmSBoYlpK119Zrg7VJXor672o1tsxFdrg9QFM6nTFavbXJGcsd/ir1atI2mvTXwkAADC2u6vnqSDPNDDo9KPYhdw55qRiRfw7Z2lOCEETOSnT/sqldkH1nVFdeenv9LwX1eqisyPGbHEVyd7KJr9WB9wy+fLGrYP+SgAAgMzMieSrfsVc/fzIZZ1q7dGxS1e1qnJO2GWlx4xmRgiayDnp+isf27Zb6+xY8vzKei+qhXZlxDgx/7b4+ZV+jRrcWp12lRreX+norwQAALgp7qup1M+PXJYk/Sh6MUeCZqJHs+OMNDgg5RWEW0+WI2giK2XSX1msXq3o3Kf9L35NpQUxveZimhPpTRmv3+Vpt6tOzFjWaLtfoysa+9wm+isBAABunvtqK/X4dw9Kivdp/u67VoZc0QTMWSDlF0nXeqX2U9LcVWFXlNUImsg64/dX7tc8tWtzcjfYqNbbceWbf30AkzrcrOSGPU1+jXa5avWpMO216a8EAAC4+WoXlmhRaZFaOnr1y6OX1TswqKKC8Y+EC51ZfPnspVh8+SxBc1wETdxymfZXrrAW1XdF1fXS3+prXlSrilpGjHnWzVWjvzbZXxl1S+TorwQAAMhKQ8ecfL3plPqu+frl0cu6v3ZB2GWlNxQ0OeIkLYImbqnx+isfvH2hPrttl+6wI8nZyjovqkrrSBnDd6aoW5KcrWz01+qs5ifvLy8ukOivBAAAyGr318aDphTv08yNoJno02RDoLQImghcJv2Vs9Wj6s49im55XuV5Mf1AhzQr0pcyXp8r0E5XnQiVtWr216hDYzeM018JAACQ/e5ZPV95nmnQd/pJ7GLY5UxMxVDQZEYzHYImApWuv3K+2lTnxXSXd1B1XlTr7ITyzKWMccXNTvZXNvhrtdetVL/S7+pFfyUAAEDuKCsu0OZlFWo43qqjl67q5OVuLZs3K+yyxscRJxNG0ERGMu2vrLazquuKqXvrU9piUa0oOj9izNNuvhr8tckzLA+7xfRXAgAAzAD31Vaq4XirJOnHsQv6yN0rwi0onaGls/RopkXQxISl66987Js79XY7nDy7ss6Lap51pozhO9NBtyxxxEg8WJ7T9VlF+isBAABmjvtqKvXX34tKih9zkjNBs6tFGuiVCorCrSeLETSRIpP+yjnq1prOXYpt+aoq8mL6gQ6rONKfMl6vK9BOtzq5G2yzv0adGntJBP2VAAAAM8e6qlLNnxPRpa4+/fzIZfVdG1QkP4uPOZk1VyqcI/V3xc/SnL8m7IqyFkETSeP1Vz72rf1aoNbkTGW9F9VaOzmiv7LVzUnZDXavW6mBCbzM6K8EAACYeTzP9J6a+fpG8xl19w/qC28c0rtWV2bvJMLQWZoX9sc3BCJojomgOYNk0l9p8lVtZ1XfFVXP1if1fy2qpUUjdwM74S9Qk4svgW30a3XUVdFfCQAAgAmrmFWY/PPTPzyip394JNmeNbwtKmuUL48HTfo0x0XQnCHS9Vd+7ps7tNFiw86vjKnCulLGGHSm/W55YsYyHiwvqCJ5P/2VAAAAyMSre8/pf//s2IjbW9p79cjzzXrm4U3ZFzbZeXZCCJrTSCb9laW6qtrOHTqy5R81Ly+mH+qIIpHUkNjtItrhr07OWO7wV+uqise8Pv2VAAAAmKhB3+mxV/bLjXKfU/w95GOv7NdD6xZl13tFztKcEILmNJGuv7JKl1J2g6210/Ju6K+85EqH9VfWap9boWv0VwIAAOAmaDjWmvKe8UZO0rn2XjUca82u945zV8W/n9ku+b7kjd82NlMRNHNEpv2VNXZa9V1R9W/9vL7hRXVb0eURYx71F8WXwLr4UthjbpGuzz+Ojv5KAAAABOFC59ghczKPu2VW3icVlsSXzp54U1r57rArykoEzRwwXn/l+9Yu1Oe+2aw6O5jSX1lq3SljXHOe9roVyf7KJr9Wl1SWvL+8uEDWMzDq0gX6KwEAABC0BSUTO4Nyoo+7ZQpnSet/XWp+Ttr5AkFzDATNLJFJf2WZunR7Z7OObXlWu72ofmhHFYlcSxnvqouo2V+TnLHc6a9Wt8b+j5T+SgAAANxKd62cq6qyIrW094472XHXyrm3urT0Nj4cD5r7X5Y+8D+lotKwK8o6BM0sMG5/5Sv7tFgXVedFdddQf6V3esQYF12ZGhMzlQ1+rQ645RpU+sNu6a8EAABAGPI806MfXKdHnm8eMdmhxM9D7VlZZ0m9NG+NdPmQtO+b0uZ/H3ZFWcecG+3zg+nDzEoltbe3t6u0NJxPGjLpr/Tkq9ZOpWzcs9haR4x5xK+KB8vEjrAn3EJNtr9y+LbR49UKAAAABG20SRdJWldVqu98MouXpf7sb6TXPystfYf0H78fdjU3TUdHh8rKyiSpzDnXMdHnETRvsvH6Kx9Yu1Dvffw7WtJ9MBksN3kxlVpPyhgDLk973Uo1+rXantgRtlXX/7eUFxeofRL9lVl9EC4AAABmjKHJjpb2Hv2P7xzQpa5+SdK3/+BdetvisjTPDknHOelv1knOlz7RJM1fE3ZFN8VkgyZLZ6cok9lKSapQh9Z3NunElq9onxfTD+2oCiODKWN2umI1+2uSM5Y7/Wr1KjJmDfRXAgAAIJcNb8/q7LumP3t5nyTpSz8+qi/8641hlja20ipp9YPSoe/HNwV68LNhV5RVmNGcgvFmKx9at0h3/+XrKrp6SvUWTc5YrvbOjhinxVXEN+1JfB10y+Qr/Xk8Q7OVP/v0A3ptf8uYtTBjCQAAgFzR0z+oe5/4gVqv9ssz6cf/9b1aOndW2GWNbv/L0kv/Tiqpkv5wn+Sl3yMl1zCjeZNkshtsngY1t+OAfvniNnkFh/SKO6CFkSsjxoz5t10Plq5Wp12lOL8SAAAAkIoL8/Q796zQk6/F5Dvpyz89qs99aH3YZY2u5gNS8Vyp85x05AfSmofCrihrEDTHMd5usJ/71gEVqVcbvCPJGctN3iHNsWFNzCb1uzztdtWJsytr1OTX6IpKkg/h/EoAAAAg1b+7e7n+/sdH1N0/qK83ntIfvG+N5s8Zu5UsNPmF0p2/Kb3199KO/0PQHGZGL53NtL9yvtq1edhusOvtuPLNTxmzw81KBMr4jOVut0p9Khyzvj98sEZPvR6TNP6OsOwGCwAAgJnkz7+1X1/52TFJ0n95YLX++J/VhlzRGM7tlr70bimvUPrjqDQrC8/9nAKWzmYoXX/ln23bqxV2Lt5bmZixXOW1jBjnjJuXPL+y0a9VzC2RS/RXlhcXqL9nYNTrc34lAAAAMLb/+K6Veu7nx3XNd/rqL07oo/dVa04kC+NL1Z3Sojuklj3Snq3SOz4adkVZIQv/n7o5Go626r13lozZX5mva6rs2KfGF7+h/IKYvu0OqjKSGth9Z4q6JclQ2eTX6qzmj3nNiewGS38lAAAAMNLi8mJ9aMNt+qfm02rvGdCWhpP6vXevCrus0W14WHr109KO56dV0Bz0nRqOtk7quTMmaP7uc426bcGRZH/lLPVoo3c4vgzWotroHdYs67v+BJP6XIF2uurkUtjt/hp1aE7aa2U6WykxYwkAAADc6GP3rdI/NZ+WJH35J0dVu7BErd392Tcxc+dvSq/9d6lld3wpbdWdYVc0ZUMrQM9cIGiOa76uaGPnbp1/6Uv6shfVusgJ5Vlqf+oVNzulv3KPW6V+FYw7LrOVAAAAwM2xZmGJHrx9oV4/cF7nO/v0kf/dkLwvq47ymzVXqv1A/LiTnS/kTNDM5ISNTM2czYD+pESlkdRwd8qvVKOr1Xa/Rg3+Wh12i1P6K9snsRtsVr3gAQAAgBz39A8P6a+/Fxtx+42bZ4Yu9n3pa/8qftzJH0fjO9JmsfFP2Nivro42bfIOaf3AHn3681skNgMane9M+/zlyd7KJr9GLRp7qepE+ivfv75Kv7Ke2UoAAADgZhj0nZ7/5clR73OKvzd/7JX9emjdovDfg1c/IM1ZJHW1SLHvSus+FGo5mZ6wsVCt2tQZ0/mXvqSveAe1NnJSeebU4Zw+PYnrz5igeU/fF9Q9zsY9Q9gNFgAAAMgODcdaU96L38hJOtfeq4ZjreG/J8/Ll97+29KbT0k7XrjpQTNdkBzrhI0Hb1+oP9u2R9V2OnlsY71Ftcy7OOIaJ/1K/XhwlaTXMq5vxgTNLs1KLIq9jv5KAAAAIHtd6Bw7ZE7mcTfdqMhssQAAETpJREFUxofjQfPwa1LHOan05izpHS9IShoxW1moAS3uiGrHi1sVyY/pe4qqItKVMuagMx1wqStAz2uu/IFuETQnYLz+SnaDBQAAALLHgpKiQB93081fIy19h3TqLWn3FuldfzipYTJd9ipJLe29+tjzzSovLlCZOrXZi6nOi2mzF9Pb7agiNpDy+B5XqB3+ajW5+Gaozf4adWnWpOodzYwKmvRXAgAAALnjrpVzVVVWpJb23jF3QK0qi7+Pzxob/m08aO54Qbr3U5Jlli/Gm618aN0iPfbK/hv+LpxWWEs8VFpMdYMxrSk6M2Lci640OVPZ5Ndqn1uha2ni4FSS0YwKmvRXAgAAALkjzzM9+sF1euT55hFtb0MWlxcpfpJGlkwYve3D0nc/LV0+JD37AWnjR+L9mpE5yYdkeqxIS3uvHnm+WZ983xp1trfqbu+YNtphbfAOa6N3SJU2cjPYI36VGv1abXc1avRrddwt0tDfUXlxgQZ7BkY8Rxq5AvTMhe5J/TXMmONNXttxTO+9czkzlgAAAECOGW2Wb7gPb7xNT/zGndp+oi07Viu++QXp9Ucl58d/LpwTD6AbH9ar7cv02LcOjHqsyJ9/O/X2fF1TrZ3WBu+wNthhvd07otV2Vp6lZrg+l6/dbpW2J2Yst/tr1KbSMcv7wwdr9NTr8SNjRtuzZujImEHf6Ye7T+ihjSulDI83mTFBs729XaWlY/9lAwAAAMheN84Cdvdd08de2K6BwXieKS7IU8/AYPLxYZ1vP1Rnx4UTWnv+21p28huy1qPJ+4/6VXpp8D59Y/DduqxSLVSbbrNLqrLLWmKXtNguabFd1mK7rBXWoiIbOfN42s3XTr9aO/3V2uGv1h63Sv0qSFvb0Gzlzz79gF7b3zLmEt3hf2cdHR0qKyuTCJqpCJoAAADA9PTa/vP66P9pkj9KpLlxdu5WGLW/sjSip+7uUf2V76h31zc0S/H7Bp3JyZRv/rhjdrhZ2uWv0k63Wjv9au231To3WDbqY01S2awCtXfHw+l4s5XS+JsOJa9P0BwdQRMAAACYngZ9p81/8ZqudI/fb/izTz+gPM8mFKwmcs1M+iuHRv+371imb74V1a/mvaXfzPuR6r340tUBl6dzbq7Oar7OuPk64+bprJuvs26eTroFOuEWyg07qHEiy14lTWi2ciImGzRn1GZAAAAAAKaPhmOtY4ZMKR7EzrX3quFYq9p7+tOGr3RBdKwdYYf6K0ebwhu67fm3Tkoq1tbB+7V18H4tVKucTJdUJn9YkBzLUGj+xAOrVbtozog6btz49KF14Z6wwYwmAAAAgJz08s4z+uSWnWkfd9eKCjUcbxtxeyazgGPNWN4MN+6wO9llr0Fg6ewYCJoAAADA9PSLI5f1r7/8yymNMbyvcbRlr07SI/ev0gu/PKWO3rFnTyeiqMBT78DoPZk3HisSxLLXIBA0x0DQBAAAAKanQd/pXU/8QC3tvbdkpnGqMjlWJMxlr8NNNmimXwwMAAAAAFkozzM9+sF1kq6HtSGW+PqNTbfd6rJGMMVnJT/xwGo98/AmLSorSrl/UVlRyrLYPM90d/U8fWjDbbq7el5454FOATOaAAAAAHLaWJv0PPrBdSorLpzy8tpMZFN/ZRDYdRYAAADAjPT+9VVj7rI66DtVlRVNeXnt3NmFarvaP+oY4/VX3rgbrHR9xnI6Y0YTAAAAwLQ2tGOsNHK20UkqH2MzoKHHDIXIj39t9DGk7OyvDMK07NE0s8+YWaOZdZrZBTPbZma1YdcFAAAAIHe8f33VmL2Rf//wJv3Vr98hafQ+T0l69IPr9Kt3jj3GdOuvDEJWz2ia2auStkhqVHyZ719KWi9pnXPu6gTHYEYTAAAAwLizjeP1eQ5f9jrdZizTmRHHm5hZpaQLku5zzv1kgs8haAIAAABIa6aFyImYKZsBlSW+t471ADOLSIoMu6nkplYEAAAAYFqYCZv03CpZ3aM5nJl5kp6S9KZzbu84D/2MpPZhX6dvQXkAAAAAgIScCZqSnla8P/O30zzuccVnPoe+ltzkugAAAAAAw+TE0lkz+6KkX5P0HufcuDOUzrk+SX3DnnuTqwMAAAAADJfVQdPiKfHvJH1Y0v3OuWMhlwQAAAAASCOrg6biy2X/jaQPSeo0s0WJ29udcz3hlQUAAAAAGEu292g+onif5Y8knRv29Vsh1gQAAAAAGEdWz2g652iwBAAAAIAck+0zmgAAAACAHEPQBAAAAAAEiqAJAAAAAAgUQRMAAAAAECiCJgAAAAAgUARNAAAAAECgCJoAAAAAgEARNAEAAAAAgSJoAgAAAAACRdAEAAAAAASKoAkAAAAACBRBEwAAAAAQKIImAAAAACBQBE0AAAAAQKAImgAAAACAQBE0AQAAAACBImgCAAAAAAJF0AQAAAAABIqgCQAAAAAIFEETAAAAABAogiYAAAAAIFAETQAAAABAoAiaAAAAAIBAETQBAAAAAIEiaAIAAAAAAkXQBAAAAAAEiqAJAAAAAAgUQRMAAAAAECiCJgAAAAAgUARNAAAAAECgCJoAAAAAgEARNAEAAAAAgSJoAgAAAAACRdAEAAAAAASKoAkAAAAACBRBEwAAAAAQKIImAAAAACBQBE0AAAAAQKAImgAAAACAQBE0AQAAAACBImgCAAAAAAJF0AQAAAAABIqgCQAAAAAIFEETAAAAABAogiYAAAAAIFAETQAAAABAoAiaAAAAAIBAETQBAAAAAIEiaAIAAAAAAkXQBAAAAAAEiqAJAAAAAAgUQRMAAAAAECiCJgAAAAAgUARNAAAAAECgCJoAAAAAgEARNAEAAAAAgSJoAgAAAAACRdAEAAAAAASKoAkAAAAACBRBEwAAAAAQKIImAAAAACBQBE0AAAAAQKAImgAAAACAQBE0AQAAAACBImgCAAAAAAJF0AQAAAAABIqgCQAAAAAIFEETAAAAABAogiYAAAAAIFAETQAAAABAoAiaAAAAAIBAETQBAAAAAIEiaAIAAAAAAkXQBAAAAAAEiqAJAAAAAAgUQRMAAAAAECiCJgAAAAAgUARNAAAAAECgCJoAAAAAgEARNAEAAAAAgSJoAgAAAAACRdAEAAAAAASKoAkAAAAACBRBEwAAAAAQKIImAAAAACBQBE0AAAAAQKAImgAAAACAQBE0AQAAAACBImgCAAAAAAJF0AQAAAAABIqgCQAAAAAIFEETAAAAABAogiYAAAAAIFAETQAAAABAoHIiaJrZx83suJn1mtlbZnZX2DUBAAAAAEaX9UHTzH5L0pOSHpO0SdIuSd8zswWhFgYAAAAAGFXWB01JfyTpy865Z51z+yV9TFK3pN8NtywAAAAAwGiyOmiaWaGkzZJeH7rNOecnfr47rLoAAAAAAGPLD7uANOZLypN0/obbz0taO9oTzCwiKTLsphJJ6ujouBn1AQAAAMC0Ndkcle1BczI+I+nRG29cunRpCKUAAAAAwLRQImnCqTPbg+YlSYOSFt5w+0JJLWM853HFNw8aUiLptKQlkjqDLhA5j9cHxsJrA+Ph9YHx8PrAWHhtYDzZ/PookXQ2kydkddB0zvWb2XZJ75O0TZLMzEv8/MUxntMnqW/oZzMb+mOnc471s0jB6wNj4bWB8fD6wHh4fWAsvDYwnix/fWRcT1YHzYQnJT1nZk2SGiR9StJsSc+GWhUAAAAAYFRZHzSdc183s0pJn5O0SNJOSe93zt24QRAAAAAAIAtkfdCUJOfcFzXGUtkJ6JP0mIYtpwWG4fWBsfDawHh4fWA8vD4wFl4bGM+0en2Ycy7sGgAAAAAA04gXdgEAAAAAgOmFoAkAAAAACBRBEwAAAAAQKIImAAAAACBQ0zpomtnHzey4mfWa2VtmdlfYNSE7mNl7zOwVMztrZs7M/kXYNSE7mNlnzKzRzDrN7IKZbTOz2rDrQnYws0fMbLeZdSS+fmFmHwi7LmQfM/uTxL8vT4VdC8JnZp9NvB6Gfx0Muy5kDzO7zcyeN7PLZtZjZnvMrC7suqZi2gZNM/stSU8qvkXwJkm7JH3PzBaEWhiyxWzFXxMfD7sQZJ37JD0t6Z2SHpJUIOn7ZjY71KqQLU5L+hNJmyXVSfqBpJfN7G2hVoWsYmb1kj4qaXfYtSCr7JNUNezrXeGWg2xhZhWS3pQ0IOkDktZJ+mNJbWHWNVXT9ngTM3tLUqNz7hOJnz1JpyT9nXPur0ItDlnFzJykDzvntoVdC7KPmVVKuiDpPufcT8KuB9nHzFol/Vfn3FfCrgXhM7M5kpol/WdJ/03STufcp8KtCmEzs89K+hfOuQ1h14LsY2Z/Jele59y7w64lSNNyRtPMChX/tPn1oducc37i57vDqgtATipLfG8NtQpkHTPLM7PfVnyFxC/CrgdZ42lJ33bOvZ72kZhp1iRado6a2QtmtizsgpA1/rmkJjPbmmjb2WFmvx92UVM1LYOmpPmS8iSdv+H285IW3fpyAOSixEqIpyS96ZzbG3Y9yA5mdoeZdUnqk/T3iq+I2B9yWcgCiQ8eNkn6TNi1IOu8Jel3JL1f0iOSVkr6qZmVhFkUssYqxV8XhyT9iqRnJH3BzP59qFVNUX7YBQBAFnta0nrRR4NUUUkbFJ/t/peSnjOz+wibM5uZLZX0t5Iecs71hl0Psotz7rvDftydaPE6Iek3JbHsHp6kJufcnyZ+3mFm6yV9TNJz4ZU1NdN1RvOSpEFJC2+4faGklltfDoBcY2ZflPRrkt7rnDsddj3IHs65fufcYefcdufcZxTfWOyTYdeF0G2WtEBSs5ldM7Nrim8u9geJn/PCLQ/ZxDl3RVJM0uqwa0FWOCfpxg8rD0jK6eXV0zJoOuf6JW2X9L6h2xJL4N4n+mgAjMPivijpw5IecM4dC7smZD1PUiTsIhC6NyTdofhs99BXk6QXJG1wzg2GWBuyTGLTqGrFAwbwpqQbj1KrUXzWO2dN56WzTyq+nKlJUoOkTym+YcOzoVaFrJD4BT/8U8SVZrZBUqtz7mRIZSE7PC3p30j6kKROMxvq6253zvWEVxaygZk9Lum7kk5KKlH8tXK/4j01mMGcc52SUnq5zeyqpMv0eMPMPi/pFcWDw2LFj98blPRimHUha/yNpJ+b2Z9KeknSXZL+U+IrZ03boOmc+3riWILPKb4B0E5J73fO3bhBEGamOkk/HPbzk4nvzynerI+Z65HE9x/dcPt/kPSPt7QSZKMFkr6q+Bl47Yqfk/grzrnXQq0KQLZbonionCfpoqSfSXqnc+5iqFUhKzjnGs3sw5Iel/Rnko5J+pRz7oVwK5uaaXuOJgAAAAAgHNOyRxMAAAAAEB6CJgAAAAAgUARNAAAAAECgCJoAAAAAgEARNAEAAAAAgSJoAgAAAAACRdAEAAAAAASKoAkAAAAACBRBEwAAAAAQKIImAAAAACBQBE0AAG4xM6s0sxYz+9Nht91jZv1m9r4wawMAIAjmnAu7BgAAZhwz+1VJ2yTdIykqaaekl51zfxRqYQAABICgCQBASMzsaUkPSmqSdIekeudcX7hVAQAwdQRNAABCYmbFkvZKWipps3NuT8glAQAQCHo0AQAIT7WkxYr/e7wi3FIAAAgOM5oAAITAzAolNSjemxmV9ClJdzjnLoRaGAAAASBoAgAQAjP7a0n/UtLbJXVJ+rGkdufcr4VaGAAAAWDpLAAAt5iZ3a/4DOZHnHMdzjlf0kckvdvMHgm1OAAAAsCMJgAAAAAgUMxoAgAAAAACRdAEAAAAAASKoAkAAAAACBRBEwAAAAAQKIImAAAAACBQBE0AAAAAQKAImgAAAACAQBE0AQAAAACBImgCAAAAAAJF0AQAAAAABIqgCQAAAAAIFEETAAAAABCo/wfP4M1rLjmzDAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(11, 7), dpi=100)\n", "plt.plot(x,u, marker ='o', lw=2, label='Computational')\n", "plt.plot(x, u_anal, label='Analytical')\n", "plt.xlim([0, 2* np.pi])\n", "plt.ylim([0,10])\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.title('Burgers Equation at t=10');\n", "plt.legend();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Results\n", "\n", "Looks pretty Cool! But I would like to know how this evolves over time, not just at its final timestep. As always we will attempt to animate the wave and take a look at the results.\n", "\n", "## Animating the wave moving" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5oAAAJmCAYAAAA90v9dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzs3Xu4VmWdP/73DYhKCKiBWKFomUlqZdqMdFDLwjyMNebYN02tqSkt+zZN37LpoGb9KsvKDtbYoKY1dtbM0jIUG800My1LHTMJG/GMKBqmsH5/rGfTw2bvzQZu2Bt8va7rufC51+mz1rP29nnv+15rlaZpAgAAALWMGOoCAAAAWL8ImgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAyolNKUUo4f6jq6DceahkIp5cjOsZhaeb2zSymza65zJba9XSnlJ6WUBZ19e9VQ1AHA6hE0gSq6vvB2v+4upVxaSnnlUNc3lPo4Lt2vLw91fUlSStl3uAW34VjT+qSUMq2UcnztkFrBV5PslOT9SV6f5JqhLWdgpZR/rx2GSykTSimnlVLuKaU83Pk9ussglz26lHJkzXo66626n6WULUspH+/s20Od34d7DjD/9FLK5aWUR0opd5ZSPldKGVurHqC+UUNdALDe+VCS25KUJFskOTLJj0opBzRNc8FQFjbELk5yVh/t/7O2C+nHvkneluT4PqZtnOTxtVpNazjWtD6ZluS4JLOTzOk17RVru5gkKaVsnGT3JB9tmuYLQ1HDKvj3JN9Jcl6NlZVSRiT5YZLnJPlkknuTHJ1kdinl+U3T3LKCVRzdWebMGvV0qbqfSbZP8t4ktyT5bdrPvU+llOcmmZXkxiTvSvK0JO9Osl2SJ/QfMmE4EzSB2i5smmZpD0QpZWaSu5L8nyRVgmYp5UlN0zxcY12D2NaIJKObplm0mqv6n6ZpvlajprWtwr5XNxxrWp80TfPXIdr0xM6/D6xoxrX5e2Ate02S6UkObprmO0lSSvlW2j9KnZDkdUNYW02/SrJ50zT3l1Jek+TbA8z7/yWZn2TPpmkeTJJSypwkXymlvKJpmp+s8WqBlWboLLCmPZDkL+nqfSql7NnXMKlSytRO+5FdbWeWUhaWUp5eSvlRKeWhJF/vmv62UsofSyl/KaVcXUp5cV/Xl5VSNiylnFBK+UMp5dFSyu2llJNKKRv2mq8ppXyhlHJoKeV3SR5Nsk9n2mtLKb/qDPN6sJTy21LK/611oDrb+JdSyq0D7U9/1+X1dVw7y3+7lDK3a78/0+k56pnnzLQ9h8sM8+11TI7vta3nlVIu7ByHhaWUWaWUv+81T0+dLyylfLprGOC5pZSJGcDK1tQZAtqUUp5ZSvlaaa/vu6eUcmJpTSmlfL9T752llH/rY5uDOkcGqPnvSikXdbb9SCnlslLKC7umv6ZT4x59LPuWzrQdu9peWkr5784xe6BT/w6DqKPP61dLKXM6xzWdn7GeL/aXdh3jPTvT+/oZmlRKmVlKuauUsqiUcn0p5Yhe8/T8DL+761x+tJTyy1LKbiuo+/gkf+q8/WRnPXN6pnXeTyul/FcpZX6Sy1fmWNU4R/qouUnypCRHdB3DM1e03Aq8Ju0f577X09A0zT1JvpXkwIHOx87xenaSPbrqmd01fdvS/j64v3OO/qKUst+KCloT+9k0zUNN09w/iG2PS/LyJF/rCZkdZyVZmOSfVqcOYM3RownUNr6U8uS0Q2cnJTkmydgkq9ObNyrJj9N+sXx3kkeSpJRyVJIvJPnvJJ9JMjXtsK75Sf7cs3BpeyXPT/KiJKelHX61U5J/TfLMJL2vO3pp2i8vX0g7BG1OKeXlSc5JO3zrvZ35dkjywiSnDGIfNuocl94e7Ok9KqX8c5L/SPLzJJ9Nsm2n7vuT3D6IbfTl4CRjknwpyX1JXpD2M3laZ1o623xK2i9zr1/RCkspz057zB9MclKSx5K8Je3Qvj2aprmq1yKfT/uZnJD2M3pn2mN7yACbWamaunwz7ed7bJL9knwg7fF7S5JL0n52hyb5VCnll03T/KyzTyt7jiyjlPLSJBem7aU5IcmSJG9Ickkp5cVN01yddjhkzxfjy3qt4pAkv2ua5obO+vburO+PaYcOb5z2c7uilLJL0zRzVuKY9OVnST6X5B1pe4tu7LTf2NfMpf3DxOwkz0j72d2W9vw5s5QyoWma3j8Dr0uySdrPsUnyniTfK6Vs2zTNY/3U9L20f5j6TNqftR+lPV7dvp12qOW/p/0dsyrHapXOkX68Psl/Jrk67XmTJLd26togyfgBlu12f9M0Szr//bwk13a973F1kn9Jez7+tp/1vDPtz9vCJB/ttN3VqWeLtL9bxqT97O9LckSS80spr2ma5twB6lsT+zlYO6X9f8Ay1+o2TfPXUsp1aY8XMBw1TePl5eW12q+012I2fbwWJTmi17x7dqbt2at9aqf9yK62MzttH+s17+i0IfDqJKO62o/ozD+7q+2wJIuTvKjXOt7SmXd6V1vTmXdar3k/m2RBkpGrcGz6Oi49r9d25tkg7RfCX6cdqtuz7Jv72J+eYz11Rcc1ycZ91HNs2iC0VVfbF9r/JfRb//Fd789N29O7bVfblmmD52V91HlxktLV/um0PdzjV3DcVqam4ztt/9HVNjJtQF+S5L1d7RPS/rHizFU5R/qopaQd1nhRr/3cOG34+UlX2391PueRXW2TO9v+YFfbrzvzbdbVtnNnvq8OdC70PjZd7XN67fNrep8vXdNm9zrn/m9n3kO72jZIG1weSrJJr5/he5Ns2jXvP3Ta91/BZ96z/Lt7tfd8vv/VxzKDPVardY4MUPPCvubL334eB/Oa2mt9M/tY376deWesoJ4buj+7rvbPdJZ/UVfb2M45eluSEWtzP3utY6BzsWfai/uY9q0k81b0GXl5eQ3NS48mUNvb8rcb3GyR9gv8f5ZSHmqa5nv9L7ZCX+r1ftckmyd5X9M03TeF+XraL1TdDk7bg3FTr17FSzr/7pX2C3OPy5qm+X2vdTyQdujYy9MGipX1/bTBqbeenold0/YAf6hZ9vq4M9PeEGSVNE3zl57/LqU8KW34+XnacPS8JHNXZn2llJFpbxRzXtM0f+zazrxSyn8leXMpZVyz7BC305qmabre/3fansKtk/xmJXdpRf6zq6bFpZRr0vbezuxqf6CUcnPaHuMeK3uOdHtu2puSfCTJ5qWU7mmzkry+lDKiaXtyvpn2euU9O9OS9ov0iM60lFK27KzzpKZraGHTNL8ppVycNnCsbfsmuTNtT2NPPY+VUj7Xadsjy16D/c2maeZ3vf/vzr/dx3xVLHOX5lU8Vqt6jqys69P+vhiMO7v+e+O0f8jpbVHX9FWxb5Krm6ZZOuS4aZqFpZTTknws7c2hbliF9a7qfg5Wz/72d0xW9XgAa5igCdR2dbPszYDOSdvj8IVSygXNqt1k5PF0DYXt2Lrz7x+6G5umebznuq4u26Ud5npPP+uf1Ov9bX3Mc2raIY8XllL+N8lPknyraZrBhs4/N03z0wGm9+zPMneU7HyZ/2Mf8w9KKWWrJB9O26O0aa/Jgx3u1m1i2qF3N/cx7ca0gWlKkt91tfcOsz0BpHc9NfTe1oIki5qmubeP9s273q/sOdJtu86/Xx1gnvFp9/uizrYPyd+C5iFJrmuapucPND3nQn/HeEZZ+zfC2TrJLc3ywx5v7JrebZnPoWma+Z0Avrqfee+fzVU5Vqt6jqyUTtAe6Ge+P39J0td1mBt1TV8VWyfpPaw9WfYzXOmguRr7OVg9+9vfMVnV4wGsYYImsEY1TbOklHJp2qF326UNIE0/s4/sp/3RPr7growRaXsO39XP9N7XPy73xaVpmrtLe4v9GWlvp//KJG8opZzVNM0Rq1HbqhjU8ev0Pl6cZLMkn0hyU5KHkzw1bU/p2roh3OJ+2ks/7bW3NZjtr+w50q3nOP6/JNf1M8/CJGma5tFSynlJXl1KOTptr/8L015zuKb19/O1Jqypz7xGqFjVc2SllFJGp/3ZG4x7mqbpqWFe2qHovfW03bGqNa0Jq7GfgzWv829/x2RYHQ/gbwRNYG3o+V3T83Dtnh6tCb3m690rMpCeu1M+I8mlPY2llFFpr/PqHpJ5a9pn0s3qNYRzpXR6Y3+Q5Aedm8ecmuQtpZQTm6b5w8BLr1DP/myXvw3X7LnRxjZph6f1GOzx2yntjUOOaJpm6TM8Ozc26m2wx+WetNeubd/HtGelvdZtVW9ctKo11bA658itnX8fXEGvdY9vpr2W+GVpe1FLp61Hz7nQ3zG+dwW9mfPT69zohIHeX9RXZj//lGTnriHA3fX0TB8Kq3usaujvOE5P1++mFdgmf3uW6XVJXtzHsf67tD97K3r2bn/1/Cn9H6ee6auy3lXdz8G6Ie2oll3TXpOZZOk5/dzuNmB48XgTYI3qBKVXJPlr/jZE609pexBe0mv2o1di1dekvWvimzvhssehWX543rfS9uK9uY/6Nu5cuzigUsoyQ+g6XwB7wuygHn+xAtekDXFv7XyB6nFklg+UPcFm6fHr9F7+S6/5enoOStd8JW3vcm8Pd6b33tYyOr0RP0n7mIWpXevdIu2dRi/vdX3m6hhUTZWszjnyq7SfybtLKWN7TyzLP8rlp2nvcnpI53V10zRLh4Q2TTMvbdg4onvfS/vok1ekvRvrQG7N8j9b/5LlezR7Athgju+P0t60aOmdgjs/d8ek7a3tfRfdtaLCsarh4fR9DHuuXRzMq/vaxe+k7en+x56GznXDByf5QdM0fV2rOJh6fpTkBaWU3bvW+6S058acJL2vSx/seld1PwelaZoFaX9mDiulbNI16fVp/3g50PM3gSGkRxOo7ZWllJ6/kE9KGz62S/LxngDSNM2CUsq3kxzTeT7brUn2z8DXwS2jaW9tf3zaW/lfUtoHmk9NG8xuzbJ/fT877fWVXy6l7JXkirRfup/VaZ+RXrfO78N/llI2S9vb+Oe0vYfHpP2S2+cjIXp5ZinlsD7a72qa5uLOtZgfSPs4iEtKKd9M+9f/N6S9K2T3vv+ulPKLJB/r1HR/ktdm+d/pN6U9Fp8qpTw17V1hD0rf18n9qvPv50opP06yuGmab/SzLx9I+6Xx8lLKqWl7G96SNnC/p/9DsNJWpqbVtcrnSGd4+JvSPmLjd6WUM5L8b9rgulfa435A1/yPlVK+l/Yze1LaR/b09v8667uylDIzf3tkx4K0d08dyH929uO7aYdOP6dTf+9rEK9L+8eI95ZSxqe92colTdPc3cc6T0v7GZ9ZSnl+2mDymrTDft/ZNM1DK6hpTVqdY1XDr5LsXUp5V9phnLc1TXPValy7+J0kv0hyRillWtrP7ei05+Nxg6znqM7vkz8kubtpmkuSfDztjagu7NzE6f60PevbJDloEJcn1N7PdGpM2md/Ju2Ns16UJE3TfKRr1venvRnXZZ2bFz0tyb+lvaPzqtycDVgbhvq2t15eXuvHK30/3uQvaW8E9NZ0PfahM/+T036hejjtF54vp/2y0WT5x5ssHGC7x6T90rso7Y0upqcNBBf2mm+DtCHohs6893fm+1CScV3zNUm+0Md2Dkr7LM+70n4h/1On5smDODYD3e5/dq95j0obLBcl+WWSF6fXoyY6822bNkQsSttL8NEke2f5x5vs0JnvobQ9pqelffRD7+M8Mu2z9e5OO/y16VX/8b22/7y0N7Z5qPMZXpJk937OiV17te/Zu85+jtuga8rfHl3x5F7r6PP86RzTG1blHBmg3ucm+W7aYLCoc15+M8lL+5i357NakuRp/azvZWmfHftI2tB0fpId+jnGU7vaRqQNFfd0PpuLkjw9vR5v0pn3TWn/GPF492fSzzk3KcnpnfU+mrZH/8he80xNH48n6e886mOePpfv7/NdyWO12udIP9vePm2P7iOd9Z+5omUGsc5N0/7B4N7OZzg7vX6OBlh2i7R3AH4wyz8aadu0PYDz0/5+virJfoNc75rYz35/N/Yx74vS/gHoL2l/J3whncfqeHl5Dc9XaZq1eQkMwJrVuXbyniTfa5pmuWGQ66JSyuwkaZpmz6GtBABgcFyjCayzSikblV4PLUxyeNo7IM5e+xUBAJC4RhNYt/19ks90rve8L8kuSf457dBHN4gAABgigiawLpuT9lEa70jbi3l/krOSHNu0jyIBAGAIDOk1mqWUl6S9W9zz0z7f69VN05zXNb0kOSHt7eYnpL0I/KimaW4ZgnIBAAAYhKG+RvNJaZ+/9LZ+pr8nbU/FW9M+qPjhJD8upWy0dsoDAABgZQ2bu852nqW3tEez05t5R5KTm6b5VKdtfNpHCxzZrLlnqQEAALAahvM1mtskmZyuhwA37UPer0qye5I+g2YpZcO0Dw3v1nPtFgAAACtnkyR3NCvRSzmcg+bkzr939Wq/q2taX96X5Lg1UhEAAMAT09OS/O9gZx7OQXNVfSzJp7veb5Lkz7fffnvGjRs3RCUBAACsex588MFMmTIlSR5ameWGc9C8s/PvFknmdbVvkeS6/hZqmubRJI/2vO95lvu4ceMETQAAgLVgqO86O5Db0obNl/U0lFLGpb377JVDVRQAAAADG9IezVLK2CTP6GrappTy3CT3N00zt5Ty2SQfKKXckjZ4npj2TrTnLb82AAAAhoOhHjq7a5JLu973XFv51SRHJjkp7bM2T0syIcnlSfZpmmbRWqwRAACAlTBsnqO5pnSG2y5YsGCBazQBAFjnLV68OI899thQl8F6ZPTo0Rkxou+rKh988MGMHz8+ScY3TfPgYNc51D2aAADAIDRNkzvvvDMPPPDAUJfCembEiBHZZpttMnr06GrrFDQBAGAd0BMyJ02alDFjxix9ugKsjiVLluSOO+7IvHnzstVWW1U7rwRNAAAY5hYvXrw0ZG6++eZDXQ7rmYkTJ+aOO+7I448/ng022KDKOofz400AAIBk6TWZY8aMGeJKWB/1DJldvHhxtXUKmgAAsI4wXJY1YU2cV4ImAAAAVQmaAAAAa8mee+6Zd77znWt9u6WUnHfeeWtte4ImAAA8QSxe0uTKW+/L96/731x5631ZvKRZK9u98847c8wxx2TbbbfNhhtumClTpuSAAw7IrFmz1sr2V8eZZ56ZCRMmrPRys2fPTillucfRfO9738uJJ55Yq7xhy11nAQDgCeCiG+blhB/8PvMWLFratuX4jXLcAdOyz45brrHtzpkzJy984QszYcKEfPKTn8xOO+2Uxx57LD/+8Y/ztre9LTfddNMa2/ZwtNlmmw11CWuFHk0AAFjPXXTDvBz1tWuXCZlJcueCRTnqa9fmohvmrbFtH3300Sml5Oqrr85BBx2UZz7zmXn2s5+dd73rXfnFL36RJJk7d24OPPDAjB07NuPGjcs//dM/5a677lq6juOPPz7Pfe5zc/rpp2errbbK2LFjc/TRR2fx4sU56aSTMnny5EyaNCkf/ehHl9l2KSVf+tKX8spXvjIbb7xxtt1223znO99ZOr2vXsfrrrsupZTMmTMns2fPzhve8IYsWLAgpZSUUnL88ccnSc4+++zsuuuu2WSTTTJ58uS87nWvy913352kDdd77bVXkmTTTTdNKSVHHnlkkuWHzs6fPz+HH354Nt1004wZMyavfOUrc8sttyyd3tOj+uMf/zg77LBDxo4dm3322Sfz5v3tM/vlL3+Zl7/85Xnyk5+c8ePHZ4899si11167Oh/bahM0AQBgPbZ4SZMTfvD79DVItqfthB/8fo0Mo73//vtz0UUX5W1ve1ue9KQnLTd9woQJWbJkSQ488MDcf//9ueyyy3LxxRfnj3/8Yw455JBl5r311ltz4YUX5qKLLso555yTmTNnZr/99suf//znXHbZZfnEJz6RD3zgA7nqqquWWe6DH/xgDjrooFx//fU59NBD89rXvjY33njjoOqfPn16PvvZz2bcuHGZN29e5s2bl3e/+91J2kfOnHjiibn++utz3nnnZc6cOUvD5JQpU/Ld7343SXLzzTdn3rx5OeWUU/rcxpFHHplrrrkm559/fq688so0TZN999136SNtkuSRRx7Jpz71qZx99tn52c9+lrlz5y6tI0keeuihHHHEEbn88svzi1/8Itttt1323XffPPTQQ4PazzXB0FkAAFgHHfD5y3PPQ4+ucL5HH1+c+Y881u/0Jsm8BYuy60cuzoajRq5wfRM32TA/OOZFg6rxD3/4Q5qmybOe9ax+55k1a1Z++9vf5rbbbsuUKVOSJGeddVae/exn55e//GV22223JMmSJUty+umnZ5NNNsm0adOy11575eabb86PfvSjjBgxIttvv30+8YlP5NJLL83f/d3fLV3/wQcfnDe96U1JkhNPPDEXX3xxPv/5z+fUU09dYf2jR4/O+PHjU0rJ5MmTl5n2xje+cel/b7vttvnc5z6X3XbbLQsXLszYsWOXDpGdNGlSv9d43nLLLTn//PNzxRVXZPr06UmSr3/965kyZUrOO++8HHzwwUnaUPvlL385T3/605Mkb3/72/PhD3946Xpe+tKXLrPe0047LRMmTMhll12W/ffff4X7uSYImgAAsA6656FHc+eDi1Y84yC1YbT/QLoqmmbFvaQ33nhjpkyZsjRkJsm0adMyYcKE3HjjjUuD5tSpU7PJJpssnWeLLbbIyJEjM2LEiGXaeoav9th9992Xe3/dddet0v50+9WvfpXjjz8+119/febPn58lS5YkaYcBT5s2bVDruPHGGzNq1KhlgvHmm2+e7bfffple1zFjxiwNmUmy5ZZbLrOfd911Vz7wgQ9k9uzZufvuu7N48eI88sgjmTt37uru5ioTNAEAYB00cZMNBzXfino0e2w6ZoNB92gO1nbbbZdSSpUb/mywwQbLvC+l9NnWE/gGoyekdgfi7iGr/Xn44YczY8aMzJgxI1//+tczceLEzJ07NzNmzMhf//rXQW9/sPraz+6ajzjiiNx333055ZRTsvXWW2fDDTfM7rvvvkZqGSxBEwAA1kGDHb66eEmTF33ikty5YFGf12mWJJPHb5TL3/vSjBxRqta42WabZcaMGfniF7+Yd7zjHctdp/nAAw9khx12yO23357bb799aa/m73//+zzwwAOD7hkcyC9+8Yscfvjhy7x/3vOelySZOHFikmTevHnZdNNNk2S53s7Ro0dn8eLFy7TddNNNue+++/Lxj398ac3XXHPNcsslWW7ZbjvssEMef/zxXHXVVUuHzt533325+eabV2rfr7jiipx66qnZd999kyS333577r333kEvvya4GRAAAKzHRo4oOe6ANrT0jpE97487YFr1kNnji1/8YhYvXpwXvOAF+e53v5tbbrklN954Yz73uc9l9913z957752ddtophx56aK699tpcffXVOfzww7PHHntk1113Xe3tf/vb387pp5+e//mf/8lxxx2Xq6++Om9/+9uTJM94xjMyZcqUHH/88bnlllvywx/+MCeffPIyy0+dOjULFy7MrFmzcu+99+aRRx7JVlttldGjR+fzn/98/vjHP+b8889f7tmYW2+9dUopueCCC3LPPfdk4cKFy9W23Xbb5cADD8yb3/zmXH755bn++utz2GGH5alPfWoOPPDAQe/jdtttl7PPPjs33nhjrrrqqhx66KHZeOONV+Fo1SNoAgDAem6fHbfMlw7bJZPHb7RM++TxG+VLh+2yRp+jue222+baa6/NXnvtlX/7t3/LjjvumJe//OWZNWtWvvSlL6WUku9///vZdNNN85KXvCR77713tt1223zzm9+ssv0TTjgh3/jGN7LzzjvnrLPOyjnnnLO0t3CDDTbIOeeck5tuuik777xzPvGJT+QjH/nIMstPnz49b33rW3PIIYdk4sSJOemkkzJx4sSceeaZ+fa3v51p06bl4x//eD71qU8ts9xTn/rUnHDCCTn22GOzxRZbLA23vZ1xxhl5/vOfn/333z+77757mqbJj370o+WGyw5k5syZmT9/fnbZZZe8/vWvzzve8Y5MmjRpJY9UXWUwF+iuy0op45IsWLBgQcaNGzfU5QAAwEpbtGhRbrvttmyzzTbZaKONVrxAPxYvaXL1bffn7ocWZdImG+UF22y2xnoyh4NSSs4999y86lWvGupShrWBzq8HH3ww48ePT5LxTdM8ONh1ukYTAACeIEaOKNn96ZsPdRk8ARg6CwAAQFV6NAEAgPXS+n6Z4HCmRxMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQCAddrUqVPz2c9+drXWMXv27JRS8sADD1Spac6cOSml5LrrrquyvnWNoAkAAKxRV155ZUaOHJn99ttvqEtJkuy555555zvfuUzb9OnTM2/evIwfP36Iqlq/CJoAAMAaNXPmzBxzzDH52c9+ljvuuGOoy+nT6NGjM3ny5JRShrqU9YKgCQAArDELFy7MN7/5zRx11FHZb7/9cuaZZy6d1jNcddasWdl1110zZsyYTJ8+PTfffPPSeW699dYceOCB2WKLLTJ27Njstttu+elPf9rv9t74xjdm//33X6btsccey6RJkzJz5swceeSRueyyy3LKKaeklJJSSubMmdPn0Nkrrrgie+65Z8aMGZNNN900M2bMyPz585MkF110UV70ohdlwoQJ2XzzzbP//vvn1ltvrXTU1n2CJgAAsMZ861vfyrOe9axsv/32Oeyww3L66aenaZpl5nn/+9+fk08+Oddcc01GjRqVN77xjUunLVy4MPvuu29mzZqVX//619lnn31ywAEHZO7cuX1u701velMuuuiizJs3b2nbBRdckEceeSSHHHJITjnllOy+++5585vfnHnz5mXevHmZMmXKcuu57rrr8rKXvSzTpk3LlVdemcsvvzwHHHBAFi9enCR5+OGH8653vSvXXHNNZs2alREjRuTVr351lixZUuOwrfNGDXUBAADAKviPPZKFd6/97Y6dlLzlskHPPnPmzBx22GFJkn322ScLFizIZZddlj333HPpPB/96Eezxx4ksrUGAAAcBElEQVR7JEmOPfbY7Lffflm0aFE22mijPOc5z8lznvOcpfOeeOKJOffcc3P++efn7W9/+3Lbmz59erbffvucffbZec973pMkOeOMM3LwwQdn7NixSdphsmPGjMnkyZP7rfukk07KrrvumlNPPXVp27Of/eyl/33QQQctM//pp5+eiRMn5ve//3123HHHwR6e9ZagCQAA66KFdycPDc/rHXvcfPPNufrqq3PuuecmSUaNGpVDDjkkM2fOXCZo7rzzzkv/e8stt0yS3H333dlqq62ycOHCHH/88fnhD3+YefPm5fHHH89f/vKXfns0k7ZX87TTTst73vOe3HXXXbnwwgtzySWXrFTt1113XQ4++OB+p99yyy350Ic+lKuuuir33nvv0p7MuXPnCpoRNAEAYN00dtKw3+7MmTPz+OOP5ylPecrStqZpsuGGG+YLX/jC0rYNNthg6X/33IynJ7i9+93vzsUXX5xPfepTecYznpGNN944r3nNa/LXv/613+0efvjhOfbYY3PllVfm5z//ebbZZpu8+MUvHnTdSbLxxhsPOP2AAw7I1ltvna985St5ylOekiVLlmTHHXccsK4nEkETAADWRSsxfHUoPP744znrrLNy8skn5xWveMUy0171qlflnHPOybOe9awVrueKK67IkUcemVe/+tVJ2ms258yZM+Aym2++eV71qlfljDPOyJVXXpk3vOENy0wfPXr00mst+7Pzzjtn1qxZOeGEE5abdt999+Xmm2/OV77ylaUB9vLLL1/hvjyRCJoAAEB1F1xwQebPn59//ud/Xu7ZlAcddFBmzpyZT37ykytcz3bbbZfvfe97OeCAA1JKyQc/+MFB3XDnTW96U/bff/8sXrw4RxxxxDLTpk6dmquuuipz5szJ2LFjs9lmmy23/Pve977stNNOOfroo/PWt741o0ePzqWXXpqDDz44m222WTbffPOcdtpp2XLLLTN37twce+yxK6zpicRdZwEAgOpmzpyZvffee7mQmbRB85prrslvfvObFa7n05/+dDbddNNMnz49BxxwQGbMmJFddtllhcvtvffe2XLLLTNjxoxlhu4m7XDckSNHZtq0aZk4cWKf13s+85nPzE9+8pNcf/31ecELXpDdd9893//+9zNq1KiMGDEi3/jGN/KrX/0qO+64Y/71X/91UKH5iaT0vrXw+qaUMi7JggULFmTcuHFDXQ4AAKy0RYsW5bbbbss222yTjTbaaKjLWScsXLgwT33qU3PGGWfkH//xH4e6nGFtoPPrwQcf7PljwfimaR4c7DoNnQUAANYbS5Ysyb333puTTz45EyZMyD/8wz8MdUlPSIImAACw3pg7d2622WabPO1pT8uZZ56ZUaNEnqHgqAMAAOuNqVOnZn2/PHBd4GZAAAAAVCVoAgDAOkJPHWvCmjivBE0AABjmNthggyTJI488MsSVsD7661//miQZOXJktXW6RhMAAIa5kSNHZsKECbn77ruTJGPGjEkpZYirYn2wZMmS3HPPPRkzZkzVGycJmgAAsA6YPHlykiwNm1DLiBEjstVWW1X944WgCQAA64BSSrbccstMmjQpjz322FCXw3pk9OjRGTGi7lWVgiYAAKxDRo4cWfVaOlgT3AwIAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqGtZBs5QyspRyYinltlLKX0opt5ZSPlhKKUNdGwAAAH0bNdQFrMB7kxyV5Igkv0uya5IzkixI8rkhrAsAAIB+DPegOT3J95um+WHn/ZxSyv9J8oIhrAkAAIABDOuhs0l+nuRlpZRnJkkp5TlJXpTkwv4WKKVsWEoZ1/NKssnaKRUAAIBk+PdofjzJuCQ3lVIWJxmZ5P1N03x9gGXel+S4tVEcAAAAyxvuPZr/lOTQJK9LskvaazXfXUo5YoBlPpZkfNfraWu6SAAAAP5muPdofjLJx5um+Ubn/W9LKVun7bX8al8LNE3zaJJHe967QS0AAMDaNdx7NMckWdKrbXGGf90AAABPWMO9R/MHSd5fSpmb9vEmz0vyriSnD2lVAAAA9Gu4B81jkpyY5NQkk5LckeQ/knx4KIsCAACgf8M6aDZN81CSd3ZeAAAArANc6wgAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUN+6BZSnlqKeVrpZT7Sil/KaX8tpSy61DXBQAAQN9GDXUBAymlbJrkiiSXJnllknuSbJdk/lDWBQAAQP+GddBM8t4ktzdN84auttuGqhgAAABWbLgPnf2HJNeUUr5dSrm7lPLrUsqbh7ooAAAA+jfcg+a2SY5KckuSGUm+lORzpZQj+luglLJhKWVczyvJJmunVAAAAJLhP3R2RJJrmqb59877X5dSdkzy1iRf7WeZ9yU5bm0UBwAAwPKGe4/mvCS/79V2Y5KtBljmY0nGd72etmZKAwAAoC/DvUfziiTb92p7ZpI/9bdA0zSPJnm0530pZc1UBgAAQJ+Ge4/mZ5L8fSnl30spzyilvC7JvyT54hDXBQAAQD+GddBsmuaXSV6d5P8kuSHJB5O8s2marw9pYQAAAPRruA+dTdM0FyS5YKjrAAAAYHCGdY8mAAAA6x5BEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgqlGrslAp5UMDTW+a5sOrVg4AAADrulUKmkle3ev9Bkm2SfJ4kluTCJoAAABPUKsUNJumeV7vtlLKuCRnJjl3NWsCAABgHVbtGs2maR5MclySE2utEwAAgHVP7ZsBje+8AAAAeIJa1ZsBvaN3U5Itk7w+yYWrWxQAAADrrlW9GdC/9nq/JMk9Sb6a5GOrVREAAADrtFW9GdA2tQsBAABg/VD7Gk0AAACe4ARNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKoSNAEAAKhK0AQAAKAqQRMAAICqBE0AAACqEjQBAACoStAEAACgKkETAACAqgRNAAAAqhI0AQAAqErQBAAAoCpBEwAAgKoETQAAAKpap4JmKeXYUkpTSvnsUNcCAABA39aZoFlK2S3JW5L8ZqhrAQAAoH/rRNAspYxN8vUkb04yf4jLAQAAYADrRNBM8sUkP2ya5qcrmrGUsmEpZVzPK8kma748AAAAeowa6gJWpJTy2iS7JNltkIu8L8lxa64iAAAABjKsezRLKVOSnJLk0KZpFg1ysY8lGd/1etoaKg8AAIA+DPcezecnmZTk2lJKT9vIJC8ppbw9yYZN0yzuXqBpmkeTPNrzvms5AAAA1oLhHjRnJdmpV9sZSW5K8oneIRMAAIChN6yDZtM0DyW5obutlPJwkvuaprmh76UAAAAYSsP6Gk0AAADWPcO6R7MvTdPsOdQ1AAAA0D89mgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFXDOmiWUt5XSvllKeWhUsrdpZTzSinbD3VdAAAA9G9YB80keyT5YpK/T/LyJBsk+Ukp5UlDWhUAAAD9GjXUBQykaZp9ut+XUo5McneS5yf52VDUBAAAwMCGddDsw/jOv/f3N0MpZcMkG3Y1bbJGKwIAAGAZw33o7FKllBFJPpvkiqZpbhhg1vclWdD1+vNaKA8AAICOdSZopr1Wc8ckr13BfB9L2/PZ83raGq4LAACALuvE0NlSyheS7J/kJU3TDNhD2TTNo0ke7Vp2DVcHAABAt2EdNEubEj+f5NVJ9mya5rYhLgkAAIAVGNZBM+1w2dclOTDJQ6WUyZ32BU3T/GXoygIAAKA/w/0azaPSXmc5O8m8rtchQ1gTAAAAAxjWPZpN07jAEgAAYB0z3Hs0AQAAWMcImgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAABVCZoAAABUJWgCAABQlaAJAABAVYImAAAAVQmaAAAAVCVoAgAAUJWgCQAAQFWCJgAAAFUJmgAAAFQlaAIAAFCVoAkAAEBVgiYAAPD/t3f3oXrWdRzH35/NJmEjIt3MXGT2RDkaTkUtm7VEDSkFyQikB6hcRg0jcBKl/jOhWFYbBRG2SML6ZyIhlj1nMp06H7KsSDMx51O5Wbrl/PbHdR84O+wcH87Vfte5937BzX2ui3PD548v5/C5r9/vuqReWTQlSZIkSb2yaEqSJEmSemXRlCRJkiT1yqIpSZIkSeqVRVOSJEmS1Ks5UTSTnJ/kviRPJ9mc5LjWmSRJkiRJezf4opnkHGAdcAlwNHA7cF2SRU2DSZIkSZL2avBFE7gA+HZVXVFVdwPnAf8BPtY2liRJkiRpbwZdNJMsAJYD10+cq6pnR8cntMolSZIkSZreAa0DPIeDgfnAtinntwFv3tsHkhwIHDjp1EKA7du3/z/ySZIkSdLYerE9auhF88VYA3xp6sklS5Y0iCJJkiRJY2Eh8Lxb59CL5qPAbmDxlPOLgYem+cxaupsHTVgIPAAcDuzoO6DmPOdD03E2NBPnQzNxPjQdZ0MzGfJ8LAQefCEfGHTRrKpdSW4BVgKbAJLMGx2vn+YzO4GdE8dJJn7cUVWun9UenA9Nx9nQTJwPzcT50HScDc1k4PPxgvMMumiOrAM2JtkC3ASsBg4CrmiaSpIkSZK0V4MvmlV1VZJDgEuBQ4GtwGlVNfUGQZIkSZKkARh80QSoqvVMs1T2edgJXMKk5bTSJM6HpuNsaCbOh2bifGg6zoZmMlbzkapqnUGSJEmSNEbmtQ4gSZIkSRovFk1JkiRJUq8smpIkSZKkXlk0JUmSJEm9GuuimeT8JPcleTrJ5iTHtc6kYUjyziTXJHkwSSU5s3UmDUOSNUluTrIjycNJNiV5U+tcGoYkq5LckWT76HVjktNb59LwJLlw9P/l8tZZ1F6Si0fzMPn1x9a5NBxJXp3k+0keS/JUkjuTHNM612yMbdFMcg6wju4WwUcDtwPXJVnUNJiG4iC6mTi/dRANzgpgA3A8cArwEuAnSQ5qmkpD8QBwIbAcOAb4OXB1krc2TaVBSXIs8EngjtZZNCi/B1416fWOtnE0FEleAdwA/Bc4HXgL8Dngny1zzdbYPt4kyWbg5qr69Oh4HvB34BtVdVnTcBqUJAWcVVWbWmfR8CQ5BHgYWFFVv26dR8OT5HHg81X1ndZZ1F6SlwG3Ap8CvgBsrarVbVOptSQXA2dW1bLWWTQ8SS4D3l5VJ7XO0qexvKKZZAHdt83XT5yrqmdHxye0yiVpTnr56P3xpik0OEnmJ/kg3QqJG1vn0WBsAH5cVdc/529qf/OG0Zadvya5MslrWgfSYLwP2JLkR6NtO7cl+XjrULM1lkUTOBiYD2ybcn4bcOi+jyNpLhqthLgcuKGq7mqdR8OQZGmSJ4GdwLfoVkTc3TiWBmD0xcPRwJrWWTQ4m4GPAKcBq4AjgN8kWdgylAbjdXRz8WfgVOCbwNeTfLhpqlk6oHUASRqwDcBRuI9Ge7oHWEZ3tftsYGOSFZbN/VuSJcDXgFOq6unWeTQsVXXtpMM7Rlu8/gZ8AHDZveYBW6rqotHxbUmOAs4DNraLNTvjekXzUWA3sHjK+cXAQ/s+jqS5Jsl64AzgXVX1QOs8Go6q2lVVf6mqW6pqDd2NxT7bOpeaWw4sAm5N8kySZ+huLvaZ0fH8tvE0JFX1L+BPwOtbZ9Eg/AOY+mXlH4A5vbx6LItmVe0CbgFWTpwbLYFbiftoJM0gnfXAWcC7q+re1pk0ePOAA1uHUHM/A5bSXe2eeG0BrgSWVdXuhtk0MKObRh1JVzCkG4Cpj1J7I91V7zlrnJfOrqNbzrQFuAlYTXfDhiuaptIgjP7AT/4W8Ygky4DHq+r+RrE0DBuADwHvB3YkmdjX/URVPdUuloYgyVrgWuB+YCHdrJxMt6dG+7Gq2gHssZc7yb+Bx9zjrSRfAa6hKw6H0T1+bzfwg5a5NBhfBX6X5CLgh8BxwCdGrzlrbItmVV01eizBpXQ3ANoKnFZVU28QpP3TMcAvJh2vG71vpNusr/3XqtH7L6ec/yjw3X2aREO0CPge3TPwnqB7TuKpVfXTpqkkDd3hdKXylcAjwG+B46vqkaapNAhVdXOSs4C1wBeBe4HVVXVl22SzM7bP0ZQkSZIktTGWezQlSZIkSe1YNCVJkiRJvbJoSpIkSZJ6ZdGUJEmSJPXKoilJkiRJ6pVFU5IkSZLUK4umJEmSJKlXFk1JkiRJUq8smpIkSZKkXlk0JUmSJEm9smhKkrSPJTkkyUNJLpp07sQku5KsbJlNkqQ+pKpaZ5Akab+T5L3AJuBE4B5gK3B1VV3QNJgkST2waEqS1EiSDcB7gC3AUuDYqtrZNpUkSbNn0ZQkqZEkLwXuApYAy6vqzsaRJEnqhXs0JUlq50jgMLr/x69tG0WSpP54RVOSpAaSLABuotubeQ+wGlhaVQ83DSZJUg8smpIkNZDky8DZwNuAJ4FfAU9U1RlNg0mS1AOXzkqStI8lOZnuCua5VbW9qp4FzgVOSrKqaThJknrgFU1JkiRJUq+8oilJkiRJ6pVFU5IkSZLUK4umJEmSJKlXFk1JkiRJUq8smpIkSZKkXlk0JUmSJEm9smhKkiRJknpl0ZQkSZIk9cqiKUmSJEnqlUVTkiRJktQri6YkSZIkqVcWTUmSJElSr/4Hy72uPkOCy5MAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "#Imports for animation and display within a jupyter notebook\n", "from matplotlib import animation, rc \n", "from IPython.display import HTML\n", "\n", "#Generating the figure that will contain the animation\n", "fig, ax = plt.subplots()\n", "fig.set_dpi(100)\n", "fig.set_size_inches(9, 5)\n", "ax.set_xlim(( 0, 2*np.pi))\n", "ax.set_ylim((0, 10))\n", "comp, = ax.plot([], [], marker='o', lw=2,label='Computational')\n", "anal, = ax.plot([], [], lw=2,label='Analytical')\n", "ax.legend();\n", "plt.xlabel('x')\n", "plt.ylabel('u')\n", "plt.title('Burgers Equation time evolution from t=0 to t=10');\n", "\n", "#Resetting the U wave back to initial conditions\n", "u = np.asarray([ufunc(0, x0, nu) for x0 in x])" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "#Initialization function for funcanimation\n", "def init():\n", " comp.set_data([], [])\n", " anal.set_data([], [])\n", " return (comp,anal,)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "#Main animation function, each frame represents a time step in our calculation\n", "def animate(j):\n", " un = u.copy() #copy the u array to not overwrite values\n", " for i in range(1,grid_points-1):\n", " u[i] = un[i] - un[i] * dt/dx * (un[i]-un[i-1]) + nu * (dt/dx**2) * (un[i+1]- 2*un[i] + un[i-1]) \n", " \n", " u[0] = un[0] - un[0] * dt / dx * (un[0] - un[-2]) + nu*(dt / dx**2) *(un[1] - 2* un[0] + un[-2])\n", " u[-1] = u[0]\n", " u_anal = np.asarray([ufunc(j * dt, xi, nu) for xi in x])\n", " comp.set_data(x, u)\n", " anal.set_data(x, u_anal)\n", " return (comp,anal,)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "anim = animation.FuncAnimation(fig, animate, init_func=init,\n", " frames=nt, interval=20)\n", "anim.save('../gifs/1dBurgers.gif',writer='imagemagick',fps=60)\n", "#HTML(anim.to_jshtml())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Conclusion\n", "\n", "This concludes our examination of 1D sims and boy oh boy was this cool! This last model in particular shines in the animation showing the behavior and properties of the burghers equation quite well. \n", "\n", "Next, we will start our move to 2D but before this a quick detour on array operations on NumPy." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.5" } }, "nbformat": 4, "nbformat_minor": 2 }