{ "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": [ "