{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Step 3: Diffusion Equation in 1-d\n", "\n", "Next, we tackle the one-dimensional diffusion equation:\n", "\n", "$$\\frac{\\partial u}{\\partial t} = \\nu \\frac{\\partial^2 u}{\\partial x^2}$$\n", "\n", "First obvious difference is that — unlike our previous two simple equations — this equation has a second-order derivative. Before continuing we must learn how to discretize it.\n", "\n", "## Discretizing $\\frac{\\partial^2 u}{\\partial x^2}$\n", "\n", "Any second order derivative can be understood geometrically as the line tangent to the curve given by a first derivative. To discretize this second order derivative we can use a Central difference scheme, a special combination of the two methods presented earlier Forward Difference and Backward Difference. \n", "\n", "Consider the Taylor series expansions of $u_{i+1}$ and $u_{i-1}$ around $u_i$:\n", "\n", "$$u_{i+1} = u_i + \\Delta x \\frac{\\partial u}{\\partial x}\\biggr\\rvert_i + \n", " \\frac{\\Delta x^2}{2} \\frac{\\partial^2 u}{\\partial x^2}\\biggr\\rvert_i + \n", " \\frac{\\Delta x^3}{3!} \\frac{\\partial^3 u}{\\partial x^3}\\biggr\\rvert_i + O(\\Delta x^4)\n", "$$\n", "\n", "$$u_{i-1} = u_i - \\Delta x \\frac{\\partial u}{\\partial x}\\biggr\\rvert_i + \n", " \\frac{\\Delta x^2}{2} \\frac{\\partial^2 u}{\\partial x^2}\\biggr\\rvert_i - \n", " \\frac{\\Delta x^3}{3!} \\frac{\\partial^3 u}{\\partial x^3}\\biggr\\rvert_i + O(\\Delta x^4)\n", "$$\n", "\n", "If we add together both these expansions, the odd ordered derivative terms will cancel out. If we neglect any terms of $O(\\Delta x^4)$ or higher because they are so small we can rearrange the sum of these two expansions to solve for our second-derivative.\n", "\n", "$$u_{i+1} + u_{i-1} = 2u_i + \\Delta x^2 \\frac{\\partial^2 u}{\\partial x^2}\\biggr\\rvert_i \n", "$$\n", "\n", "$$\\frac{\\partial^2 u}{\\partial x^2} = \\frac{u_{i+1} - 2u_i + u_{i-1}}{\\Delta x^2}\n", "$$\n", "\n", "## Discretizing the full equation\n", "\n", "We now have all the tools necessary to write the discretized version of the diffusion equation in 1D:\n", "\n", "$$\\frac{u^{n+1}_i - u^n_i}{\\Delta t} =\n", "\\nu \\frac{u^{n}_{i+1} -2u^n_i + u^n_{i-1}}{\\Delta x^2}\n", "$$\n", "\n", "With only one unknown $u^{n+1}_i$ we can now rearrange this equation and obtain our final result.\n", "\n", "$$u^{n+1}_i = u^n_i + \\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", "## Libraries & Initial Conditions" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEWCAYAAACJ0YulAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvFvnyVgAAHZ1JREFUeJzt3Xm0bGV55/Hvr4ZbJYMM3hvDICodIx2iOFzU2K72qokCiVFX2tXSKgExxLTaZmknGrUlK9ptm8S0bRxYBAmaAbWVOC21pRsIbWjQiwsZNBAEMRcVLiDI0HWo4ek/9q596sxVdWpXnffU77PWWbdO7V1Vz91n7/3U+7z73a8iAjMzM4DKrAMwM7Otw0nBzMwKTgpmZlZwUjAzs4KTgpmZFZwUzMys4KRgc03SDZL25I8l6S8l/UTSN/LnfkfSHZIekPSoTXzOA5KOnVDYZqVxUrAtQdIbJO2VtCDpgmXL9kjq5SfWByTtk/RpSSeu836PkxQDr7lD0pck/crgehFxfERclv/6HOBXgKMj4hmS6sCfAS+MiIMi4u5x/3/5628Z9/VrkfR9Sb885mufIulqSQ/l/z5l0vFZepwUbKv4IfAe4Py1lkfEQcDBwLOAfwT+j6QXbPC+h+avOwG4GPg7Saevse5jge9HxIP5748GmsANQ/8vEiFpB/B54K+Bw4CPA5/Pn7c55qRgW0JEXBQRnwPW/TYemX0R8S7gPOB9Q77/jyPivwN/CLxPUgUWv2lLOjN/v1/KWxYXAjfmL79X0iUDrY9a/30lXSbptfnjn5P095Luk3SXpE8NrBeSfi5/fIikT0jaL+k2Se8ciOd0SV+X9Kd5GetWSSev9n+S9FfAMcAX85h/f5htkdsD1IAPRMRCRHwQEPD8Ed7DtiEnBUvZRcDTJB044mt+Bnji4JMR8THgdcD/zUs9pwLH54sPjYhhTpbvBr5G9s37aODP11jvz4FDgGOB5wKnAWcMLH8mWULaCfwx8DFJWv4mEfFq4AfAi/OY/xhA0r3r/Lwtf/nxwLWx9D431w78n21O1TZexWzL+iHZt9tDgQc3WHfwNQCHlxBPm6wEdWRE7AO+vnwFSVXgFcBTIuJ+4H5J7wdeDXwsX+22iPiLfP2PAx8hK2X9eJggIuLQIVY7CLhv2XP3kZXnbI65pWApOwoI4N4RXwNwz+TD4ffJktQ38quaXrPKOjuBOnDbwHO3DcQFAyf/iHgof3jQhGN9AHjksuceCdw/4c+xxDgpWMpeBnxroGN42NfcyWJ/wSj6n3PAwHM/23+Q91v8VkQcCfw28JF+P8KAu1hsUfQdA9w+RjyQJcUlBq64Wu3n7flqNwBPXlaWejLbsFPdRuOkYFuCpJqkJlAFqpKagx26A+tJ0lGSzgZeC7x9+TprvP+jJb0BOBv4g4jojRpjROwnO3m/SlI1bwn8i4HPeLmko/Nff0J2wu4te48u8GngP0s6WNJjgTeTXQU0jjvI+iYGP+OgdX7+S77aZUAX+A+SGvm2AbhkzDhsm3BSsK3incD/A94GvCp//M6B5UdKeoCs7PFN4EnAnoj42gbve6+kB4HrgFOAl0fEWpe9DuO3gN8ju0rqeOCKgWUnAlflcX4BeNMaYxPeSNbquIWs3+FvWftS3I28F3hn3on8H4d9UUQ8DLyUrJP7XuA1wEvz522OyZPsmJlZn1sKZmZWKC0pSHqMpEslfSe/EuNNq6wjSR+UdLOkayU9rax4zMxsY2WOU+gAb4mIb0k6GLha0sUR8Z2BdU4GnpD/PBP4aP6vmZnNQGkthYj4UUR8K398P/Bdll6LDfAS4BP5rQuuBA6VdERZMZmZ2fqmMqJZ0uOApwJXLVt0FPDPA7/vy5/70bLXnwWcBXDggQc+/bjjjisrVMs99HCH7+0f5fJ/Azjm8AM45BH1WYdhtsLVV199V0Ts2mi90pOCpIOAzwK/GxE/Hec9IuJc4FyA3bt3x969eycYoa3m0hvv5Iy//CZ/cdpufvGo5QNfbbk7frrASz/8D7znN57Evz3xmFmHY7aCpNs2XqvkpJDfj/6zwN9ExEWrrHI78JiB349m/JGdNkEL7S4ARx7a5IhDHjHjaLa+HdWsEttqjzwmzmxLKfPqI5Hd4Ou7EfFna6z2BeC0/CqkZwH3RcSP1ljXpqh/cmvWqzOOJA397dTKk6lZqspsKfwrsjs/Xifpmvy5t5Pd54WIOAf4Mtko05uBh1h6+2Cbof7JzUlhOItJwS0FS1tpSSEivk52x8j11gng9WXFYOMrkkLN4xuHUa2IelW0Om4pWNp8xNuqWh2Xj0bVrFVdPrLkOSnYqlw+Gl2jXnX5yJLnpGCrarV71KuiWlm3AmgDmvVKcdWWWaqcFGxVrXaXZs2thFE061X3KVjynBRsVQudLg2XjkbSrFdcPrLkOSnYqlrtHs26d49RuKPZtgMf9baqVrvrTuYRNetOCpY+JwVbVZYUvHuMwuUj2w581NuqWu2eO5pH1HBHs20DTgq2qlbH5aNRNWtVFtxSsMQ5Kdiq3NE8uqx85JaCpc1Hva1qoe1LUkfljmbbDpwUbFUevDa6Zr1S3DPKLFVOCraqVsflo1E1a1W6vaDddWKwdPmot1V5nMLoPNGObQdOCrZCRHicwhj628tjFSxlPupthXY36AXuUxhRwy0F2wacFGyFhY7nUhhHf3steACbJcxJwVbolz9cPhpNf+pSl48sZT7qbYV++cPjFEbjjmbbDpwUbAWXj8azWD5yS8HS5aRgKxTlo5p3j1EsXn3kloKly0e9reDy0XgatX75yC0FS5eTgq3glsJ43FKw7cBHva3QP6m5T2E0RUezL0m1hDkp2AotdzSPpenykW0DTgq2gscpjKfh8pFtAz7qbQWXj8bTqFWQsrkozFLlpGArFEnB9z4aiSQaNc+pYGlzUrAV+oOvGi4fjcyzr1nqfNTbCq12Fykrh9homjUnBUubj3pbodXu5vVxzTqU5DTrFV99ZElzUrAVWu2eO5nH5PKRpc5JwVZotbvuZB5To151R7MlzUnBVmh1eh6jMKZmreKWgiWttCNf0vmS7pR0/RrLD5H0RUnflnSDpDPKisVGk83P7JbCOJr1qscpWNLK/Dp4AXDSOstfD3wnIk4A9gDvl7SjxHhsSK1213dIHZM7mi11pSWFiLgcuGe9VYCDlV3iclC+bqeseGx4C+2e75A6pma96hviWdJmeeR/CPiXwA+B64A3RcSqX7EknSVpr6S9+/fvn2aMc6nVcfloXB6nYKmbZVJ4EXANcCTwFOBDkh652ooRcW5E7I6I3bt27ZpmjHMp61NwS2EcLh9Z6mZ55J8BXBSZm4FbgeNmGI/lFjoepzAuj1Ow1M0yKfwAeAGApEcDTwRumWE8lvM4hfE16lUWOj0iYtahmI2lVtYbS7qQ7KqinZL2AWcDdYCIOAd4N3CBpOsAAW+NiLvKiseGl41odvloHP3t5taWpaq0pBARp26w/IfAC8v6fBufxymMb3H2NW9DS5O/DtoSEcFCp+dxCmMq5ml2Z7MlyknBlujPpeDy0XgWy0fubLY0+ci3JTzr2ua4pWCpc1KwJfonM8+6Np7+xES+LNVS5SPflnBLYXMWWwpOCpYmJwVbon/fHl85M55+n4LnVLBUOSnYEv3ykTuax9OouaVgafORb0sU5SO3FMbi8pGlzknBllhMCt41xlFckuqrjyxRPvJtieLqI3c0j6VoKXicgiXKScGWWHBH86a4fGSpc1KwJVw+2pxmMU7B5SNLk498W2Lx6iO3FMZRq1aoVeSWgiXLScGW8NVHm5dNtOOWgqXJScGWKFoKNe8a42rWK+5otmT5yLclWp0utYqoVb1rjKtR85Scli4f+baEJ4fZvGa94nEKliwnBVvCU3FuXtan4JaCpclHvy2x0O564NomNetV9ylYspwUbIlswnnvFpvRrFd89ZEly0e/LeE+hc1ruqPZEuakYEu0Ok4Km+U+BUuZk4It4Y7mzWu4fGQJ89FvS7TaXU/FuUnNerW4saBZapwUbAn3KWxe1qfgloKlyUnBlmi1ezRcPtqU7OojtxQsTT76bYkFdzRvWrNepdMLOl23Fiw9Tgq2RKvdc5/CJhVTcnacFCw9Tgq2RNan4N1iMzz7mqXMR78VOt0enV74Nheb1OjPvuaWgiXIScEK/ZOYWwqb45aCpcxHvxU869pk9FtaTgqWIicFKywmBe8Wm9Hffh6rYCny0W+FYipOtxQ2pb/9FtxSsASVlhQknS/pTknXr7POHknXSLpB0t+XFYsNp99ScEfz5hR9Cr7VhSWozJbCBcBJay2UdCjwEeDXI+J44OUlxmJD6N+vx+WjzXH5yFJW2tEfEZcD96yzyr8DLoqIH+Tr31lWLDYcl48mo+mOZkvYLL8S/jxwmKTLJF0t6bS1VpR0lqS9kvbu379/iiHOF199NBmLl6S6pWDpmWVSqAFPB34VeBHwnyT9/GorRsS5EbE7Inbv2rVrmjHOlcWWgstHm7FYPnJLwdJTm+Fn7wPujogHgQclXQ6cANw0w5jmWtFScEfzprij2VI2y6+EnweeI6km6QDgmcB3ZxjP3Gt1XD6ahOI2Fy4fWYJKaylIuhDYA+yUtA84G6gDRMQ5EfFdSV8FrgV6wHkRseblq1Y+l48mQxKNWsXjFCxJpSWFiDh1iHX+BPiTsmKw0Sy4pTAxzXrVfQqWJH8ltEK/pdAvf9j4stnXXD6y9Pjot8JCu0ujVkHSrENJXrNedUezJclJwQrZBDsuHU1Cs+bykaXJScEKrXbPncwT4vKRpcpnACu0Om4pTErDHc2WKCcFK7TaXQ9cm5CsT8EtBUuPk4IVXD6anKbHKViifAawQqvdpeHy0UR4nIKlyknBCq1Oz30KE+KOZkuVk4IVFtpdmh64NhHNerUYIW6WEp8BrOBxCpOTlY/cUrD0OClYodXu+RYXE9KoVWh1ukTErEMxG4nPAFbwOIXJadarRMDDXbcWLC1OClbIykfeJSbBcypYqnwGMAAiIh+n4JbCJPS3o8cqWGqGmk9B0rtWez4i/miy4disLHT6E+w4KUxCMSWnWwqWmGEn2Xlw4HET+DU8dea2suC5FCaqX4bz7bMtNUMlhYh4/+Dvkv4U+J+lRGQz4fmZJ6t/DymParbUjPu18ADg6EkGYrPVP3k5KUyGy0eWqmH7FK4D+hdcV4FdgPsTtpH+yctXH01GUT5yS8ESM2yfwq8NPO4Ad0REp4R4bEaKloJvnT0Riy0FJwVLy7B9CreVHYjNlstHk7XY0ezykaXFtQIDBi9J9S4xCQ13NFuifAYwwC2FSfPgNUuVk4IBi2UOtxQmY7Gj2eUjS4vPAAYsthQa7mieCHc0W6qcFAxYLHO4fDQZ9WqFakUe0WzJcVIwwOMUytCseUpOS4/PAAa4o7kM2exrbilYWpwUDMjufVStiHrVu8SkeEpOS5HPAAZk5aOm75A6UY16xX0KlhyfBQzoz7rm0tEkNWtVj1Ow5DgpGIBnXStBs+6OZkuPk4IBWZ9Cw1ceTZQ7mi1FPgsYkI1T8B1SJ6tZr7pPwZJTWlKQdL6kOyVdv8F6J0rqSPo3ZcViG8vKR/6OMEnNeqWY5tQsFWWeBS4ATlpvBUlV4H3A10qMw4bgjubJa9bcUrD0lJYUIuJy4J4NVnsj8FngzrLisOG0Ol0aviR1ohruaLYEzewsIOko4GXAR4dY9yxJeyXt3b9/f/nBzSFffTR5jZo7mi09s/xq+AHgrRGx4VepiDg3InZHxO5du3ZNIbT54/LR5DXrVfcpWHKGnaO5DLuBT0oC2AmcIqkTEZ+bYUxzyx3Nk9esV3i426PbC6oVzTocs6HMLClExOP7jyVdAHzJCWF2Ftpdz6UwYcXsa50uB+yY5fcvs+GVtqdKuhDYA+yUtA84G6gDRMQ5ZX2ujafVcflo0vr3kmq1exywY8bBmA2ptKQQEaeOsO7pZcVhG+v2gnY3XD6aMM++ZinyWcBY6HguhTI4KViKnBRscdY1j1OYqH7Ly2MVLCU+C5hnXStJo99S8KhmS4iTgjkplKR/g0GXjywlTgq2WD5yR/NE9benB7BZSnwWsKK80XBLYaLc0WwpclKwxfKRB69NVNN9CpYgJwUryhsuH02Wrz6yFPksYO5oLok7mi1FTgpWlDecFCZrsU/BLQVLh5OC+eqjkjSKex+5pWDp8FnA3NFckkpF7KhV3NFsSXFSsIGWgpPCpDVrFY9TsKQ4KVjRUvAczZPXrHtKTkuLzwJGq9NlR61CxbODTZyTgqXGScFYaPd8h9SSNOsVX31kSfGZwGi1PetaWZr1qjuaLSlOCuakUKJmreqOZkuKk4LRavc8RqEkjbovSbW0+ExgtDpdGh6jUIpGreo+BUuKk4Ll5SPvCmVo1iss+OojS4jPBJaXj9xSKIMvSbXUOCkYrbbLR2Vp1iu0Oi4fWTqcFIyHO+5oLkuz5paCpcVnAvMlqSXql48iYtahmA3FScFouaVQmma9Qi+g3XVSsDT4TGBZS8F9CqXwPM2WGieFORcRLh+VqFH3lJyWFieFOdfuBr3wrGtl6d9o0Le6sFT4TDDnPD9zuZpuKVhinBTmXDHBjpNCKRaTglsKlgYnhTnXL2t4PoVy9Mty7mi2VPhMMOf6LQWXj8rh8pGlxklhzvXLGk4K5ehf6uvykaWitKQg6XxJd0q6fo3lr5R0raTrJF0h6YSyYrG1LXY0+/tBGYrykVsKlogyzwQXACets/xW4LkR8STg3cC5JcZia3D5qFwuH1lqamW9cURcLulx6yy/YuDXK4Gjy4rF1laUjzyiuRSNoqPZ5SNLw1apGZwJfGWthZLOkrRX0t79+/dPMaztb7GlsFV2he2l31LwRDuWipmfCSQ9jywpvHWtdSLi3IjYHRG7d+3aNb3g5oDLR+Va7Gh2UrA0lFY+GoakJwPnASdHxN2zjGVe9csaDbcUSlGviop89ZGlY2ZnAknHABcBr46Im2YVx7xbcEuhVJI8JaclpbSWgqQLgT3ATkn7gLOBOkBEnAO8C3gU8BFJAJ2I2F1WPLa6onzkjubSNOtVj2i2ZJR59dGpGyx/LfDasj7fhtNq96goK3NYOZq1iu+SaslwIXnO9edSyFtrVoKspeCkYGlwUphzrU6Xhm+GV6odtYr7FCwZPhvMuVa7507mkrmj2VLipDDnFjpOCmVr1t2nYOlwUphzrbbLR2Xz1UeWEp8N5ly/o9nK06y5fGTpcFKYcwvtnu97VLJmveIRzZYMnw3mXKvjlkLZ3NFsKXFSmHOtdtejmUvmpGApcVKYcy2Xj0rXqFc8eM2S4bPBnHNHc/matSoPd3r0ejHrUMw25KQw55wUyldMtOPWgiXASWHOtTo9z6VQsn55zv0KlgKfDeZYrxc83Om5o7lk/ZaCB7BZCpwU5li/nOHyUbkWWwouH9nW56QwxxbnZ/ZuUCbP02wp8dlgjvXLGW4plKsoHzkpWAKcFOZYv5zhlkK5Gi4fWUJ8Nphjnp95OtzRbClxUphji30KTgpl6ifdBZePLAFOCnOsX87wOIVy+eojS4nPBnPMHc3T4Y5mS4mTwhxbcJ/CVDgpWEqcFOaYrz6ajqJ85HsfWQJ8Nphj7mieDg9es5Q4KcwxJ4XpqFTEjmrFd0m1JDgpzLHFex95Nyhbo15xS8GS4LPBHCsuSXVHc+kataovSbUkOCnMsVanS70qqhXNOpRtr1mvePCaJcFJYY612l1fjjolzXrVt7mwJDgpzLFWu0fDncxT0axXXD6yJDgpzLGFdtedzFPSrFXd0WxJ8BlhjrU6XV+OOiXNupOCpcFJYY612j23FKbE5SNLhc8Ic8wdzdPTcEezJaK0pCDpfEl3Srp+jeWS9EFJN0u6VtLTyorFVtdqu3w0Lc1alQW3FCwBZbYULgBOWmf5ycAT8p+zgI+WGIutwuWj6Wl6RLMlolbWG0fE5ZIet84qLwE+EREBXCnpUElHRMSPyojnq9f/mLd8+poy3jpZD7W7HLvrwFmHMRcO2FHl7gcf5vh3fXXWoVjCznzO43nzC59Y6meUlhSGcBTwzwO/78ufW5EUJJ1F1poAeEDSjWN+5k7grjFfW6aZxfVh4MOvXHcVb7PROK7ROK4RvAV2vmX8uB47zEqzTApDi4hzgXM3+z6S9kbE7gmENFFbNS7YurE5rtE4rtHMc1yzLCjfDjxm4Pej8+fMzGxGZpkUvgCcll+F9CzgvrL6E8zMbDillY8kXQjsAXZK2gecDdQBIuIc4MvAKcDNwEPAGWXFMmDTJaiSbNW4YOvG5rhG47hGM7dxKbv4x8zMzCOazcxsgJOCmZkVtk1SkHSSpBvz22a8bZXlDUmfypdfNTiwTtIf5M/fKOlFU47rzZK+k9/q439LeuzAsq6ka/KfL0w5rtMl7R/4/NcOLPtNSf+U//zmlOP6bwMx3STp3oFlZW6vsW/bUvL22iiuV+bxXCfpCkknDCz7fv78NZL2TjmuPZLuG/h7vWtg2br7QMlx/d5ATNfn+9Th+bJStpekx0i6ND8P3CDpTausM739KyKS/wGqwPeAY4EdwLeBX1i2zr8HzskfvwL4VP74F/L1G8Dj8/epTjGu5wEH5I9/px9X/vsDM9xepwMfWuW1hwO35P8elj8+bFpxLVv/jcD5ZW+v/L3/NfA04Po1lp8CfAUQ8CzgqrK315BxPbv/eWS3lrlqYNn3gZ0z2l57gC9tdh+YdFzL1n0xcEnZ2ws4Anha/vhg4KZVjsep7V/bpaXwDODmiLglIh4GPkl2G41BLwE+nj/+DPACScqf/2RELETErWRXQz1jWnFFxKUR8VD+65Vk4zXKNsz2WsuLgIsj4p6I+AlwMevf46rMuE4FLpzQZ68rIi4H7llnleK2LRFxJXCopCMod3ttGFdEXJF/Lkxv/xpme61lM/vmpOOayv4VET+KiG/lj+8Hvkt2d4dBU9u/tktSWOuWGauuExEd4D7gUUO+tsy4Bp1J9m2grylpr6QrJb10QjGNEtdv5E3Vz0jqDzTcEtsrL7M9Hrhk4Omyttcw1oq9zO01quX7VwBfk3S1slvJTNsvSfq2pK9IOj5/bktsL0kHkJ1cPzvwdOnbS1lZ+6nAVcsWTW3/SuI2F/NA0quA3cBzB55+bETcLulY4BJJ10XE96YU0heBCyNiQdJvk7Wynj+lzx7GK4DPRMTgrUdnub22NEnPI0sKzxl4+jn59voZ4GJJ/5h/k56Gb5H9vR6QdArwObI7Jm8VLwb+ISIGWxWlbi9JB5Elod+NiJ9O6n1HtV1aCsPcMqNYR1INOAS4e8jXlhkXkn4ZeAfw6xGx0H8+Im7P/70FuIzsG8RU4oqIuwdiOQ94+rCvLTOuAa9gWdO+xO01jLVin/ntXCQ9mexv+JKIuLv//MD2uhP4OyZXNt1QRPw0Ih7IH38ZqEvayRbYXrn19q+Jby9JdbKE8DcRcdEqq0xv/5p0p8ksfshaPLeQlRP6nVPHL1vn9SztaP50/vh4lnY038LkOpqHieupZB1rT1j2/GFAI3+8E/gnJtThNmRcRww8fhlwZSx2bN2ax3dY/vjwacWVr3ccWaefprG9Bj7jcazdcfqrLO0I/EbZ22vIuI4h6yd79rLnDwQOHnh8BXDSFOP62f7fj+zk+oN82w21D5QVV778ELJ+hwOnsb3y//cngA+ss87U9q+JbehZ/5D1zt9EdoJ9R/7cH5F9+wZoAv8jP0C+ARw78Np35K+7ETh5ynH9L+AO4Jr85wv5888GrssPiuuAM6cc13uBG/LPvxQ4buC1r8m3483AGdOMK//9D4H/uux1ZW+vC8lu694mq9ueCbwOeF2+XGR3Iv9e/vm7p7S9NorrPOAnA/vX3vz5Y/Nt9e387/yOKcf1hoH960oGktZq+8C04srXOZ3s4pPB15W2vchKegFcO/B3OmVW+5dvc2FmZoXt0qdgZmYT4KRgZmYFJwUzMys4KZiZWcFJwczMCk4KZmZWcFIwM7OCk4LZJkk6Mb9xYFPSgfk98X9x1nGZjcOD18wmQNJ7yEbNPwLYFxHvnXFIZmNxUjCbAEk7gG8CLbJbNnQ3eInZluTykdlkPAo4iGzmrOaMYzEbm1sKZhOQzwn9SbK7ex4REW+YcUhmY/EkO2abJOk0oB0RfyupClwh6fkRcclGrzXbatxSMDOzgvsUzMys4KRgZmYFJwUzMys4KZiZWcFJwczMCk4KZmZWcFIwM7PC/wd4rqyi1lkoNQAAAABJRU5ErkJggg==\n", "text/plain": [ "