{ "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": "iVBORw0KGgoAAAANSUhEUgAAAv8AAAHPCAYAAADEabUSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAIABJREFUeJzt3XuclWW9///XBxAUOWogpCiQeMBDpugOssTU8ETqVrOdJ2xreci22/yWlhVm/VLTUvPQto2S5jazPGUeMhTbmoJoWCYSqYRuR/GAA6igDNfvj3utac1iDmuGGdYa7tfz8VgPWNd9WJ/7MDPvdd/Xfd+RUkKSJEnS+q9HtQuQJEmStG4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlhOFfkiRJygnDvyRJkpQThn9JkiQpJwz/kiRJUk4Y/iXlSkSkiJha7TpK1WJN1RARUwrrYmQnz3dmRMzszHm247PHRMTvIqK+sGyHVqMOSSoy/EvNKAkhpa/FEfFgRBxQ7fqqqZn1Uvr6SbXrA4iIA2stTNdiTeuTiBgbEVM7+4tDJ/gZsBPwDeBYYE51y2ldRHy9s7+gRMSgiLgmIl6LiLcLv0d3rXDaUyNiSmfWU5hvpy5nRAyPiAsKy7as8PtwYivjT4iIhyPinYh4JSIuj4h+nVWP1Jpe1S5AqnHfAl4AAtgMmALcHRGTU0p3VbOwKrsfuL6Z9r+t60JacCBwGjC1mWEbAavWaTWZWqxpfTIW+DYwE1hYNuxT67oYgIjYCBgPfC+ldEU1auiArwO/Am7vjJlFRA/gt8CHgR8ArwOnAjMjYreU0oI2ZnFqYZrpnVFPiU5dTmBb4GvAAuAvZNu9WRGxCzADmAecCWwBnAWMAXJ9cEnrhuFfat09KaXGI3URMQ14Ffg3oFPCf0RsnFJ6uzPmVcFn9QB6p5RWrOWs/pZS+nln1LSudcKyd7parGl9klJ6r0ofPaTw71ttjbgufw+sY0cAE4AjU0q/AoiIX5IdKDgP+FwVa+tMTwCbppTejIgjgFtaGff/A5YAE1NKSwEiYiHw04j4VErpd11erXLNbj9S+7wFvEvJUdqImNjcKd6IGFlon1LSNj0ilkfEhyLi7ohYBtxYMvy0iHg+It6NiNkR8fHm+itHRJ+IOC8i/h4RKyPixYi4KCL6lI2XIuKKiDg6Iv4KrAT2Lwz7bEQ8UThFvTQi/hIR/9FZK6rwGV+IiOdaW56W+nk3t14L098SEYtKlvtHhSOsxXGmkx1hb9JFqWydTC37rI9ExD2F9bA8ImZExEfLxinW+bGI+GFJF4bbImIIrWhvTYXuKykitomIn0fWX/y1iDg/MiMi4o5Cva9ExFea+cyK9pFWav6XiLi38NnvRMRDEfGxkuFHFGrcq5lpv1gYtmNJ2ycj4n8L6+ytQv3bV1BHs9dDRMTCwnql8DNWDFsPlqzjiYXhzf0MDY2IaRHxakSsiIinIuL4snGKP8NnlezLKyPi8YjYvY26pwL/KLz9QWE+C4vDCu/HRsT/RMQS4OH2rKvO2EeaqTkBGwPHl6zD6W1N14YjyA6Y3FpsSCm9BvwSOKS1/bGwvnYA9iqpZ2bJ8NGR/T54s7CPPhYRB7VVUFcsZ0ppWUrpzQo+ewCwH/DzYvAvuB5YDnxmbeqQKuGRf6l1AyPiA2TdfoYCpwP9gLU56t0LuI/sj/1ZwDsAEXEKcAXwv8CPgJFkp6SXAC8VJ47s6P2dwJ7ANWSnjncC/hPYBijvx/pJsj8oV5CdPl8YEfsBN5Gdev5aYbztgY8Bl1WwDBsW1ku5pcWjrBHx78B/AX8ELgVGF+p+E3ixgs9ozpFAX+Bq4A1gD7JtskVhGIXP/CDZH9hj25phROxAts6XAhcB7wNfJOuWsFdKaVbZJD8m2ybnkW2jM8jW7VGtfEy7aipxM9n2PRs4CDiXbP19EXiAbNsdDVwcEY+nlP5QWKb27iNNRMQngXvIjmaeB6wGTgAeiIiPp5Rmk3XlKIaVh8pmcRTw15TS04X57VuY3/Nk3Z42Ittuj0TErimlhe1YJ835A3A58GWyo6rzCu3zmhs5si+LM4GtybbdC2T7z/SIGJRSKv8Z+BzQn2w7JuCrwK0RMTql9H4LNd1KdrDgR2Q/a3eTra9St5B1E/k62e+YjqyrDu0jLTgW+G9gNtl+A/Bcoa4NgIGtTFvqzZTS6sL/PwI8WfK+aDbwBbL98S8tzOcMsp+35cD3Cm2vFurZjOx3S1+ybf8GcDxwZ0QckVK6rZX6umI5K7UT2d+AJtd+pJTei4i5ZOtL6lopJV++fJW9yPr2p2ZeK4Djy8adWBg2sax9ZKF9Sknb9ELb98vG7U0WzGcDvUrajy+MP7Ok7RigAdizbB5fLIw7oaQtFcYdWzbupUA90LMD66a59VJ8fbYwzgZkf6T/RNbNqDjtSc0sT3Fdj2xrvQIbNVPP2WThdMuStiuyX28t1j+15P1tZGdERpe0DSf7MvBQM3XeD0RJ+w/JzgQNbGO9taemqYW2/ypp60n2pWk18LWS9kFkXyCnd2QfaaaWIOuScW/Zcm5EFkh/V9L2P4Xt3LOkbVjhs79Z0vanwniblLTtXBjvZ63tC+XrpqR9YdkyH1G+v5QMm1m2z/1HYdyjS9o2IAuTy4D+ZT/DrwODS8b9dKH94Da2eXH6s8rai9v3f5qZptJ1tVb7SCs1L29uPP7581jJa2TZ/KY1M78DC+NOaqOep0u3XUn7jwrT71nS1q+wj74A9FiXy1k2j9b2xeKwjzcz7JdAXVvbyJevtX155F9q3Wn88yLWzchC1X9HxLKU0q0tT9amq8vejwM2Bc5JKZVe+Hkj2R+5UkeSHel7tuzo+wOFf/cmCzFFD6WUnimbx1tkp733Iwt57XUHWZgtVzyCN47sTMm3UtP+1tPJLvrrkJTSu8X/R8TGZIH0j2SB9SPAovbMLyJ6kl0MentK6fmSz6mLiP8BToqIAanp6flrUkqp5P3/kh1R3wr4czsXqS3/XVJTQ0TMITvLMa2k/a2ImE92ZqWovftIqV3ILjz8LrBpRJQOmwEcGxE9UnbE82ay618mFoZBFm56FIYREcML87wolXSLSCn9OSLuJwuB69qBwCtkR+SL9bwfEZcX2vai6TU9N6eUlpS8/9/Cv6XrvCOa3B2rg+uqo/tIez1F9vuiEq+U/H8jsi/X5VaUDO+IA4HZKaXG7lIppeURcQ3wfbILwJ/uwHw7upyVKi5vS+uko+tDqpjhX2rd7NT0gt+byI7MXRERd6WOXUi4ipJuPAVbFf79e2ljSmlVsZ9wiTFkXXRea2H+Q8vev9DMOFeRdde4JyL+D/gd8MuUUqVfBF5KKf2+leHF5WlyJ49CwHq+mfErEhFbAt8hO/I6uGxwpafqSw0h6zYwv5lh88hC7AjgryXt5V8wiqGwvJ7OUP5Z9cCKlNLrzbRvWvK+vftIqTGFf3/WyjgDyZb73sJnH8U/w/9RwNyUUvFLc3FfaGkdT4p1f7HrVsCCtGaXjXklw0s12Q4ppSWFL0Vru83LfzY7sq46uo+0S+HLT2s/8y15F2iuX/+GJcM7YiugvEseNN2G7Q7/a7GclSoub0vrpKPrQ6qY4V9qh5TS6oh4kKzbwBiyUJhaGL1nC+0rmwkd7dGD7Aj7mS0ML+9Pv8Yfk5TS4shuNzeJ7NZyBwAnRMT1KaXj16K2jqho/RWO0t8PbAJcCDwLvA1sTnZGYV3dwKChhfZoob2zP6uSz2/vPlKquB7/HzC3hXGWA6SUVkbE7cBhEXEq2dmxj5H1Ye9qLf18dYWu2uadEfQ6uo+0S0T0JvvZq8RrKaViDXVk3ejKFdte7mhNXWEtlrNSdYV/W1onNbU+tH4y/EvtV/y5KT6QpXjkd1DZeOVHD1tTvCvI1sCDxcaI6EXWb7i0O8lzZPfMnlHW/aRdCmctfgP8pnCB6FXAFyPi/JTS31ufuk3F5RnDP7uaFC+mG0V2ar2o0vW3E9nFgcenlBqfMVC4eLlcpevlNbK+0Ns2M2w7sr7THb04uaM1dYa12UeeK/y7tI2zO0U3k12bsg/Z2YYotBUV94WW1vHrbRz1X0LZvlEIaOXhqT3L+Q9g55LuS6X1FIdXw9quq87Q0nqcQMnvpjaM4p/PWpgLfLyZdf0vZD97bT0bpKV6/kHL66k4vCPz7ehyVuppsrO/48j6+AON+/QupW1SV/FWn1I7FMLrp4D3+Ofp5X+QHWn7RNnop7Zj1nPI7lZxUiHwFx3Nml0Lfkl2tPukZurbqNAXvlUR0eT0f+GPcvELRkW3gmzDHLJgfXLhj1rRFNYM+cWw2bj+Ckf5v1A2XvEIW5SMF2RnYcq9XRhe/llNFI7a/Y7sloMjS+a7GdkdXh4u6++/NiqqqZOszT7yBNk2OSuaeeJorHlb09+T3V3mqMJrdkqpsTtLSqmOLAAeX7rskd0G9FNkd8FpzXOs+bP1BdY88l8MxZWs37vJLkxuvENT4efudLKzGuV3L1onOmFddYa3aX4dFvvCV/Iq7Qv/K7IzQv9abChch3Ik8JuUUnN93yup525gj4hofJhWYb/+AlkgL7/OqdL5dnQ5K5JSqif7mTkmIvqXDDqW7IBSa88HkDqFR/6l1h0QEcUjSUPJAuEY4IJiKEwp1UfELcDphftHPwccTOv9qptI2W3eppLd1u6ByB6CM5IsLD9H06NUN5D11/9JROwNPEIWhLYrtE+i7DZyzfjviNiE7Kj8S2RH2U8nCx7N3h6xzDYRcUwz7a+mlO4v9O0/l+zWiA9ExM1kR8lOILsbR+my/zUiHgO+X6jpTeCzrPn76VmydXFxRGxOdjeew2m+3/UThX8vj4j7gIaU0i9aWJZzyf6QPxwRV5Edlfsi2Zegr7a8CtqtPTWtrQ7vI4WubSeS3W7yrxFxHfB/ZF8m9iZb75NLxn8/Im4l22Ybk92+ttz/K8zv0cgelFe8fWU9zT/xuNR/F5bj12Tdvj5cqL+8T/tcsi+IX4uIgWQXVD6QUlrczDyvIdvG0yNiN7KweARZl6UzUkrL2qipK63NuuoMTwD7RsSZZF1QXkgpzVqLvvC/Ah4DrouIsfzzCb89yZ7IXEk9pxR+n/wdWJxSegC4gOxi83sKF2q/SXYGahRweAVdKzt7OSnUCNmzCSC7OH5PgJTSd0tG/QbZBfcPFS5Q3gL4CtmdtDpyAwapfap9uyFfvmrxRfO3+nyX7GLfkym5BWJh/A+Q/ZF7m+yP0E/I/gAk1rzV5/JWPvd0siCyguxitglkIe2esvE2IAumTxfGfbMw3reAASXjJeCKZj7ncLJnDbxKFpL+Uah5WAXrprVb380sG/cUsrC/Angc+Dhlt10sjDeaLNitIDua9j1gX9a81ef2hfGWkZ1ZuIbsNojl67kn2b2/F5N13Ull9U8t+/yPkF28uqywDR8AxrewT4wra59YXmcL663imvjnbRw/UDaPZvefwjp9uiP7SCv17gL8miysrSjslzcDn2xm3OK2Wg1s0cL89iF7tsU7ZEH2TmD7FtbxyJK2HmRB77XCtrkX+BBlt/osjHsi2RfEVaXbpIV9bihwbWG+K8nOfE0pG2ckzdyqs6X9qJlxmp2+pe3bznW11vtIC5+9LdmZj3cK85/e1jQVzHMw2Ze41wvbcCZlP0etTLsZ2Z2XlrLmbYJHkx0pX0L2+3kWcFCF8+2K5Wzxd2Mz4+5J9qX8XbLfCVdQuMWsL19d/YqU1mU3VEntUeiL/xpwa0ppjS4c3VEUntCZUppY3UokScof+/xLNSIiNoyym6oDx5HdeWLmuq9IkiStb+zzL9WOjwI/Klw/8AawK/DvZN02vAhMkiStNcO/VDsWkt1W8stkR/vfBK4Hzk4de5iYJElSE1Xt8x8RnyC7s8FuZPdsPiyldHvJ8ADOI7td3SCyi2NOSSktaGZ2kiRJklpR7T7/G5PdU/e0FoZ/lewo6MlkDwR5G7gvIjZsYXxJkiRJLaiZu/0U7o/eeOS/cNT/ZeCSlNLFhbaBZLcmnJK67v7YkiRJ0nqplvv8jyJ7AmPjwzZS9jClWcB4oNnwHxF9WPMJpcX+05IkSdL6oj/wcmrH0fxaDv/DCv++Wtb+asmw5pxDZU8NlCRJkrq7LciexF6RWg7/HfV94Icl7/sDL7344osMGDCgSiVJkiRJnWfp0qWMGDECsqfTV6yWw/8rhX83A+pK2jcD5rY0UUppJdmj2gEoPjNpwIABhn9JkiTlWrXv9tOaF8i+AOxTbIiIAWR3/Xm0WkVJkiRJ3VVVj/xHRD9g65KmURGxC/BmSmlRRFwKnBsRC8i+DJxPdgeg29ecmyRJkqTWVLvbzzjgwZL3xb76PwOmABeRPQvgGrKHfD0M7J9SWrEOa5QkSZLWCzVzn/+uUugqVF9fX2+ff0mS1GUaGhp4//33q12G1iO9e/emR4/me+kvXbqUgQMHAgxMKS2tdJ7VPvIvSZLUraWUeOWVV3jrrbeqXYrWMz169GDUqFH07t270+Zp+JckSVoLxeA/dOhQ+vbt23inQWltrF69mpdffpm6ujq23HLLTtuvDP+SJEkd1NDQ0Bj8N91002qXo/XMkCFDePnll1m1ahUbbLBBp8yzlm/1KUmSVNOKffz79u1b5Uq0Pip292loaOi0eRr+JUmS1pJdfdQVumK/MvxLkiRJOWH4lyRJUu5NnDiRM844Y51/bkRw++3r7vm1hn9JkqQqa1idePS5N7hj7v/x6HNv0LB63TyH6ZVXXuH0009n9OjR9OnThxEjRjB58mRmzJixTj5/bUyfPp1Bgwa1e7qZM2cSEWvcmvXWW2/l/PPP76zyapZ3+5EkSaqie5+u47zfPENd/YrGtuEDN+Tbk8ey/47Du+xzFy5cyMc+9jEGDRrED37wA3baaSfef/997rvvPk477TSeffbZLvvsWrTJJptUu4R1wiP/kiRJVXLv03Wc8vMnmwR/gFfqV3DKz5/k3qfruuyzTz31VCKC2bNnc/jhh7PNNtuwww47cOaZZ/LYY48BsGjRIg455BD69evHgAED+MxnPsOrr77aOI+pU6eyyy67cO2117LlllvSr18/Tj31VBoaGrjooosYNmwYQ4cO5Xvf+16Tz44Irr76ag444AA22mgjRo8eza9+9avG4c0dnZ87dy4RwcKFC5k5cyYnnHAC9fX1RAQRwdSpUwG44YYbGDduHP3792fYsGF87nOfY/HixUD2hWfvvfcGYPDgwUQEU6ZMAdbs9rNkyRKOO+44Bg8eTN++fTnggANYsGBB4/DimYf77ruP7bffnn79+rH//vtTV/fPbfb444+z33778YEPfICBAwey11578eSTT67NZltrhn9JkqQqaFidOO83z9BcB59i23m/eaZLugC9+eab3HvvvZx22mlsvPHGawwfNGgQq1ev5pBDDuHNN9/koYce4v777+f555/nqKOOajLuc889xz333MO9997LTTfdxLRp0zjooIN46aWXeOihh7jwwgs599xzmTVrVpPpvvnNb3L44Yfz1FNPcfTRR/PZz36WefPmVVT/hAkTuPTSSxkwYAB1dXXU1dVx1llnAdntV88//3yeeuopbr/9dhYuXNgY8EeMGMGvf/1rAObPn09dXR2XXXZZs58xZcoU5syZw5133smjjz5KSokDDzyw8fauAO+88w4XX3wxN9xwA3/4wx9YtGhRYx0Ay5Yt4/jjj+fhhx/mscceY8yYMRx44IEsW7asouXsCnb7kSRJ6kSTf/wwry1b2eZ4K1c1sOSd91scnoC6+hWM++799OnVs835Denfh9+cvmdFNf79738npcR2223X4jgzZszgL3/5Cy+88AIjRowA4Prrr2eHHXbg8ccfZ/fddweyJ9Fee+219O/fn7Fjx7L33nszf/587r77bnr06MG2227LhRdeyIMPPsi//Mu/NM7/yCOP5MQTTwTg/PPP5/777+fHP/4xV111VZv19+7dm4EDBxIRDBs2rMmwz3/+843/Hz16NJdffjm77747y5cvp1+/fo3de4YOHdriNQMLFizgzjvv5JFHHmHChAkA3HjjjYwYMYLbb7+dI488Esi+aPzkJz/hQx/6EABf+tKX+M53vtM4n09+8pNN5nvNNdcwaNAgHnroIQ4++OA2l7MrGP4lSZI60WvLVvLK0hVtj1ih7AtCy18SOiKlts8mzJs3jxEjRjQGf4CxY8cyaNAg5s2b1xj+R44cSf/+/RvH2WyzzejZsyc9evRo0lbselM0fvz4Nd7PnTu3Q8tT6oknnmDq1Kk89dRTLFmyhNWrVwNZF6axY8dWNI958+bRq1evJl9WNt10U7bddtsmZyf69u3bGPwBhg8f3mQ5X331Vc4991xmzpzJ4sWLaWho4J133mHRokVru5gdZviXJEnqREP696lovLaO/BcN7rtBxUf+KzVmzBgiolMu6t1ggw2avI+IZtuKIbwSxS8OpV9SSrvbtOTtt99m0qRJTJo0iRtvvJEhQ4awaNEiJk2axHvvvVfx51equeUsrfn444/njTfe4LLLLmOrrbaiT58+jB8/vktqqZThX5IkqRNV2vWmYXVizwsf4JX6Fc32+w9g2MANefhrn6Rnj8590usmm2zCpEmTuPLKK/nyl7+8Rr//t956i+23354XX3yRF198sfHo/zPPPMNbb71V8RH01jz22GMcd9xxTd5/5CMfAWDIkCEA1NXVMXjwYIA1zgr07t2bhoaGJm3PPvssb7zxBhdccEFjzXPmzFljOmCNaUttv/32rFq1ilmzZjV2+3njjTeYP39+u5b9kUce4aqrruLAAw8E4MUXX+T111+vePqu4AW/kiRJVdCzR/DtyVmQLI/2xfffnjy204N/0ZVXXklDQwN77LEHv/71r1mwYAHz5s3j8ssvZ/z48ey7777stNNOHH300Tz55JPMnj2b4447jr322otx48at9effcsstXHvttfztb3/j29/+NrNnz+ZLX/oSAFtvvTUjRoxg6tSpLFiwgN/+9rdccsklTaYfOXIky5cvZ8aMGbz++uu88847bLnllvTu3Zsf//jHPP/889x5551r3Lt/q622IiK46667eO2111i+fPkatY0ZM4ZDDjmEk046iYcffpinnnqKY445hs0335xDDjmk4mUcM2YMN9xwA/PmzWPWrFkcffTRbLTRRh1YW53H8C9JklQl++84nKuP2ZVhAzds0j5s4IZcfcyuXXqf/9GjR/Pkk0+y995785WvfIUdd9yR/fbbjxkzZnD11VcTEdxxxx0MHjyYT3ziE+y7776MHj2am2++uVM+/7zzzuMXv/gFO++8M9dffz033XRT41H1DTbYgJtuuolnn32WnXfemQsvvJDvfve7TaafMGECJ598MkcddRRDhgzhoosuYsiQIUyfPp1bbrmFsWPHcsEFF3DxxRc3mW7zzTfnvPPO4+yzz2azzTZr/MJR7rrrrmO33Xbj4IMPZvz48aSUuPvuu9fo6tOaadOmsWTJEnbddVeOPfZYvvzlLzN06NB2rqnOFZVc8NGdRcQAoL6+vp4BAwZUuxxJkrQeWbFiBS+88AKjRo1iww03bHuCFjSsTsx+4U0WL1vB0P4bsseoTbrsiH8tiAhuu+02Dj300GqXUtNa27+WLl3KwIEDAQamlJZWOk/7/EuSJFVZzx7B+A9tWu0ylAN2+5EkSZJywiP/kiRJWqfW927ntcwj/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZK6xMiRI7n00kvXah4zZ84kInjrrbc6paaFCxcSEcydO7dT5tfdGP4lSZJy6tFHH6Vnz54cdNBB1S4FgIkTJ3LGGWc0aZswYQJ1dXUMHDiwSlWtXwz/kiRJOTVt2jROP/10/vCHP/Dyyy9Xu5xm9e7dm2HDhhER1S5lvWD4lyRJyqHly5dz8803c8opp3DQQQcxffr0xmHFrjYzZsxg3Lhx9O3blwkTJjB//vzGcZ577jkOOeQQNttsM/r168fuu+/O73//+xY/7/Of/zwHH3xwk7b333+foUOHMm3aNKZMmcJDDz3EZZddRkQQESxcuLDZbj+PPPIIEydOpG/fvgwePJhJkyaxZMkSAO6991723HNPBg0axKabbsrBBx/Mc88910lrrfsz/EuSJOXQL3/5S7bbbju23XZbjjnmGK699lpSSk3G+cY3vsEll1zCnDlz6NWrF5///Ocbhy1fvpwDDzyQGTNm8Kc//Yn999+fyZMns2jRomY/78QTT+Tee++lrq6use2uu+7inXfe4aijjuKyyy5j/PjxnHTSSdTV1VFXV8eIESPWmM/cuXPZZ599GDt2LI8++igPP/wwkydPpqGhAYC3336bM888kzlz5jBjxgx69OjBYYcdxurVqztjtXV7vapdgCRJ0nrlv/aC5YvX/ef2GwpffKji0adNm8YxxxwDwP777099fT0PPfQQEydObBzne9/7HnvttRcAZ599NgcddBArVqxgww035MMf/jAf/vCHG8c9//zzue2227jzzjv50pe+tMbnTZgwgW233ZYbbriBr371qwBcd911HHnkkfTr1w/Iuvj07duXYcOGtVj3RRddxLhx47jqqqsa23bYYYfG/x9++OFNxr/22msZMmQIzzzzDDvuuGOlq2e9ZfiXJEnqTMsXw7La7D9fNH/+fGbPns1tt90GQK9evTjqqKOYNm1ak/C/8847N/5/+PDhACxevJgtt9yS5cuXM3XqVH77299SV1fHqlWrePfdd1s88g/Z0f9rrrmGr371q7z66qvcc889PPDAA+2qfe7cuRx55JEtDl+wYAHf+ta3mDVrFq+//nrjEf9FixYZ/jH8S5Ikda5+Q2v+c6dNm8aqVav44Ac/2NiWUqJPnz5cccUVjW0bbLBB4/+LF9wWw/RZZ53F/fffz8UXX8zWW2/NRhttxBFHHMF7773X4uced9xxnH322Tz66KP88Y9/ZNSoUXz84x+vuG6AjTbaqNXhkydPZquttuKnP/0pH/zgB1m9ejU77rhjq3XlieFfkiSCDu+HAAAXEUlEQVSpM7Wj6001rFq1iuuvv55LLrmET33qU02GHXroodx0001st912bc7nkUceYcqUKRx22GFAdg3AwoULW51m00035dBDD+W6667j0Ucf5YQTTmgyvHfv3o1991uy8847M2PGDM4777w1hr3xxhvMnz+fn/70p41fKh5++OE2lyVPDP+SJEk5ctddd7FkyRL+/d//fY175x9++OFMmzaNH/zgB23OZ8yYMdx6661MnjyZiOCb3/xmRRfVnnjiiRx88ME0NDRw/PHHNxk2cuRIZs2axcKFC+nXrx+bbLLJGtOfc8457LTTTpx66qmcfPLJ9O7dmwcffJAjjzySTTbZhE033ZRrrrmG4cOHs2jRIs4+++w2a8oT7/YjSZKUI9OmTWPfffdt9qFZhx9+OHPmzOHPf/5zm/P54Q9/yODBg5kwYQKTJ09m0qRJ7Lrrrm1Ot++++zJ8+HAmTZrUpNsRZF2JevbsydixYxkyZEiz1w9ss802/O53v+Opp55ijz32YPz48dxxxx306tWLHj168Itf/IInnniCHXfckf/8z/+s6ItMnkT5LZ3WNxExAKivr69nwIAB1S5HkiStR1asWMELL7zAqFGj2HDDDatdTrewfPlyNt98c6677jr+9V//tdrl1LTW9q+lS5cWv8ANTCktrXSedvuRJElSl1u9ejWvv/46l1xyCYMGDeLTn/50tUvKJcO/JEmSutyiRYsYNWoUW2yxBdOnT6dXL2NoNbjWJUmS1OVGjhy5xhOEte55wa8kSZKUE4Z/SZKkteQRbXWFrtivDP+SJEkdVHwC7jvvvFPlSrQ+Kj6VuGfPnp02T/v8S5IkdVDPnj0ZNGgQixcvBqBv375ERJWr0vpg9erVvPbaa/Tt27dTL442/EuSJK2FYcOGATR+AZA6S48ePdhyyy079Qul4V+SJGktRATDhw9n6NChvP/++9UuR+uR3r1706NH5/bSN/xLkiR1gp49e3Zq32ypK3jBryRJkpQThn9JkiQpJwz/kiRJUk4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlhOFfkiRJygnDvyRJkpQThn9JkiQpJwz/kiRJUk4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlRE2H/4joGRHnR8QLEfFuRDwXEd+MiKh2bZIkSVJ306vaBbTha8ApwPHAX4FxwHVAPXB5FeuSJEmSup1aD/8TgDtSSr8tvF8YEf8G7FHFmiRJkqRuqaa7/QB/BPaJiG0AIuLDwJ7APVWtSpIkSeqGav3I/wXAAODZiGgAegLfSCnd2NIEEdEH6FPS1L9rS5QkSZK6h1o/8v8Z4Gjgc8CuZH3/z4qI41uZ5hyyawKKr5e6ukhJkiSpO4iUUrVraFFEvAhckFK6sqTtXOCYlNJ2LUzT3JH/l+rr6xkwYECX1itJkiStC0uXLmXgwIEAA1NKSyudrta7/fQFVpe1NdDKGYuU0kpgZfG9dwWVJEmSMrUe/n8DfCMiFpHd6vMjwJnAtVWtSpIkSeqGaj38nw6cD1wFDAVeBv4L+E41i5IkSZK6o5oO/ymlZcAZhZckSZKktVDrd/uRJEmS1EkM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKSdqPvxHxOYR8fOIeCMi3o2Iv0TEuGrXJUmSJHU3vapdQGsiYjDwCPAgcADwGjAGWFLNuiRJkqTuqKbDP/A14MWU0gklbS9UqxhJkiSpO6v1bj+fBuZExC0RsTgi/hQRJ7U2QUT0iYgBxRfQf92UKkmSJNW2Wg//o4FTgAXAJOBq4PKIOL6Vac4B6kteL3V1kZIkSVJ3ECmlatfQooh4D5iTUppQ0nY5sHtKaXwL0/QB+pQ09Qdeqq+vZ8CAAV1aryRJkrQuLF26lIEDBwIMTCktrXS6Wu/zXwc8U9Y2Dzi8pQlSSiuBlcX3EdE1lUmSJEndTK13+3kE2LasbRvgH1WoRZIkSerWaj38/wj4aER8PSK2jojPAV8ArqxyXZIkSVK3U9PhP6X0OHAY8G/A08A3gTNSSjdWtTBJkiSpG6r1Pv+klO4C7qp2HZIkSVJ3V9NH/iVJkiR1HsO/JEmSlBOGf0mSJCknDP+SJElSThj+JUmSpJww/EuSJEk5YfiXJEmScsLwL0mSJOWE4V+SJEnKCcO/JEmSlBOGf0mSJCknDP+SJElSThj+JUmSpJww/EuSJEk5YfiXJEmScsLwL0mSJOWE4V+SJEnKCcO/JEmSlBOGf0mSJCknDP+SJElSThj+JUmSpJww/EuSJEk5YfiXJEmScsLwL0mSJOWE4V+SJEnKCcO/JEmSlBOGf0mSJCknDP+SJElSThj+JUmSpJzo1ZGJIuJbrQ1PKX2nY+VIkiRJ6iodCv/AYWXvNwBGAauA5wDDvyRJklRjOhT+U0ofKW+LiAHAdOC2taxJkiRJUhfotD7/KaWlwLeB8ztrnpIkSZI6T2df8Duw8JIkSZJUYzp6we+Xy5uA4cCxwD1rW5QkSZKkztfRC37/s+z9auA14GfA99eqIkmSJEldoqMX/I7q7EIkSZIkdS0f8iVJkiTlhOFfkiRJygnDvyRJkpQThn9JkiQpJwz/kiRJUk4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlhOFfkiRJygnDvyRJkpQThn9JkiQpJwz/kiRJUk4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlhOFfkiRJygnDvyRJkpQThn9JkiQpJwz/kiRJUk4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlhOFfkiRJygnDvyRJkpQT3Sr8R8TZEZEi4tJq1yJJkiR1N90m/EfE7sAXgT9XuxZJkiSpO+oW4T8i+gE3AicBS6pcjiRJktQtdYvwD1wJ/Dal9PtqFyJJkiR1V72qXUBbIuKzwK7A7hWO3wfoU9LUvyvqkiRJkrqbmj7yHxEjgMuAo1NKKyqc7BygvuT1UheVJ0mSJHUrkVKqdg0tiohDgduAhpLmnkACVgN9UkoNZdM0d+T/pfr6egYMGNDFFUuSJEldb+nSpQwcOBBgYEppaaXT1Xq3nxnATmVt1wHPAheWB3+AlNJKYGXxfUR0aYGSJElSd1HT4T+ltAx4urQtIt4G3kgpPd38VJIkSZKaU9N9/iVJkiR1npo+8t+clNLEatcgSZIkdUce+ZckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJyoqbDf0ScExGPR8SyiFgcEbdHxLbVrkuSJEnqjmo6/AN7AVcCHwX2AzYAfhcRG1e1KkmSJKkb6lXtAlqTUtq/9H1ETAEWA7sBf6hGTZIkSVJ3VdPhvxkDC/++2dIIEdEH6FPS1L9LK5IkSZK6iVrv9tMoInoAlwKPpJSebmXUc4D6ktdL66A8SZIkqeZ1m/BP1vd/R+CzbYz3fbIzBMXXFl1clyRJktQtdItuPxFxBXAw8ImUUqtH8lNKK4GVJdN2cXWSJElS91DT4T+y5P5j4DBgYkrphSqXJEmSJHVbNR3+ybr6fA44BFgWEcMK7fUppXerV5YkSZLU/dR6n/9TyPrtzwTqSl5HVbEmSZIkqVuq6SP/KSU77EuSJEmdpNaP/EuSJEnqJIZ/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScoJw78kSZKUE4Z/SZIkKScM/5IkSVJOGP4lSZKknDD8S5IkSTlh+JckSZJywvAvSZIk5YThX5IkScqJbhH+I+K0iFgYESsiYlZE7FHtmiRJkqTupubDf0QcBfwQOA/YFXgKuC8ihla1MEmSJKmbqfnwD5wJ/DSldF1K6RngZOAd4PPVLUuSJEnqXmo6/EdEb2A34PfFtpTS6sL78dWqS5IkSeqOelW7gDZ8AOgJvFrW/iqwXXMTREQfoE9JU3+ApUuXdkV9kiRJ0jrX0Wxb6+G/I84Bvl3eOGLEiCqUIkmSJHWp/kDF3wRqPfy/DjQAm5W1bwa80sI03ye7QLioP/ASsAWwrLMLVJdy23Vfbrvuy23Xfbntui+3XfdV7W3XH3i5PRPUdPhPKb0XEU8A+wC3A0REj8L7K1qYZiWwsvg+Ior/XZZSsu9PN+K2677cdt2X2677ctt1X2677qsGtl27P7Omw3/BD4GfRcQcYDZwBrAxcF1Vq5IkSZK6mZoP/ymlmyNiCPAdYBgwF9g/pVR+EbAkSZKkVtR8+AdIKV1BC918KrCS7AFhK9saUTXHbdd9ue26L7dd9+W2677cdt1Xt9t2kVKqdg2SJEmS1oGafsiXJEmSpM5j+JckSZJywvAvSZIk5YThX5IkScqJ9Tr8R8RpEbEwIlZExKyI2KPaNaltEfGJiPhNRLwcESkiDq12TapMRJwTEY9HxLKIWBwRt0fEttWuS22LiFMi4s8RsbTwejQiDqh2XWqfiDi78Hvz0mrXorZFxNTC9ip9PVvtulSZiNg8In4eEW9ExLsR8ZeIGFftutqy3ob/iDiK7AFh5wG7Ak8B90XE0KoWpkpsTLa9Tqt2IWq3vYArgY8C+wEbAL+LiI2rWpUq8RJwNrAbMA54ALgjInaoalWqWETsDnwR+HO1a1G7/BUYXvLas7rlqBIRMRh4BHgfOAAYC3wFWFLNuiqx3t7qMyJmAY+nlL5UeN8DeBH4cUrpgqoWp4pFRAIOSyndXu1a1H6FB/QtBvZKKf2h2vWofSLiTeD/pZSmVbsWtS4i+gFPAqcC5wJzU0pnVLcqtSUipgKHppR2qXYtap+IuAD4WErp49Wupb3WyyP/EdGb7OjV74ttKaXVhffjq1WXlEMDC/++WdUq1C4R0TMiPkt2Fu7RatejilwJ/Dal9Ps2x1StGVPo5vp8RNwYEVtWuyBV5NPAnIi4pdDN9U8RcVK1i6rEehn+gQ8APYFXy9pfBYat+3Kk/CmcbbsUeCSl9HS161HbImKniFhO9qTKn5CddXumymWpDYUvarsC51S7FrXbLGAKsD9wCjAK+N+I6F/NolSR0WTbbAEwCbgauDwijq9qVRXoVe0CJK23rgR2xP6r3cl8YBeyMzZHAD+LiL38AlC7ImIEcBmwX0ppRbXrUfuklO4pefvnQpflfwCfAexuV9t6AHNSSl8vvP9TROwInAz8rHpltW19PfL/OtAAbFbWvhnwyrovR8qXiLgCOBjYO6X0UrXrUWVSSu+llP6eUnoipXQO2YX3/1HtutSq3YChwJMRsSoiVpFdeP/lwvue1S1P7ZFSegv4G7B1tWtRm+qA8gMj84Ca77a1Xob/lNJ7wBPAPsW2QheEfbD/qtRlInMFcBjwyZTSC9WuSWulB9Cn2kWoVTOAncjO2BRfc4AbgV1SSg1VrE3tVLhw+0NkwVK17RGg/FbW25Cdualp63O3nx+SnbKeA8wGziC7eO26qlalNhV++ZUe9RgVEbsAb6aUFlWpLFXmSuBzwCHAsogoXmNTn1J6t3plqS0R8X3gHmAR0J9sO04k68uqGpVSWgY0uaYmIt4G3vBam9oXERcDvyELjB8kuz15A3BTNetSRX4E/DEivg78EtgD+ELhVdPW2/CfUrq5cJvB75Bd5DsX2D+lVH4RsGrPOODBkvc/LPz7M7ILo1S7Tin8O7Os/QRg+jqtRO01FLie7D7j9WT3ip+UUrq/qlVJ67ctyIL+psBrwMPAR1NKr1W1KrUppfR4RBwGfB/4FvACcEZK6cbqVta29fY+/5IkSZKaWi/7/EuSJElak+FfkiRJygnDvyRJkpQThn9JkiQpJwz/kiRJUk4Y/iVJkqScMPxLkiRJOWH4lyRJknLC8C9JkiTlhOFfkiRJygnDvyRprUTEkIh4JSK+XtI2ISLei4h9qlmbJKmpSClVuwZJUjcXEQcCtwMTgPnAXOCOlNKZVS1MktSE4V+S1Cki4kpgX2AOsBOwe0ppZXWrkiSVMvxLkjpFRGwEPA2MAHZLKf2lyiVJksrY51+S1Fk+BHyQ7G/LyOqWIklqjkf+JUlrLSJ6A7PJ+vrPB84AdkopLa5qYZKkJgz/kqS1FhE/AI4APgwsBx4C6lNKB1e1MElSE3b7kSStlYiYSHak/9iU0tKU0mrgWODjEXFKVYuTJDXhkX9JkiQpJzzyL0mSJOWE4V+SJEnKCcO/JEmSlBOGf0mSJCknDP+SJElSThj+JUmSpJww/EuSJEk5YfiXJEmScsLwL0mSJOWE4V+SJEnKCcO/JEmSlBOGf0mSJCkn/n+zGkaslC5IlgAAAABJRU5ErkJggg==\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 }