{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example of stiff ODE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider the following ODE:\n", "\n", "$$\\frac{df}{dt} = f^2 - f^3$$\n", "\n", "We solve this ODE with the initial condition $f(0) = \\delta$ over the time interval from $t = 0$ to\n", "$t = 2/\\delta$. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true, "execution": { "iopub.execute_input": "2026-06-11T00:11:07.728589Z", "iopub.status.busy": "2026-06-11T00:11:07.728426Z", "iopub.status.idle": "2026-06-11T00:11:09.178694Z", "shell.execute_reply": "2026-06-11T00:11:09.177848Z" } }, "outputs": [], "source": [ "## IMPORTS ##\n", "from __future__ import division, print_function\n", "import matplotlib.pyplot as plt, numpy as np\n", "%matplotlib inline\n", "import PDSim.core.integrators as integrators # The abstract integrators with callback functions\n", "from PDSim.misc.datatypes import arraym # An optimized list-like object with rapid element-wise operators" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true, "execution": { "iopub.execute_input": "2026-06-11T00:11:09.180654Z", "iopub.status.busy": "2026-06-11T00:11:09.180399Z", "iopub.status.idle": "2026-06-11T00:11:09.185328Z", "shell.execute_reply": "2026-06-11T00:11:09.184496Z" } }, "outputs": [], "source": [ "class TestIntegrator(object):\n", " \"\"\"\n", " Implements the functions needed to satisfy the abstract base class requirements\n", " \"\"\"\n", "\n", " def __init__(self):\n", " self.x, self.y = [], []\n", " \n", "\n", " def post_deriv_callback(self): \n", " \"\"\" Don't do anything after the first call is made to deriv function \"\"\"\n", " pass\n", "\n", " def premature_termination(self): \n", " \"\"\" Don't ever stop prematurely \"\"\"\n", " return False\n", "\n", " def get_initial_array(self):\n", " \"\"\" The array of initial values\"\"\"\n", " return arraym([self.delta])\n", "\n", " def pre_step_callback(self): \n", " if self.Itheta == 0:\n", " self.x.append(self.t0)\n", " self.y.append(self.xold[0])\n", "\n", " def post_step_callback(self): \n", " self.x.append(self.t0)\n", " self.y.append(self.xold[0])\n", "\n", " def derivs(self, t0, xold):\n", " \n", " dfdt = xold[0]**2 - xold[0]**3\n", " \n", " return arraym([dfdt])\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we define the actual concrete implementations of the integrators, which are formed of the common functions and the abstract integrator" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2026-06-11T00:11:09.186910Z", "iopub.status.busy": "2026-06-11T00:11:09.186751Z", "iopub.status.idle": "2026-06-11T00:11:09.190254Z", "shell.execute_reply": "2026-06-11T00:11:09.189569Z" } }, "outputs": [], "source": [ "class TestEulerIntegrator(TestIntegrator, integrators.AbstractSimpleEulerODEIntegrator):\n", " \"\"\" Mixin class using the functions defined in TestIntegrator \"\"\"\n", " pass\n", "\n", "class TestHeunIntegrator(TestIntegrator, integrators.AbstractHeunODEIntegrator):\n", " \"\"\" Mixin class using the functions defined in TestIntegrator \"\"\"\n", " pass\n", "\n", "class TestRK45Integrator(TestIntegrator, integrators.AbstractRK45ODEIntegrator):\n", " \"\"\" Mixin class using the functions defined in TestIntegrator \"\"\"\n", " pass\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As an example, we assume $\\delta = 0.01$. After the transition from 0 to 1, the time step is very small and the computation goes slowly. For smaller values of $\\delta$, the situation is even worse resulting in overflow. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2026-06-11T00:11:09.192000Z", "iopub.status.busy": "2026-06-11T00:11:09.191841Z", "iopub.status.idle": "2026-06-11T00:11:09.310662Z", "shell.execute_reply": "2026-06-11T00:11:09.309874Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcmJJREFUeJzt3Ql4U1X6P/Bvkjbpvre0Zd9EdgRFcV8RVFxncHTcxZ+7ouOMor8ZRf2POs4orrjruOJvZtxlUFxQtC7IolCQfSm0pXRf02z3/7wnuWmSJm1aSpu038/zhCQ3J7dJb8l9855z3mPQNE0DERERUQ8x9tQPJiIiIhIMRoiIiKhHMRghIiKiHsVghIiIiHoUgxEiIiLqUQxGiIiIqEcxGCEiIqIexWCEiIiIelQMooDL5UJxcTGSk5NhMBh6+uUQERFRGKSual1dHfLz82E0GqM7GJFAZODAgT39MoiIiKgTioqKMGDAgOgORiQjor+ZlJSUnn45REREFIba2lqVTNDP41EdjOhdMxKIMBghIiKKLu0NseAAViIiIupRDEaIiIioRzEYISIioh4VFWNGwp0+5HA44HQ6e/qlUBczmUyIiYnhtG4iol6qVwQjNpsNJSUlaGxs7OmXQgdIQkIC8vLyYDab+TsmIuploj4YkYJo27dvV9+epaiKnKxYGK33kIyXBJv79u1Tx3nkyJFtFs4hIqLoE/XBiJyoJCCReczy7Zl6n/j4eMTGxmLnzp3qeMfFxfX0SyIioi7Ua75i8tty78bjS0TUe0V9ZoSIqDOcLg0/bq9EWZ0VOclxmDo0AyYj174i6gkMRoioz1myrgTzP1yPkhqrd1teahzunjUGM8bl9ehrI+qLek03TVd8S/puawXeX7NHXcv9SLds2TI1WLe6urqnXwpRVAUi176+yi8QEaU1VrVdHiei7sVgxPPhdPRDX+CC57/HzYvWqGu5fyA/lC677DIVSAReZsyYgUhktVrVax4/fryq+XH22WcHbffVV19hypQpapDpsGHD8Mwzz7Rq85///AdjxoyBxWJR1++++243vAMi95cOyYgE+6qhb5PH9S8j0fglhSga9fluGv1bUuBHjP4taeFFkw9Y2lYCj5dfftlvm5ygu5PMTgmndocUk5NZLTfddJMKJoKRqbennXYarrrqKrz++uv49ttvcd111yE7OxvnnXeeavPdd9/h/PPPx3333YdzzjlHBSKzZ8/GN998g8MPP7zL3x+RLxkjEpgR8SWfA/K4tKtpsrErh6ibGHtjXYpGmyOsS53Vjrs/KGzzW9I9H6xX7cLZn/zsjpDAIzc31++Snp6uHtuxY4fKlKxZs8bbXrpjZJt0z4RSUFCAY489VgUOMt1ZgoeGhgbv40OGDMH999+vshypqakqcAhHYmIiFi5cqNrL6wxGsiCDBg3CggULMHr0aMyZMwdXXHEF/v73v3vbyGOnnHIK5s2bh4MPPlhdn3TSSWo70YEmg1XDsXR9adhdOcyeEO2/XpcZabI7MeYvn3TJviS0KK21Yvw9n4bVfv29pyLB3HO/0rVr1+LUU09VWYcXX3xRFQq74YYb1MU3A/Pwww/jz3/+M/73f//Xu02CHGkjQUpnSdZj+vTpftvk9chrsdvtqlaItLnllltatWEwQt1BZs2E4701xSG/pBg8XTmnjMlVQQsHwhLtv16XGYkmH330EZKSkvwuEkh0lgQZF154IebOnasqlR555JF4/PHH8eqrr6oxH7oTTzwRt912G0aMGKEuYtSoUSpTsj9KS0vRr18/v21yX9YMKi8vb7ONbCc60GT6rsyaCTWBV7ZnJMaissHWblfOk19sYfaEqIv0usxIfKxJZSjCIf3Cl728ot12r1x+mPoQC+dnd8QJJ5yguj58ZWS0/3NCWblyJbZs2YI33njDu026jvSS+dJ1Ig499NBWz/3111/RFQJL8etdV77bg7VhCX/qDlJHRKbvSleL/BX6Zj/0v8pzJvXHi9/uaHdfL3+7vUuzJ+HWPWF9FOqNel0wIie1cLtKjhmZrT4UpB842IeKfAzkpsapdgeiGJKMw9AzE6EqjvqOQ5GujrZI0HH11VercSKBZCyH7889EGQsSWCGo6ysTM2+yczMbLNNYLaE6ECRAEAGpgcGCLmeACE13hxWMFLdZA8re7Lgs03tDpAPt+5JuO0YsFC06XXBSFd/S5LHe6Iqo8xAEbIa8SGHHKJu+w5mDWby5MkoLCwMGeAcaNOmTcOHH37ot+3TTz9VmRgZL6K3Wbp0qd+4EWkjXUpE3UVO3JK1GPOXJWh2uJBoNuGb209U/9flRN7el5TU+Ng2g5GOZE9cLuD6N9uf0RfuzL+uDFgY1FB36dPBSDjfkg5kNcbm5uZWWQLJImRlZanZMEcccQQefPBBNQNGxlz4DjgN5vbbb1fPuf7669WsF8mAbNiwQZ38n3jiiTafKzNbHnjgATXdNpT169erqcCVlZWoq6vzBkeTJk1S19dccw2efPJJ3Hrrrerny2BVGbz61ltvefdx8803q9k+Dz30EM466yy8//77+Oyzz9TUXqLu5NI0FYiIBpsTDpcLJqMprC8plx81BI9+trlLsif/+/66dgOWEw/u12Z9lM4ENu0FLN0d1DDw6dv6fDDi+y2pu9epWLJkCfLy/IMdGUiqj9946aWX1NRYySzI9r/97W+tZqv4mjBhgio6dtddd+GYY45RXTzDhw9XdT3as3HjRtTU1LTZRmqIyMq5Oj1jo3clDR06FIsXL1ZZj6eeegr5+flqAK1eY0RIBmTRokUqsJIZPfL63n77bdYYoW5XGxAolNfb0D8t3u9Lyj0fFKK0trnVlxT5vFi0oqhLsifhDJZ97bsdYdVHCSewCSdgEV2VhYnWwIfBUfcyaB0tjtEDamtr1UwPOVmmpKT4PSazRGRwppwIubR878XjTF1te3kDTvh7S82e964/CpMGpvm12bavHif+4yvv/V/uno6UeHeXo95tghDZk7knjwwrexKOS6YNxqvftXwR2B8ZieaQAZC89n4pUnjRoMoahGojQdmfTx8TNKjR339bQU1H23R34BOtwZEzAhd/bOv87YuZESLqk2oCshZlQU6+FQEn7bK6Zm8womdP/vJ+odqu65cah3vCzJ6kq2nE7WdPBmckoKu0l4nxzQTtTxZGskpya3/bdGU2p6vaRGJwtCTKBzczM0JRgZkR6mpfb9qHS1760Xv//50zDr8/fLBfm49+KcYNb6723n/tyqlqdp2vFTsq8dtnvvPe/2TusRiVm6xuhxp0qn/0P3XhIbjv4w3tzuj76o8n4LiHv+ySwCYadUU2p6vadHdWqLuzS12NmREiorY+JK3+J+59PtkN3d6ALEGwcRt7A05cJTVN3mBEz57c9e46vyxLZpIZ9589Tj1uNBq83T3ByInCHGP0DqpFiBPO/WeNazewidaApSuyOV3VpjuzQvd0c3bpQE7YaA8rsBJRn1Tb5Gg3GAlcy0Y+uAMVVzf53d9d5X9fPuDnzTzYb9tfzmj5JqoHLPLt35fc9z1ByPWdp7kLF/qSb+rS7rQJ+SpgCRWI6AFLexVoc1MsyE1pv0ptXxZOcBQqu9ITbUraCaACV6vuCQxGiKhPChwzEjQY8XxTljoketYjUHG1tc1gROwJaFMUJGC5+eSRfttuOmlEq2+qeWn+a+tce9xwVR/FN2CZe5L/fjoasNxz5ljcc2b3BDUMfLpPZRiztmQsSU9hMEJEfbqbRp/Ou68+WDeNO4gYPyA1ZDeNnhkZkukeZLq7qrFVmyLPNn3JiF0VrdsEbttR3rrNtn0tK3ALWVkhcPBhekCG5bdTBrQKWKROSqiARR6Xy/UnDD/gQU2kBT59PTgqC3NV6wOBwQgR9ek6I8OyE9vopnFvm+iZ8lsSkOEQxZ5sib5+VbDMSFGlO7A4fJi7zc5K/6BCbatwbxubn+KdehxIphrDO9AS2FLmvu9L35bgyebISuaBAYseFOlbjxmZ5RewiESL/2TLE0ZltwpqfnfYwHaDmsuPGtxuG8nwRELgE43BUW4XBlDhrmp9IDAYIaI+qdbqHjMyIifJG4wEll3SMyOH6MFIG900U4dmthuMHD0iK2RmZIdn2wmjcjz3gwQjngDllDHutZy2eoKTYMHICQfnhAxYNu2tU9dHeV6PpPADA5ZfS9xtxvVPCdlGX/MyO9mdjRmRk9gqqDEa3KeZlDh3cDOqX1KrNrEm947iYtxtD85NbhX4nHNI/3aDmnCCI6nZ0l6bq44ZGjXB0T0dCKBCkXbyeDgLwkZMMPL1119j1qxZqrqmLEr33nvvtfscqQo6ZcoUVZRs2LBheOaZZzr7eomIunTMyPBsdzAipeHrmlsGtTbZnKjzBCwTBqR5A5iGgDZ6X/zhng/y8vpmWO1Obxubw4UST1CjTwuW+82OljYycFAPUE44ONsbwMhzdRIo6d00UsNE7KxohN3Z0kZs9gQfp3lO5PIcR0CbjZ5gZNbEPO9zAtv8Wlqrrs85ZIDnfl2rNuuL3W0uPsLd7bOnqnXmaN0ed2XnC6YOCjpeRvziaXPelAF+wZsv/XeqB3RxsUY15dk3qNHL+x86OF1dx8ca8XVAG7vTfdoe58lABWvj8AzkHO7JmgW2ketTx/ZrN6g5c0Jeu23Om9w1QdbFbQRZegAVTE+vw9bpYKShoQETJ05Ua5CEQ6qjShlxKU++evVq3HnnnWpV2f/85z/oyy677DKcffbZrbYvW7ZMBXnV1dWIxFof8rrHjx+v1tAJ9vrDDT7l+I8ZMwYWi0Vdv/vuu93wDqgv+a74O5z13lnqOtjtUtsvql21Vojk4Y/AlLDZO2BV2pz34dlqm3RpyLfGZE+3hYwb0ffx361feQe4DkiP97bRsyPS7uz3z4Ix3r2fg/olqbaSgCmqbPLuZ/GWr2FzulSGYOKANNVGzoe7Khu9bZZs+xr1zQ7I+UICH+mGkZOmBCR6m8+2L1fBkDj2oCz1M2W/O3zafLnzG/WzxUmj+6n9SNDz/sZl3t+RBEpbPYGPnHRlPw7zRsx61/24kMBEAhRx2vg8JKRshWHgw3jv1y+9x8Dl0rwBy5ABe5A0/BE0x/yK7eX1fkHW2t3uYOTcyf3VfrQBD+PdDS37ET8XuT8TT5hUrfZjj93YKnu0pqhaHbPqjP+HhNStaLK7sD2gzepdVapNU86DIdus2ulu4+r/N2+brT6vWR3jpp+RMOwRnDG1DhZPZuf5Sw71Bizyeyqw3aH2c8kRg71tnrvY3UY/Hjsb3Wt8ScBh8WSGnrloClIzdniPhz6+SQIX3zZ++2lY4w0w9QzT9TM1LNxyjWojbY8a7s7e6eS1JY94FLNP3uxtFzXByMyZM3H//ffj3HPPDau9nIhk+foFCxZg9OjRmDNnjlpv5e9//zsizp5VwCtnuK+pFafTqRbwk2Dy5JNP7nTwKQvoyXo5F198MX7++Wd1PXv2bPzwww/8rdN+0T+YC/YU4LFVj2FbzTYsWLkg6O19MRIAu/BJ8UuAuQyWnE9UFVY5OUqbovodaltOill9QZBvmaK4utG7j5cKn1ZzEfLT4lWb/unx3kGsgfvpnx6n2gzKdH/b3llR793Ps2tlIUsNAzMSEGMyYqjnG7mMEdHbPL2mpU1crMk71mVLWZ23zeOrH3O/Hgme4mIxsp8767OptNbb5tGV7jZZSRZ1GdlPaqJoeGHdU+pxabd5b53K1sj6OjLA96DcJPUedjfsUI/Le5MxLZKJkGBmaGYCEnI/hclShufXPeXt7pJgSrJN5hgD3tv1Agye3/Mvu1u+bBXXWFUNlhijAWPyUpCY59lP4ZPe/cjASmlnNGj4pORl734ksNBVNdhUkCPb91p3ITH3U/W+JLDQ1Vnt2Li31n2sm3chKXepavPTjiq/bFdhcY1qU95c5G2zwqeN1e7ADu3f6nXuMb6DwzxjgVZ5Xo+87gWrHkOzoUTt58LDB+Hw4e6Mzk87K71/G/L73mj7P7V/aXOkJ1j4cXuF93G5/m5bhTe7pHetycwX3/2saXhT7eeCwwbh6JHSRsPLG1qOqcvlDkrF7TNGYWz/ZPXatNi9+KL0DW+7nlohZr8qsMp/LPlGG+obspAVWmVBtccek/8AbvIcOfk0NjZ6l5aPiLVpFv8J+PFZ4PBrgJkP4UCSDINkPwK7uSQzcsIJJ6Cqqgppae7UcEFBAe644w6sWLFCregrK+vKCruyKm+o4yDPlQBQfs6OHTvU70cCAlm9V076I0eOVIHitGnTuvT1y8rBH3zwgVotWCer+UrQIUGIkEBEjul///tfb5sZM2YgPT3db4VfX6zASu2Rj7ILPr4AhRWFGJw8GDvr2l/LxV4/HLFJW733J6QdrwKK/25v+dtMdkzFBZMn4f01xerkOnmYE4W1LevV2GomYlDyQJw5KR8f/VKiTtLHj8pGWkqt334ytKn4zcRJWPxLCbaWN2DCEAc21n/tt5/BKQMxa2I+lqwrVV0n4wfbsalhuf/PShmIMyfm49PCUmzcW4+xg+zY0ujfZmDyQJw1KR+frd+LDaV1GD3Qhm1NLStj26onYkDyQDUO4/MNZVi/bxvMaT97Hx+XejxWbjGq38V5kwfgw/VrUewo8D4+c+hMNDWk4ZP1e1XW6PCDNL/3Ko/La9i8tx5LCkuRkVKLBvNP3scHxh6FmaPHqttby+qxeF0pspMsOGo0gu5n+74GfLS2BKnJNbBaVnofz4+ZhjPGjFe3JUMkr9Oc+rPf72JM1lCcODrH2/3zfmFAm+qJODhrqHccjmS13l37s9/vQ9oclDkEp451d4/9XLIF3+/73Pv4iPhjULhLgr9kzByXi6K6Iizevtj7+IwhM1FWmaSCiuFZiRg/1OH3PlF3COZMO1Rldr7dWoG8rAZU4ge/9xHjzMJVxw5TWaTlW8rV7K1Jw5x++3HWTsJV0w5TwdTyHRv93uex+dOx9BcHJEFz1THD8N2uTfjV5+9P98zJz+Co/kehuyuwHvBg5KCDDlInLvmGrJOT61FHHYXi4uJWq9aK5uZmdfF9MwMHDgwvGJG3Y2/d3xhSzW6gsdLdcbboIqCxHEjIAn73unvydUIGkOrux2xXbELLiK4uCkbWrl2rVrq97777cPrpp2Pfvn244YYbVFfZyy+/3KFg5OCDD1YZKQlEZGVfCW62bNmiulz0/cg+5Tmdff3hBJ+SKZOVfeWie/TRR9Xr9V0V2BeDEWrPt3u+xTWfXcNfFFEnGWDAmMwxeOv0t9T5oNeVgw98U3r8E+rNyrf++fPnd+6HSSDy13zsFwlIXprR8efdWQyY3dmKcHz00UdISnKnUX27Qnw9/PDDuPDCCzF37lx1XwKJxx9/HMcddxwWLlzYoWzQbbfdpgIaIb/fsWPHqmBEghQxatQo9UezP0pLS9Gvn//ALrnvcDhQXl6ugs9QbWQ7UWdIt8xNX96kPky1oPMK9s/xA45HWVUCfi7djtjklqyfbnDcYThqyEj1jVTS+XmZTajQVgfdT11DMr7ftSXofoYnTMXhg0aoWTLf7NgctM2IxKk4bOBwlQn4atumNttId9EXW0K0SZqKwwYMR2HZTvxS2ZL10NnrRuOIQcORnFiPZbtbVjfWpbgmorw6AWMGObGl/seg73XdLiNKGkqC/vzjBhyP/KQ8LF2/V9VqGTPQiS0NHd/PMfnHYUBKPj5evx61xp+Dvo/fHTJBldRfvH49akK0mT1pvOr+WrxhPWoMwducN2Ecau37gv4+nHWj4bSn4ejRJvy4tyULpTu2/3H4stAOl6kq6Ps4fsDx6JeYi7dX/QIkrm/1+OD4w3Dk4BHq+/bbK3bBZqgMuR8R7DXKe5iYOwQ56U1BH5f/O5JZLCgu6NLsSERM7c3NzW11kikrK1PfxjMz/QfT6ObNm6eiKP1SVFSE3kgyIGvWrPG7vPDCC35tVq5ciVdeeUUFLfrl1FNPVf1/khHqiAkTJnhv6xkpORa6X3/9VXUBdUfwGaxNV0Xi1LfI385ff/grbE5bhwORcPLCBhixr2kfzh54LYwxdYAW+LdrgN1QjTum3oGLRt6M5r2zUO+o9E5pDdzPFaPmtrGfGrWf6yfc1mabeVPn4dbJf1JttCBtHEZ3mzsOmxeyjdNQizsOuwPNWnWrx+XnyvNunHCbes2B70Xu19kr1XuV/bR63PNeK4tmqv3I/cCfv6dur/r5lbtOR/PeM+Eyht5P+a7Q+ymqK8WfptyOWntlyPdxQtb/4PZD70BNG22OTr/K3cYWus1hyVeo1xN4XOR1JyQ0qt/H7pq9QX9fFdZyTEy8VO1H/haCvc87p85T+wl2vOR5ckzvOuJOHJk2J+R+yhrL3K8x8PTu2cdtU25Xj0vgHoy81idWP9HtY0cOeGZExiR8+OGHfts+/fRTHHrooSHHi8gMC7l0inSVSIaiI0p/CZ4JuWIJkDuhYz+7A2TMx4gRI/y27d692+++BB1XX321GgQaSLo7hJzEA/9w7PbWi2H5/r71E7/sv7uDz1BtArMlROF2z4QzPiSYcOJfDS71bXFYwrswxe8Osg8NxU2b1bfJAenjYErcjGbTzpZFPwL280PWv0LuZ1fDRrWfBocjZJud9e42MqU3VJsdde42Lk0L2WZ73a94Yd0L2Fi1vvXvweB+3neV/1KvOZBLcwGWIliylqn9tHrc816t8UsRF78briA/f0vNBjWDR6ZYW5K3YFtb+0kIvR95r39fsVC9nlaH0/M+Ptq8DGWO9HbbVGmZbbZ5c+PL2NBc2DIf1uf3Icc8NnMZiq2bg/6+5H1MyfwKptrdrUJml+d9yvFoNOxo9fP146VnLPLydsG0J/h+1le2zqr4voeVNf8OekwDX2t3Z0c6HIzU19er1L5Ovp3LN/qMjAx1cpSsxp49e/Dqq696By/KNOBbb70VV111lRrE+OKLL4YcqLjf5H9VB7pKlBj3CHh3JOlquZbtHd1XF5s8eTIKCwtbBS2+srOzUVJS4r2/efNmNT6jJ4QTfEqbpUuX+o0ZkTYyNoaoIyQIf/DHB7vll/bpnjfVN1Q5MQSSb5nybfLp41+BJfvTkO3E/215pc39PL7qcfeXBfl2fADbiOd/eT5015YG/HP9iyEfl+8/5qxlbXaNWbK+bOPnG/Ds2qdkBBqS8j6DvdP7Af5vc+jfqWz/qvxV/NxgbrNNQdVr2LDO0kYbYIP1PU8E27nXubb+P23+bTzfxvHQ/8am5U3DmvpFbe6nLS+u88++B6P/rCPzj+y2jHWHu2l++uknNUBRLkKCDLn9l7/8Rd2Xk+KuXbu87WXg5OLFi9XAzEmTJqmBmDLm4bzzzkPESMwGknKA/InAGY+6r+W+bO9hMjtFArjrr79eBX0SaMhslRtvvNHb5sQTT1QB36pVq9TxkQAwnFlKgWTsSHv1PtavX69eR2VlpepC07uXdPKzZRCq/F3IjJqXXnpJBZ8yXkV38803q+DjoYceUl1Dcv3ZZ595x8UQdUdWpKNsruaQH/5y8ihtKEVsrB3G2Oo2TxLNzvb3I+MjcIDb6K8lZNeWoe3H1TnKYG+7a8zgaOPna6iwlgFGGwwx1fuxH8Cu2UL+TmV7k7MCFda9bbaxuipR0dxWG3nFknHej9fpCv06RVu/b/2YNjoaUdVc1qlARP8Z7dF/lt3VOsMeMZmR448/vs2+JBnfEEgGW8qJMmKl9gfmrgNMZvdf3JTLAacNiOlkV1EXknEeUkRMZr9I7Q753Q8fPlxNj9X94x//wOWXX65mskhlXJnJImNNOmrjxo0qwGiL1BDxnfGiB6X634QefErW46mnnlKvJzD4lAzIokWL8L//+7/485//rN7P22+/jcMPP7zDr5n6Lvmbe3z14+22S4pJwjUTr4HVZVW3U+NS8c2Wnfj3qj0YmZWD608cATgScPPbq2EwNeGR305CQmwS/ueVTer5//c/RyDBEgOHy4HfP7cCDXYnFv5+CprtTsz9vzXISIjFq1cejoy4DCRbkpFbfwe2Ve3F/DPHwO7QcP/iDRiZk4RHz5+k9if7ue/DjVi5qwo3nuDOeD7x5RZMGZSO+WeNVfsRf/noe3y2oQwXHT4IMUYjXvluB44/KBu3nTrK2+Yfn6/Ef1btwRnj8xBvNuFfK3erqaXXnzDC2+b5b3/BywU7VO2JtLhYNUX2vEP64/Kjh6rXEmOMUYNIH/t8Myb0T0VOikX93AunDsLsw/LV4+KHbZW47+P1akrpgPQEfLOlHJcdORC/mdJS+VMKnP3pP78gM9GMYdlJWLGjHNceNwKnT2iZVCCVZq97cxXiY4wYmJ2PfS4N/zPiMRw9Ss9QAxX1Nlz68o/q2/LgrERsL6/DbaccjOM9Je5FY7MD5z/3vTp1ZyXGoLzBgTtmHOypseHmcGqY/WwBmpsTYYk1wKbV4a7TRmOaT/Ev+Tu68PnvUdsQD6tBg8vUgLvPGIPDAkqjX/nKCuyta/J+h7//zLGY5Kn0Kgq2luOh/67HUK0EV8T8Fy85ZmJf3BBcf+JwHJlcDvywEKvzZ+Pe74HBrt1B2zh+eBrr+1+Ie753YpDTp038ELVw4RmZNiS9+TvcNuha3FZQi0GOljbl8UNw3QnDcZB5KwxrXsOe/hfh3h9cGOAo8mtzzXFDcPTeRSj/9V18pE3CO65jca7xaxxr+gVfOyfgHfsMaM4k3HXmETDLObGbdMtsmqjgG3hIQHKAA5FgQVuoYO+www5TmYRQ5IT/ySef+G3zreA6ZMiQVvuUqb+B28IZsCTThNsTTvD5m9/8Rl2IOkv6tDdUtp5NoLtz6p2YmDNRnZRzE931IXSVe3fAUV2IIQNyccawKarA102NTZBK4Uf0OwmNNidc1gZV0GtK3nhvqjo3vgqb6+qRgEHQXNJmHwZkp6npkLohafnYUhIDg20gNIcLLmstRqTm+bUZnWnEik070dyYr76FuqxNGJs1xK/NxJxx+HT1RtTV5qqCYC6rHRNzDsKYzJHeNoflT8C/CoCKqkwkmGPgsmo4LG8MxmS2rK1y9CAjXvzCjtKyJDQnmOGyGnHEgIkYk9lSssAxNB+PWhuxq9QMe1McXNZYHD1oEiZkt/zeUoyNmG+tQVGpAc0N0saCE4cehjGZLSf/QUkO/LG5AvusgLU+Bq7meJw4bArGZLactEela7A4y9FgdWJNg/sz5+ghIzEm02faZyaQGbMP++pt2KqGvKThlBFTMDJTCrS1GJZSic1lDSjzVKGfedChqiicbsm6Eris/THOtRXzDG/iAceFeHZpM/onpGNGegmw9C8oGHYTrA35GOfcgnkxb+KB5gvx9KdWzD/Tv011dQ7GOjxtHBfisSVNmH/mMNWm4r078MHuWbBrw3BBzBf4HTaiyTEE91ZPxf97pw7Lxn2OwTt+xA9b0tHsuBQXxCwN0WYFVm7JgDWwTZW7zSnjFgM7lmPfloSQbZaN+wmDd6zEy1uy0KjafKraGB398Fp1Fj54vxCXpnyKMc3NGK+tg9U+BXfGrkOaoQlHaBtQaDsTRjTj7SWF+O2k8d1WIp7BCBFFHQmcpU+7rf7197e+j98d/Lugfd76ir0pce7uTPnAzUyyqMXyZKVeff2ZnGSL3/Pz0uJVMTIpCV/vWbcmP81/er1kDfTiWfp6KgMyWr71i0GeE+auygboY8gHZ/oPgB+W5anCWt6gghGhV2bV6Yv8yWJ4iWb3x/mIHP8Ttl6FVYqx7TW7U/SjcgPa5CSr72Cyzk51o3utndF5/m2kEqssdifr8+ypbgraJskSg8EZCarSp1RelVctJfB9ye9anrdqV7UqeR9rNHiryarq156Tf63VifGGbd6T/8UvxqkF4XwDhN1VLr82v302DvfMGuMNEJ7ePQs2bRjOjVmOI03rca62HPfWDcO1r6/CsnGLMXjHcmzckqACBL82ta3bNIVo8+XYxRiy7wfMMRnwrHMWZpncxR3PMn2LNa7h6neQuuV9te1s0zdY7RqOc0zfhN3mbNM32OrKRSxcSN/i7kb/jekrVGmJ+K3JXXzvPNNXcGoGmOBC9hZ3IbPzTV8iSWvEGSZ38bTZMV+ri9A8PTXphgY8bpYqwm6ZqMXHlrvcd2zAd9un+2WRDiQGI0QUdaQvW/q02+tfl3bBUs36ir2pCS1jq7I9wci++mZvoJGT4h9o5Hnul9ZYVWlxkZ/qH2jIGjVCanxIaXHf4EM32FsSvlFlZcQQzzadHnhs31fv/XaqByg6fZG/vbXNMBqa/YIP3yBC1rppsDnVzBUJOvQgRiddPHoQIS9H1j8JfF8SlEkQo5dFT0+IRVqCuVUQUVrrH0RMf/RrtQibbxDxa6l/EHHC35e52+xapL75dzpAqOmKAOEbbHT1hxkOv5N/jZbgPfn/xrQMMZoTsXAgZ6t732fGfK8ueoI5w1CPx3xO9MFO/uG0STc04H7zP/3aJBusuDX2He/9FIMVc2KX+LVJMNjw29iWyrx+x7Kd2WV2zYTb7FfjxLrWCx8eKAxGiCjqSIDxh0P/gDu/ubNVt4xOumdC9XnXNOqZkZaPwOxkC1ACFZDomZN+gcGIJwtSUtOE2iY9MxIqGGlCo83dZqAnW6LTsyAyfkJfITYwM6IHJ3rgJLwZBA9ZO0Zet7xm2Y0s1CfZnMAgYnhOEn7xLEYnj8ea/OcuSHeGZHt8A4RjHza3CiIKi32CiKYLcfRDplZBhDVEgOAbRDT6tLlY+xSv1QJPvrEBJ6T8C5YQGYR1rsGIg90vQKjXLN4A4bemZUjUGmGBA/22rtyPAKEBD5pfanXynxvbMrg/2dCMi2M/82vT8vtG1Dvbdi8KtaH4XfJ+LrHSAQxGiCgqu2ne2PBG2N0ygfRVUFPifTIjnpO4bzASeGKXdViEnLirPQFN62BE76ZpREOzOzPiO47BN1MiXRlCEh/683RSDVSyGnqXiPxsGRcSaFhWgnrN6vWmWFRQohaI9clWbN7rE0TUSRBh8AYaendGc2B3RhcEEXJdrGUgDjZkbnGv1TLbEzTMMn0fdveBBBH/MD/XKkC4MfYD7/0kQzNmh8oE9IIAoSu5pACaDNbV3H97Lffd1wbP39vUgEG8BxKDESKKysGrvoWbOlrG2huMeMaMBAYj+riJfimBwYg78CiptqLK00YChmCZkfJ6m/dEGDiuRAIN2bd0r7ifk6DKlQcanNkSjMgMFenSUV02PoHG2j0+gUa5BBrOVtmKznRnrHMNQSKsyNji7g6QDESsZsd5nkAjnCAizdCAu2L9a0olthU0IPoFntjVNlUrNfTJvzvbVGpJKDXmwjnxIphWvYRRKMJGbSBed5yixpnkoQIVWor6G+quwauCwQgRReXgVSlbrSqBBpSxDqdQk97FIt0cvmNGhIwZqfQEEq26aTyZEcl6NHoGp+pdNzrZpwzkrPdkPXJT4mCJMbV6DYPSE7zBSGp8jDvQKFmtggycci+WVOVhTVFNS6BRIoGGrVWgEU624mfXMJxrWu4d7AjNhRRYkbv1h5DdGf8wP+v3epMMVlwU+0VUBhHdefLfouXjFccMdWLvDym7DpgyBmPP0N/A8FPrk393tslDBX7TfDfu/P1MzBifjyXDz8eNH6zE1jpVNAZvOk/EoJQY3P3biZgxrvUitgcSgxEiiuqsSGfKWMtATpES3/IRKF0cYl9tMyoamv2yJb6zaYQMBhWSzZCMhS8JhCQ78mtpXdAuGn2Mxrri2pZAo9gdaCwa+I7qEtnx5Uu4dt1MNTy3vUBjp5btHTcRKluxwLzQb7Dj5bFLI7I7o73ug0gNENKOuhK1374IZ1URLrHdjlJkeU/sd50+GqdOGIQ0gwFLhgY/+d/VjW3u9Ak0JCA5ZWweftxeibI6K3KS3V0z3ZkR0TEYIaJeNaU3nOxI0G4an8xIuWcMRmBmRDIeMkhUH+shXTTBfk7/tDhvMBIXY/TLekjXyrWL7X6BxlXaR3iu1oDEzR+oNEPalvdwpjEJSWjCeaav2ww07ol9Peh7NERZENFe90EkBwgSxaVNuRxOezMeLWoIeWIP5+Q/oxvbCLnfXdN328JghIj6zJRe1UbTvANU/bppPFmQPVVNsDldQYMRvVumbm+9uh04FkTPeny3rbIl67FFsh713qxHxa4YjDWcBBkmqI/RCOwmSUMdHvPJZnRXoNFTQUS43QeRHCDIzzOZ4zBteNszUMI5+Zu6sU2kYDBCRFFDAoxFZyxCpbUSr6x7Bf/d8V/MGjYLF425KKwpvXoXi2c2bdDZNHogItVXJRMSKDc1Hpv0YESvxeEzoDR41gNI3iyLrAGzXF9iluVL9bTAosdd2U3SmWxFR7szujKI6Gj3QSQGCNR5DEZ60GWXXYZ//tNdzMZkMqmy7qeffjr++te/Ij093VvKXRaQ0xeRk291sujcc889h/fff18tkufr6quvVo89+uijfgvPSZl5WePGl6xvI2vEEEUTKe0ulzq7uxtkcr/JfmXU26OPFzGbjKq4l04Cj7hYI6z20FkR9fNTLN6sx/eNN8PpmgDTz+4BpZL1GGc4Ef0N+1QdDN+sRzCdDT7aCzQ6kq3Y3+6MAxJEMEDocxiM+Piu+Du1HPkdU+/AtPxp3XIAZsyYgZdffhkOh0OtiHvFFVeodWXeest/OpxwOp246qqr8OGHH+KLL75Qa9b4eu+99/DDDz+ooCYYee69997rvR8f7z8lkSia7Kx1L9g4OKVlobZweEvBx8f4jfeQ25IdKaps8q8x4sl66DNcFq8txR9M7qxH+eaXcf1f9+EJw5uQHItv1qMzvEGFJ2Pht60DUzM7kq3oyu4MwSwDdQaDEQ/JODy26jFsq9mmro/IOyKs4kn7y2KxIDfXvRjVgAEDVLYi2CJ6zc3NuOCCC7BixQp8/fXXGD16tN/je/bswQ033KAWzJPsSjAJCQnen0UUzexOO/bU71G3h6QM6WQw0tJF4zuI1RuM6JkRT9aj6NPH8dSmKRgBJ841L2/JejiCZz06k+GQbpJ3Y8/A1UnfwlG1U2UvdmvZeNtxQocCjQ7PmOjC7gyizojpjUFFk8P9YdIR3xd/750uKNdf7voSR+Qf0aF9xMcEH1kfrm3btmHJkiWIjfX/kKyvr1cBRlFREb799lsMGjTI73GXy4WLL74Yf/zjHzF27NiQ+3/jjTfw+uuvo1+/fpg5cybuvvtuJCf7L3RFFA0+3PqhmsprMVqQFd+yamw49PLqvjNpdJIZkS6Ye2JeQVHzhUCxASh0F/0auPNdfGhpKQne2cBD88wFkgxHkZbjE2iU41Lb7bh79nSkjX0In/yyC//v4w3YVevY76mZDCIo0vW6YEQCkcPfPHy/93Pzsps7/JwfLvwBCbGtawq05aOPPkJSUpLqgrFa3YsSPfLII35t7rvvPhU0SDdOTk5Oq3089NBDiImJwU033RTy5/z+97/H0KFDVWZk3bp1mDdvHn7++WcsXRq83gBRJH/heGHdC+p2Z4L/lhojrYORzCQzjjAtxxTTFkzZcS/wXEu3ZldlPf6U8yOaK4owxzUfq+rSQhabOnXiYJw8fhADDeoTel0wEm1OOOEELFy4EI2NjXjhhRewadMm3HjjjX5tpk+fjs8++0wNbF2wYIHfYytXrsRjjz2GVatWtfnBLONFdOPGjcPIkSNx6KGHqudNnjz5ALwzogNDipoV1RWp21anNewS8Dq/ab36eJAjrkPBvlhsXrMZt3sqlcpMl3BjncCxHm1lPQxjc2Fx2vAvo5mBBlFvDUakq0QyFB35lnX5J5djY9XGVqWlR6WPwsunvhz2ty/52R2VmJiIESNGqNuPP/64Ck7mz5+vsiG6k046SWU9zjrrLJVBeeKJJ7yPLV++HGVlZX5dN9LmD3/4gwpcduzYEfTnSgAi3UGbN29mMEJRW/Qs3CJnwQuexQA/v6bGg8jlSABH+g+j6FDWQ6bE3pm7ArH1xe1mPRBjgRSI5/gLol4ajMgHUke6Sr7d8y02VG5otV0CE9m+Zt+aDn3r2l8yjkPGc1x77bV+s2JOOeUU1aUza9YsNUbkySefVO9VxoqcfPLJfvs49dRT1fbLL7885M8pLCyE3W5HXl73rj9A1JML5Ekm5LxVN2CH8ViMdB0BrPu32ixfQ1ovUwe/DEl7WQ+kDsB91/0dJs0eVtaDiHpxMNITpaW7ktQDkUGo0iUjAYcvqSny8ccf44wzzlCv/amnnkJmZqa6+JKMh4wNGTVqlLq/detWNXj1tNNOQ1ZWlhp7IpmTQw45BEcd1X2BFlFPL5AnM2MGNhZigbkQWNtS4TRUICI2a/1V1uOm9O9grNujZrPsQj9v1sMCB2yIxUJZ5dQke2LWg6ij2vo/2Ot1pLR0d7r11lvx/PPPq9kzwYKVxYsX47XXXlPZE/mAbo/ZbMbnn3+uMiYSoEiXjz4ORYqtEUVTVsQ3EAlcIC+o6l1A8WqgeI13Zkxb/2306qzSFSNusV+LN50n4YeT/oXV5y6HPVWmEutBjwEZqclYeNHkbl/llKg3MWjhnM16WG1tLVJTU1FTU4OUlBS/x2QGyvbt29VMkbi4tufIByPBhpSWDkVKS0u1R+pZ+3ucKbrJx9QFH1+A9RXrQ2YxpQrrW6e/1To7ck9q2D/H2wXj9HTBGCpwZvP9KEUm3rrqCDXGQxa9YxcM0f6fv3316W4a39LSRNTLFsjTZ8qccCfw1d8Al7u+SDD6eJDLbX/Ez9oIbxeMGQ7YEYu8VPe4D8GaHURdr88HI0QUXQvk7W3Yi5u+vAlmoxmvnfZa6AXyPJVTZWApUvLd3TWBskcBh1+j1maxVRahTMvw64KRQETcLeNBOACV6IBhMEJEUZXFdHgyHJnxma0XyJOAo7HCHVB4xoeooMSbUZFAw2fl2rOfBfofotZmkYqnWLwFqLG2/MzUOBWIcDwI0YHFYISIokp1c7W6To9zr2ztZ8H4IM/w7drRcKf9SjUeZHxiPYxJnorGBkPYFU+JqOsxGCGiqKIPOE+3BAlGzn0eeO/a4ONDDCY0zHwCb76TpsaDbL75JBgt/oUKOR6EqGf0mqm9UTApiPYDjy/pqqxVoTMjE2YDl30U/Jd11ReoGH6OuplgjkFsQCBCRD0n6jMj+gq3srZLfDw/XHorOb4icEVj6nuqmkMEI/rsmVZZEfnO5QooBc+/I6JIEvXBiBTtSktLU+uziISEhG6rlkrdkxGRQESOrxxnFmkjb2YksJtGnz2jyxwJTLsOWPUqULsHSMxGzT59xd6o/+gj6lV6xf9IKX0u9ICEeh8JRPTjTH2bXzeN7+yZtf/X0mjK5cDkS4DELPdtp00tTle7q0Q9zMwIUWTpFcGIZEJkwbecnBy1+Bv1LtI1w4wIBe2mCTp7BsDKl90XcU+NCkR8u2lS49lNQxRJekUwopMTFk9aRH2om6at2TPGGODslsXwRG2Tu10KgxGiiNKrghEi6mPdNDJ7Jusg4LnjWjec8zmQP8lvU02TPoCVH31EkaTXTO0lot7P7rSj3l7vLf+uNOwL+2ON3TREkYnBCBFF3XgRk8GEZHOye+POb93X5kTgjEeB/ImAVFZNzG71/Fo9M8JuGqKIwlwlEUVdF02qJRVGg1HmfgO/fux+8NQHgSmX+M2eCVRr9YwZYZ0RoojCzAgRRV1mxNtFs2clUL4JiE0Axp3jXWcmWCDiN2aEdUaIIgqDESKKGj+W/KiujfpH15o33NejzwQsnm6bEJwuDaU1Ter27qomdZ+IIgODESKKmmq8H277UN0ubSyFZmsC1v7H/eCkC9t87pJ1JTj6oS+wp9qq7t//8QZ1X7YTUc9jMEJEUaGguAClDaXqdq2tFgUvHQM01wCpA4Ehx4R8ngQc176+CiU17kBEV1pjVdsZkBD1PAYjRBQVWZEnVj8Bg5R9dxd/xxPGWqiOlokXAMbgH2XSFTP/w/XudoH79FzL4+yyIepZDEaIKCqyIoUVhdA8IYT8W2ixoCA+zl3YrHi1e52aAD9ur2yVEfEl+5HHpR0R9RxO7SWiqMiKGDUNLp8VueX+E+mpOHLRhZ58iWcdGh9ldaEDkc60I6IDg5kRIoqKrIhvICLkvjc7IuvQyDo1AXKS48L6GeG2I6IDg8EIEUXNWJFABk92RLvyM/c6NQGmDs1AXmpciGe7x57I49KOiHoOgxEiilh2l13NoNHHigTSDAaUmmJg14Ks2itl440G3D1rTNDH9ABFHpd2RNRzOGaEiCKW2WTGojMWodLqHmB6+ZLL0OhowqNllci3NwNpQ5BhrYE5OT/kPmaMy8PCiybjjv+sRbWnAqvITY1TgYg8TkQ9i8EIEUW03MRcdZEuG6uzWW2baNeQbbMD578NpA4IWf5dJwHHvvpm/Pm9Qozvn4I7TxujumaYESGKDAxGiCgqWJ1WuDSXup1ob3RvjE9vNxDRNTY71fXInGRMG5554F4oEXUYx4wQUVRosDd416WJl9V6RVxq+M9vdo8rSbTwOxhRpGEwQkRRod5Wr64TYzyzYyypgNEU/vM9mREGI0SRh8EIEUVVZiTRFNfhrIhvZiTJEn4AQ0QRHIw8/fTTGDp0KOLi4jBlyhQsX768zfZvvPEGJk6ciISEBOTl5eHyyy9HRUVFZ18zEfVB9XZ3ZiTJaHZviO9YMFJvYzcNUa8JRt5++23MnTsXd911F1avXo1jjjkGM2fOxK5drdeFEN988w0uueQSXHnllSgsLMS//vUvrFixAnPmzOmK109EfSwYSZRqqyIurUPP55gRol4UjDzyyCMqsJBgYvTo0ViwYAEGDhyIhQsXBm3//fffY8iQIbjppptUNuXoo4/G1VdfjZ9++qkrXj8R9bFumiT9Yyu+c8FIEgewEkV3MGKz2bBy5UpMnz7db7vcLygoCPqcI488Ert378bixYtVnYC9e/fi3//+N04//fSQP6e5uRm1tbV+FyLq27wDWPXaqR3MjHAAK1EvCUbKy8vhdDrRr18/v+1yv7S0NGQwImNGzj//fJjNZuTm5iItLQ1PPPFEyJ/zwAMPIDU11XuRzAsR9W3ezIhL28/MCAewEvWKAayGgNUzJeMRuE23fv161UXzl7/8RWVVlixZgu3bt+Oaa64Juf958+ahpqbGeykqKurMyySi3jhmxOXcr9k0nNpLFHk6VP0nKysLJpOpVRakrKysVbbEN8tx1FFH4Y9//KO6P2HCBCQmJqqBr/fff7+aXRPIYrGoCxFRq8yIw9HJbhqOGSHqFZkR6WaRqbxLly712y73pTsmmMbGRhiN/j9GAho9o0JE1KHMiMPWUgo+TA6nC80Odyl5DmAl6gXdNLfeeiteeOEFvPTSS9iwYQNuueUWNa1X73aRLhaZyqubNWsW3nnnHTXbZtu2bfj2229Vt83UqVORnx96pU0iIl8NNk9mxG7tcGakwVN9VbCbhijydHiRBhmIKgXL7r33XpSUlGDcuHFqpszgwYPV47LNt+bIZZddhrq6Ojz55JP4wx/+oAavnnjiiXjooYe69p0QUd/IjNiaOjyAVS94Zo4xItbEwtNEkaZTK0Zdd9116hLMK6+80mrbjTfeqC5ERPs9ZqS5vsOZkXorx4sQRTJ+RSCi6CoH7w1GUjs8eDWR03qJIhKDESKKroXyOlFnxDut19ypZDARHWAMRogoqiqwJrlcQGwiYIoN+7ksBU8U2RiMEFHEszvtsLncU3oTNVeHq6+2dNMwM0IUiRiMEFHUdNF4u2k6uWIva4wQRSYGI0QUNYNX441mmDpTCt7mrjPCAaxEkYnBCBFFz7Reo2ecCLtpiHoVBiNEFD0FzwyeMR/spiHqVRiMEFH0ZEb0jyxmRoh6FQYjRBQ103oT9bU1O5kZ4WwaosjEYISIoqf6qtQY6URmRF8oL4kVWIkiEoMRIoqi6que1Xc7OJvGW2eEFViJIhKDESKKngGsDrt7AwewEvUqDEaIKOJtqd6irmtczZ3spvEUPYtjBVaiSMRghIgimqZpWLV3lbq90uSC1onMSB0HsBJFNAYjRBTRCooLUN1crW7vNRlREB/XocyIBDMsB08U2RiMEFHEkkDiidVPeO8b5H56KjRL+ANYrXYXZDkbwam9RJGJwQgRRXRWpLCi0HtfMxhQaLGgYJ+726YjM2lEQqxa2YaIIgyDESKK6KyI0eD/MWX0bJfHO1TwzGyC0Wg4IK+ViPYPgxEiiuisiEvzFDrzcEl2pKJQPd6hGiMWzqQhilQMRogoYrMiBgTPZMj2cLMjHLxKFPkYjBBRxLG77ChtKIXmnsjbimyXx6VdexpszIwQRTrmLYko4phNZiw6YxEqrZWwO+246L8Xqe2vFJcifvgpwEn/i4y4DNWuPfWedWkSuS4NUcRiMEJEESk3MVddqqxV3m2Tmm0wJfUHMseEvR920xBFPnbTEFFULJIXDyNM+7EuDQewEkUuBiNEFNEaHY3eYMR9o2PBCGfTEEU+BiNEFNEa7e5gJFGfWcMVe4l6HQYjRBQVwUiCXtO9w5kRzwBWM4fIEUUqBiNEFNEaHO4xI4kud1CBuNROjhlhKXiiSMVghIiiYwCr09Gpbhp9zEgSK7ASRSwGI0QUHWNGHDb3Bg5gJep1GIwQUVTMpmnppunc1F5mRogiF4MRIoqeAazGGMCc2LlgJI4DWIkiFYMRIoqKMSMJsnqvDF41BF88LxTOpiGKfAxGiCgqumlUZqSDXTSC3TREkY/BCBFFRWYk0eXq8OBVp0tDk50L5RFFOgYjRBQdmRGt45mRBptnOjDXpiGKaAxGiCg6pvZ2IjOid9HEGA2wxPDjjihS8X8nEUXJAFZtv1bsNXRw4CsRdR8GI0QUJVN7PbNpOjGThjVGiCIbgxEiip7ZNJ3spuG6NESRjcEIEUXHmBFVZ6Rz69JINw0RRS4GI0QUsRwuB6xO635nRthNQxTZGIwQUcRqcjR5b6vZNJ0dwGpmZoQokjEYIaKIn0kTo2kwy40OZkbq2E1DFBUYjBBRdAxeFTW7O9lNY+r6F0dEXYbBCBFFx+BVsXlph57f4JnaywGsRJGNwQgRRabqXWgo/bmlxoj49UOgeA1QvFo93h7OpiGKDhzVRUSRacF4NCbEA/2ykah30zRUAM8d19Lmnpo2d8HZNETRgZkRIopM5z6PBpP7+1K8lIJXPNfGGPV4e5gZIYoOzIwQUWSaMBuN1t3Ahpfd03p9zfkcyJ/U7i44gJUoOjAzQkQRq9FTZ0QtkteJjyx9AGuSJbbLXxsRdR1mRogoYjUa3VNyVWak33hAum1q9wCJ2WE9v6WbhlN7iXpdZuTpp5/G0KFDERcXhylTpmD58uVttm9ubsZdd92FwYMHw2KxYPjw4XjppZc6+5qJqI9oMJla6ozkjgeu+hKYuw5I7R/e820sB0/UKzMjb7/9NubOnasCkqOOOgrPPvssZs6cifXr12PQoEFBnzN79mzs3bsXL774IkaMGIGysjI4HO4PCSKidoueqUXyUgGDAYixdGLVXiaBiSJZh/+HPvLII7jyyisxZ84cdX/BggX45JNPsHDhQjzwwAOt2i9ZsgRfffUVtm3bhoyMDLVtyJAhXfHaiaiPlINXmZG4lA49t9nhhN3pHmvCYISoF3XT2Gw2rFy5EtOnT/fbLvcLCgqCPueDDz7AoYceir/97W/o378/DjroINx2221oampZACtYt05tba3fhYj6cAVWtUheaqcGr6rnmzlmhKjXZEbKy8vhdDrRr18/v+1yv7S0NOhzJCPyzTffqPEl7777rtrHddddh8rKypDjRiTDMn/+/I68NCLqxd00iTKbxtKxzIjeRRMXa0SMiRMHiSJZp/6HGqTf1oemaa226Vwul3rsjTfewNSpU3Haaaeprp5XXnklZHZk3rx5qKmp8V6Kioo68zKJqNd003Q8M6LPpEnieBGi3pUZycrKgslkapUFkQGpgdkSXV5enuqeSU1t+SAZPXq0CmB2796NkSNHtnqOzLiRCxH1bXo3jaoz0sExI6y+StRLMyNms1lN5V261H/lTLl/5JFHBn2OzLgpLi5GfX29d9umTZtgNBoxYMCAzr5uIupLwch+ZEYSzZxJQ9TrumluvfVWvPDCC2q8x4YNG3DLLbdg165duOaaa7xdLJdccom3/YUXXojMzExcfvnlavrv119/jT/+8Y+44oorEB8f37Xvhoh6Fe+YEVfnx4ywm4Yo8nX4K8P555+PiooK3HvvvSgpKcG4ceOwePFiVdBMyDYJTnRJSUkqc3LjjTeqWTUSmEjdkfvvv79r3wkR9SrSletfZyStQ89vqTHCmTREka5T+UuZDSOXYGRgaqCDDz64VdcOEVFbmhxNcEkQomdGOjBmxOnSsG6PuyRAk92p7puMwQfZE1HP43w3IopIelZExMXEA6bwFrtbsq4ERz/0BV77fqe6//22SnVfthNRZGIwQkQRP3jVGOZ4EQk4rn19FUpqrH7bS2usajsDEqLIxGCEiCKSd7yI6qJpfyaNdMXM/3A93AXg/enb5HFpR0SRhcEIEUV0wbNENXi1/czIj9srW2VEfEkIIo9LOyKKLAxGiCjCu2nCy4yU1YUORDrTjoi6D4MRIopIDQ6fUvBhjBnJSY4La7/htiOi7sNghIgiUpO9qWWRvDAyI1OHZiAvNQ6hJvDKdnlc2hFRZGEwQkQRqbC8UF03ShQRxpgRqSNy96wx6nZgQKLfl8dZb4Qo8jAYIaKIrL76ZdGX6vZWsxlamFN7Z4zLw8KLJqNfiv9Cm7mpcWq7PE5EkYcrSBFRxCkoLkBZU5m6XWUyocBRjaPCfK4EHBMGpOHIB79QGZE35hyOw4dlMiNCFMGYGSGiiMuKPLH6CRg8nSsGuV++Qm0PV02TXV1nJJpx5IgsBiJEEY7BCBFFXFaksKIQmqdUmWYwoLCpWG0PV1WDTV2nJ5oP2Oskoq7DYISIIoaeFTEa/D+ajDCo7eFmRyob3cFIRgKDEaJowGCEiCIuK6Kv1qtzQVPbw82OtGRGwltcj4h6FoMRIorIsSKBDB3IjlQ2tIwZIaLIx2CEiCKC3WVHaUOpd6xIINkuj0u79lR5umnS2U1DFBU4tZeIIoLZZMaiMxah0upeyO6y/16KJqcVj+4tR/4VnwEGAzLiMlS79lR6ummYGSGKDgxGiChi5CbmqovT5VSBiJgEC7KyxnZoP8yMEEUXdtMQUcQukieSzckdfj4zI0TRhcEIEUWcOluduja7NFjCWCQvEOuMEEUXBiNEFHHqbfXqOtnlAuLSOvx81hkhii4MRogo4tTaaluCkTAXydM12Zyw2t11SlhnhCg6MBghogjPjKR2KisSazIgycIx+kTRgMEIEUWcers7GEmSSqxxKZ0bL5JghsEQvIAaEUUWBiNEFLndNM5OZEZYY4Qo6jAYIaLI7aaR0u8dHDPCGiNE0YfBCBFF7NReZkaI+gYGI0TUO8eMcMVeoqjBYISIeteYEc9smgwukkcUNRiMEFEEjxnpeJ2Rqkb3qr7pie0vqEdEkYHBCBFF7JiRJJfW4cyI3k3DFXuJogeDESKKOPX6ANZOlIOv9KkzQkTRgcEIEUV2OfiODmDVx4ywm4YoajAYIaKInU2TbLQAptiwn6dpGqoaOGaEKNowGCGiiNLsbIbN5Q4okmKTOvTcBpsTNpmBw9k0RFGFwQgRReTgVYOmIamjM2k840XiYo2IN5sOyOsjoq7HYISIInJab6KmwdjZdWk4eJUoqjAYIaLILAXv6nzBM9YYIYouDEaIKKLU2fUaI50oeMYaI0RRicEIEfWezAhrjBBFJQYjRBSZpeBV9VXWGCHqCxiMEFGEloLvTGbEU2OEA1iJogqDESKKyDEjyfs1ZiT8QmlE1PMYjBBRBI8Z6eC6NJxNQxSVGIwQUYSOGenEujSsM0IUlRiMEFHkTu3t4JgRfZE81hkhii4MRogoQrtptA6NGXG5NFQ1ugewcsVeoujCYISIIribJvzMSJ3VAacEMADSEjiAlSiaMBghoohSZ6vt1JgRffBqkiUGlhgukkcUTRiMEFFEqbJWqestZjNgTup49VVO6yWKOgxGiChiOF1ONDmt6vabqalwd7qEhzNpiKIXgxEiihhfFn3pvb0x1oSC4oKwn8saI0R9LBh5+umnMXToUMTFxWHKlClYvnx5WM/79ttvERMTg0mTJnXmxxJRL6ZpGhb+vNB736gBT6x+Qm0PBzMjRH0oGHn77bcxd+5c3HXXXVi9ejWOOeYYzJw5E7t27WrzeTU1Nbjkkktw0kkn7c/rJaJeSrIgm6o2ee+7DEBhRWHY2RFmRoj6UDDyyCOP4Morr8ScOXMwevRoLFiwAAMHDsTChS3faIK5+uqrceGFF2LatGn783qJqBeS7IdkQQww+G03GoxhZ0f0zEg6p/US9e5gxGazYeXKlZg+fbrfdrlfUBD628vLL7+MrVu34u677w7r5zQ3N6O2ttbvQkS9l2Q/JAuiBQxZdWmusLMj3hV7E80H7HUSUQQEI+Xl5XA6nejXr5/fdrlfWloa9DmbN2/GHXfcgTfeeEONFwnHAw88gNTUVO9FMi9E1LeyIjrZHk52RC8Fn5HAYISoTwxgNRj8PzTkQyJwm5DARbpm5s+fj4MOOijs/c+bN0+NMdEvRUVFnXmZRBQF7C47Suv3tMqK6GR7aX2xaheKVF4trmpSt4trmryVWIkoOoSXqvDIysqCyWRqlQUpKytrlS0RdXV1+Omnn9RA1xtuuEFtc7lcKniRLMmnn36KE088sdXzLBaLuhBR72c2mbFocyEqTUY8kJGONfFxuLqqGic2uoMLkeHco9oFs2RdCeZ/uB4lte76JPd9tAEvLN+Ou2eNwYxxed32PoiomzIjZrNZTeVdunSp33a5f+SRR7Zqn5KSgrVr12LNmjXeyzXXXINRo0ap24cffvh+vHQi6i1yz3oGYxwamozuj6SJzTaMsdndF4emHg8ViFz7+iqU1LgDEV1pjVVtl8eJqJdlRsStt96Kiy++GIceeqiaGfPcc8+pab0SZOhdLHv27MGrr74Ko9GIcePG+T0/JydH1ScJ3E5EfdiE2UDWQShbcqG6m+Nwtjw253Mgv3VtIumKkYxIsA4Z2SYdx/L4KWNyYTIGH49CRFEajJx//vmoqKjAvffei5KSEhVULF68GIMHD1aPy7b2ao4QEQWyueyoMrkXuMtxSjAiAUTosR8/bq9slRHxJc+Ux6XdtOGZ/IUTRTCDFm55wx4kU3tlVo0MZpWuHyLqfYpLVuHUTy9FrKZh5Y4iGHInAvUlwFXLgNT+rdq/v2YPbl60pt39Pva7SThrUuvnE1HknL87nBkhIjoQykxGbxeNwZwMXP0V4LQBMcEHs+ckx4W133DbEVHP4UJ5RBQRyhrL1HW2dNGk5EkNgZCBiJg6NAN5qXEhqpO4O3nkcWlHRJGNwQgRRVQwkuNwACn57baXQakyfTcYPUCRxzl4lSjyMRghoohQ1uQJRiQzktx+MCKkjsjCiyYjJmC2TG5qnNrOOiNE0YFjRogoIuxr3OffTROmU8fmIi7GiHqbE/NmHowJA9JU1wwzIkTRg8EIEUVYN40EI+FlRkRNk10FIuKSaUMQb3ZPDyai6MFuGiKKrGCkA900oqjSXTY+O9nCQIQoSjEYIaKozozsqmxU1wPT4w/YayOiA4vBCBH1uAZ7AxodjT5jRjqQGanyBCMZCQfs9RHRgcVghIh63N7Gveo60eVCoiEGSMgK+7lF3swIgxGiaMVghIgiZiaN6qJJzgM8q/eGo6jKPWZkYAa7aYiiFYMRIoqswasd6KIRu/XMCLtpiKIWgxEiirDBq+HXGHG5NOzWMyPspiGKWgxGiKjH7WvSC55JKfjwV9jdW2eFzelSBc5kHRoiik4MRogowmqM5HW4xkh+WhxiPKv+ElH04f9eIoraGiOcSUPUOzAYIaKICUY6W2NkEAevEkU1BiNE1KNcmss7ZqSfPrW3g900nElDFN0YjBBRj6purobD5VC3szo8ZsSdGRnAUvBEUY3BCBFFRBdNhtOJ2MRsIMYc9nNZCp6od2AwQkRROXi12eFEaa1V3WaNEaLoxmCEiCJn8Gpy+MFIcbUVmgbEx5qQlRR+NoWIIg+DESKKnHVpOjOtNyMeBoPhgL0+IjrwGIwQUUSs2Otelyavw9N62UVDFP0YjBBRj9Kn9eZ0sBT8Li6QR9RrMBghosjppunAtN7dnhojnNZLFP0YjBBRRHTTuKuvhp8ZYfVVot6DwQgR9Ri7y45Ka6XPANaOFzxj9VWi6BfT0y+AiPqu8sZydR2jaUiPTQIsye0+x+nS8NXGMlQ12tX9/LT4A/46iejAYmaEiHpMWVNLjRFjGONFlqwrwdEPfYEr/vmTd9uMBV+r7UQUvRiMEFGPD17NDqOLRgKOa19fhZIad9VVXWmNVW1nQEIUvRiMEFGPD17t187gVemamf/hemhBHtO3yePSjoiiD4MRIoqMzEgb3TQ/bq9slRHxJSGIPC7tiCj6MBghoshYl6aNUvBldaEDkc60I6LIwmCEiHp8AGs/qb5qdxcxCyYnOS6s/YXbjogiC4MRIur5zIh00xSvCtlu6tAM5KXGIdRyeLJdHpd2RBR9GIwQUfer3gUUr8a++tKWRfK2LQOK16jt6nEfJqMBd88aE3RXeoAij0s7Ioo+LHpGRN1vwXg0GgyoHzKwpfqqvRJ47riWNvfU+D1lxrg8LLxoMm79v5/RaHN6t+emxqlARB4noujEzAgRdb9zn0dZrFndjHe5kKjJfBjPtFxjjHo8GAk4RvVLUrcvPmIw3rrqCHxz+4kMRIiiHDMjRNT9JszGsvotwOa3kOJ0+Y8FmfM5kD8p6NPsThfWl9Sp25cfNQTDst2BCRFFN2ZGiKjbaZqGRbs/U7frTEZPTqT9j6ONpXVodriQEheDIZmJB/x1ElH3YDBCRN2uoLgAe5rcBc8ajUYUTDoXyJ8IJOUAidkhn/fz7mp1PWFAGowcrErUa7Cbhoi6PSvyxOonYIABGjQY5L6xDkfO+QIGlx2IsYR87i9F7kGtEwemduMrJqIDjZkRIur2rEhhRaEKRIRmMKCwaiMKSr5rMxAJzIwQUe/BYISIuj0rYjT4f/TIfdkuj4fSaHNg01734NVJAxmMEPUmDEaIqNuzIi7N5bdd7st2eTyUdXtqIYvy9kuxoF8Ky74T9SYMRoio28eKBCPb28qO/OLpopnILhqiXofBCBF1C7vLjtKGUu9YkUCyXR6XdsGsKfIEI+yiIep1OJuGiLqF2WTGojMWodJaiT9+fhN2Ne3FH+psmPq7dwCDO1uSEZeh2gXzy27PTBpmRoh6HQYjRNRtchNzkWJOwe4m92q9p6WMQk7W2HafV9lgw67KRnV7/ABO6yXqbdhNQ0TdakPlBrigIcfhQE7+5HbbO10a3l7hXsU3LyUOSRZ+hyLqbRiMEFG3Wle+Tl2Pb7YBeRPbbLtkXQmOfugLPLRko7pfUmtV92U7EfUeDEaIqFut3fezuh7XTjAiAce1r69CSY3Vb3tpjVVtZ0BC1MeDkaeffhpDhw5FXFwcpkyZguXLl4ds+8477+CUU05BdnY2UlJSMG3aNHzyySf785qJKIqtK1ujrse5YoD0oSG7ZuZ/uD7ovBt9mzwu7YioDwYjb7/9NubOnYu77roLq1evxjHHHIOZM2di1y53n26gr7/+WgUjixcvxsqVK3HCCSdg1qxZ6rlE1LfITBp9gbyxmaO9s2gC/bi9slVGxJeEIPK4tCOiPhiMPPLII7jyyisxZ84cjB49GgsWLMDAgQOxcOHCoO3l8T/96U847LDDMHLkSPz1r39V1x9++GFXvH4iisLxIkNtdiTnHRKyXVld6ECkM+2IqBcFIzabTWU3pk+f7rdd7hcUhC7j7MvlcqGurg4ZGRkh2zQ3N6O2ttbvQkS9afBqc5vjRXKSwyv3Hm47IupFwUh5eTmcTif69evnt13ul5aWhrWPf/zjH2hoaMDs2bNDtnnggQeQmprqvUjmhYii39rd34Q1eHXq0AzkpYYONKRzRx6XdkTURwewGgL6eWUticBtwbz11lu455571LiTnJyckO3mzZuHmpoa76WoqKgzL5OIIoh8Tqyr2KBuj7c7gcwRIduajAbcPWtM0Mf0Txp5XNoRUfTrUPWgrKwsmEymVlmQsrKyVtmSQBKAyFiTf/3rXzj55JPbbGuxWNSFiHqB6l1AYwV2N+5DNRyI1TQcZLMDpWvdQ1ETMoG0Qa2eNmNcHsblp2BdsX83bW5qnApE5HEi6oPBiNlsVlN5ly5dinPOOce7Xe6fddZZbWZErrjiCnV9+umn798rJqLosmC8ulqXmADkZGGUzQazywE8d1xLm3vc6874qrXasWlvvbr9999OQKzJqMaISNcMMyJEvUuH6yrfeuutuPjii3HooYeqmiHPPfecmtZ7zTXXeLtY9uzZg1dffVXdlwDkkksuwWOPPYYjjjjCm1WJj49X40GIqJc793ngvWux1mJuGS+iM8YAZwefifflr2WwOV0Ynp2I30zhuDGi3qzDwcj555+PiooK3HvvvSgpKcG4ceNUDZHBgwerx2Wbb82RZ599Fg6HA9dff7266C699FK88sorXfU+iChSTZgNZB2Ego/PV3fjfQuVzfkcyJ8U9GmfFLq/uMwYl9s9r5OIekynVpy67rrr1CWYwABj2bJlnXtlRNRr2J0ObIuNVbeXJ8ThliojDHCFbG+1O/Hlr+7iaDPGcmwIUW/H5S+J6IB7p3wlNM+Muy1mMwr6j8ZRNfuAxGy/dlLeXaqqfr5hL5rsTuSnxmFc/xQeIaJejsEIER3wKb0vbH7be98II57IG4QjL/sChtiWWiKy8J2sN+NbBr66ya66azhzhqh346q9RHRAFRQXoNSzHo1wwYXCikIUlK1sd4XeRpuTK/QS9QEMRojogGZFHv7p4dYfPAYjnlj9hHq8rRV6dVyhl6h3YzBCRAc0K7K1emur7S7Nkx0pLuAKvUTEYISIDgzJejy++vGQjxtgUNmRvbVNYe2PK/QS9V7MjBDRAWF32VFUG3pdKQ0aShtKkZFkCmt/XKGXqPfibBoiOiDMJjPGJ/RHQc2vmFVXj4tmPgOk+VdSzYjLQHZ8P7UCb2mNNei4EYNnPRqu0EvUezEYIaIDoqKpAj/W/KpuX9FsxIjhp4ZsKwvfXfP6qlbbuUIvUd/AYISIDsgqvYt3fgKHWoumGSOsjUDxmpCr9J5wcA5S42NQ0yTPaMEVeon6BgYjRNS1FoxX3S3v5ecCFjPOqmsA7I1BV+nVK65+8PMeFYhkJ8XiH7MnoarRzhV6ifoQBiNE1LXOfR4bFt+ETRYzzC4NMxsagq7SG6ziqtWhqUJnZ03qz6NC1IdwNg0RdanvsgbiqoHuVbxPbGxEauAqvRNmh6y4Wm91sOIqUR/EYISIurS2yIJVC1DraFT3z6xvaPVR01bFVX0bK64S9S0MRoioy0hF1fUV61s2JOYAZzwK5E8EknLUKr0yRiQwIxIYkMjj0o6I+gaOGSGiLsuKSEVVnUEDnhoyFkdPuRyGKZcDThsQY0HZ9j1h7Y8VV4n6DmZGiKhLFKx9Va03o9MM8K4/A4NBBSIdqaTKiqtEfQeDESLqmrEivzwrN9pcnfe7rRUorWlCsiV0UlYKnUlFVlZcJeo72E1DRPtd4KygfC1+dda5MyBBVud9vOAjLPo6oc2xIoIVV4n6JgYjRLTfBc4e6J8HmGNDrs777Nqn0VhzvU+4ERwrrhL1TQxGiKjzzn0eNe9fh6LYmDZX5zXEVgMGJ6C1bpeRGIs/nzEWuSnurhmTse2AhYh6HwYjRNR5E2Zj4d6v4Sr+Ctl2B/6+rxxx+riRc1/AL839cOe7a6E5koIGIqKywa4CkWnDM3kkiPooBiNE1Cnf/fIa7ln1DxTDqXpf7i+vxORmm2dcvAvO5MFYvDMBLmv7pd05jZeob2MwQkSdnD3zDIql6wXAWXX1ODJ5MHDyNcCqV9FcUYTf/nMTfqlNCmt/nMZL1LcxGCGiTs2eWe+s9W4+tv+xwNSbgcQsLLHMwM1v/IhmBB/Q6svgGbTKabxEfRuDESLq3OwZGbRqMMCgaXhp33c45fn3VHAx3/JO2IGIuHvWGA5aJerjWPSMiMJ37vN4PD0dO2Uar6emiGYwoNBiwTfxifh4xPx2a4noJCOy8KLJmDEuj0eAqI9jZoSI2rdnFbD0L3hr+GF4ITW51cNGTcPc1INRsW5EWL/NG04YjltOGcWMCBEpDEaIqF3amrewsPoXLNy2M2jdMpfBAFv8PpgSN8PZcFC7+ztqRDYDESLyYjcNEQX13ab3cda/T8W3v7yG/7frQyxMT3U/ELD+jJcGWLI/dd8IgevOEFEwzIwQUSva7pV4bNkfsc1iwZ9+egC1CbEqCEnUNDQYQ3yHMaDNSqscsEpEoTAYIaJWCn56Ug1KFbUmkxoT8tC+CkyyNqPSZIRDM2GB4zws0yb5Pa+tSqtcd4aIQmEwQkR+NUSkF+bhih8Bk6Q6DCojMtDuwKkNje66IE4nTm++B4Xa0LB+czecMAJHjcjiujNEFBKDEaK+zjNTBjuWozjGhD9nZWJrfFzL4waDmspbEB+HaY3NMBpCjwkJVtDsllMO4mBVImoTgxGivu7nRajd9Q1eGDQarxvqYQ+yaq5008xPG4C/NTiRj0pUaClt7pLjQ4ioIxiMEPXBBe4eXPUP3DF8Ng4ZcgoWbXkXLwzIR42pwSeMaD11tyTOhvNiL4epYRhs7VRY5fgQIuoIgyYrXkW42tpapKamoqamBikpbX8jI6LQHE4Xzn/jaGzS6jDAbleTcPfEugOLYTYbHDCgKDZGVVUNpGkGtQJv447rQwYtafGxeOr3k3HEsEx2zRARwj1/MzNC1Es5XRoKf1qGvB//ivLxc1CBDPx9xRfYml2nHt/tCUKyHE5cX12N0+oacNrA/kEDEWEwaCGn7urPePC88WqwKhFRRzAYIeqFg1ELht2EP3xrwv80PIuLY37E5u9/waeJCdiRldgSOmgasp1OfLi7RNUPEYuKS3GZ8zpswkBPOw2xcMLu+agINXWX3TJEtD8YjBBFedZj6djzsKjmv7j9sDsw7Ou3kLdjOfZusSHVMgFb+63CCYn9UW0ytd6JwYB9MTFYbbHgaKsVLs2gpu2amnPh0gZ4mzUHPs1TY/WWk0diSFYicpLjOG2XiPYLgxGiKAs+Put/Ax7bkISrG5/FZTE/4pVNlWpw6dyP78JfK7bhxcx0LE3cjUpTCfYgXj0/3SldK0C1yejXDSOzZP6eno2Pd52C35mWIQ8V7c6UYRaEiLoagxGiCA48Sqbeid1xB+O+j9erLpfLY35EYunTyHTOwizzd/g4MUEFIsJq2Ydb81tW1E11OnFyQxNObWiA3WDA9bk5QWfJbI0D1sYNwlsN98EMR8iZMlceNQQnj8llFoSIuhyDEaIeDDTGHnq82i7bilfejQeTYtBQ9htcV7VMBR5ffvAIXnVORwYMOD7uexRY4lBqWY/hli34rSUJ5TH+/4Ul03FmfQNmNDRiapNVhRXSpXJBfj+ZOhd8cKpngbvGhpFBA5G81DjcPWsMZozLO3C/GCLq0xiMEO1HYPHj9kqU1VnVuIkpg9OxcmeV9/7UoRmqXbDuFQk0PvpgIS5ZYldtbrI9j08HlaLCZYE58d84pW4z1ljMsKWsxlhLIQotZpwZ61k1tw2S6ZBA5Kgmq/u+ZoDDoKHUFHy6buACdwYthuNBiKjbMRihPi8wqNCDiMBA49dVX7XqOsmqXY+zk1/F/Kw41O49D6MbTZgX8yb+6rgQRfGj1H5utr3QqntFnGX6FquahwAx1WjI+AmFliS13Ra/D6cPSYXLkNbq2Ay02zGu2YaxzTb8JzkJO2NjVADimxl5Ij0VRzRa1dIyW7R8vOKYgb/s/gLmmGrcYr8eFWgd1OizZDgehIh6AoMR6pPZCr3NjvJGvPXjLmTXrW8VRAxs2ujdVohh+LPpn96MxnzHparNVTFf49OsGtgtVpizP8E5u5NwpGk9LtY+xatNgAvAkfHfY02sGYbEX3BG7HosiInB7pgc7I6NQZlpscpYLIM7ENFJgJHhcGJCczPG2mwY7wlA0lyyR+Db+DhsN7fuUpHnyWq7b8YNwanWSlxiux2lyMKbTScGHQ+Sm2LBBVMHcVYMEfUoBiMUcQFCR9qEylb4BhFjsc17f1e6Aa70d9Gv9FDca/saLzguRKk2DFfHLFdBxLnacsxvHKZ+1s2ebVdpH+FZyWiYWjIaa1zD1e3cpBUotLgHjZrid2N7TiMeQDp2x25ATMwm7ImJwW+MkokI3cVidrlgMxpbbf9/5RU42tPdone5SJeKU4PKfoQaAyKVUh9Iz8N9O+6BDWbPVgMyU5Px59NHIz3R4vd7NQVZi4aIqDsxGOkDwu2G6KoAwbdNVYOt3QBBz0RkGT6HIfe/0EpnYlvjsUHbhJOt+B+fwGKtYxjO9dw/R/saf0+ph8lcBuR8gWklm1oFGmeYvsWPhnwYTQ0YYvken5ji0WBah9NNG/CM0YTKmCxUGY2oNL2LSpMRXxj9MxrvpSa0+v1L10muw4kBDoe69Ld7rj23b8jNxgazuVV3y5PpqWqVXJNBQ6WWhCItB287TsC5pi9QbLK3Uym1BjaDUQ1O5SwYIop0fTYY6ckT9IFsE2rQZDjdEG1lENoKItraz1ptWJsBQksmQkPW0C9QatbQL+cLuLYfE6RNeNkKuS2TXY+P/RY/m9Iw3lyAFSYLai0/wRTvDh72xllxY04WLNpm9DM9hjmmBFSZklFtNMJl+Fy1+bqNbEZbptc34AirFf3tTgx0OJDrcPh1jkiGw2jQ1PV3CRbVrRK6u2UwTrVW4bzmu1GEfirD8abzRJi3V8AR01KOLDnO/V+5zurwjgHJS0niLBgiigp9MhhZsq4E8z9cjxjTm6jL/hFxNaPQkFgFY9U5GFSltfsNvq2T7/60GZbwtTczUK6d1Kn9BA6atJa+gVLHpW12Q4STQQgdRITejxo34ZwODQZv0HC26Rusdg3HOaZv/IKIioQy/Bxn9QYKx6e/iDHOFfjOaEGm4XucixpoxmbUx+zEI8Y01Bs3YqxxM+42WNBgzEGD0Yh647vq+mRjOuwGd1AGfIfPkR707+CrxNZZDN8aHelOFzJcTmTIted+usuJTKcLaU4nHspMx7bY2FYZjT2xMfjNvgbvei3u7hUJPgCjz6DS2aYv8Ghacwe7W4QBGQn9W431CBaYsguGiKJBTF8MRK59fRU0uJA+agWcRqAhbaPqi7enfIxzat0DEDtz8t2/NkPViV/PDFRvP0addDqyH33QpDr5m9s++eu3g7XJT1oBU3yyNzD4JCEeRzZ/i3NjEjE+tgBrjGYMsnyH32hWNR00MWYd/s+YBKfhB5yplaDatBt/N6bBaliLQwzr0GwwYL7BjGZjtrptNfwHlxmSYTWkqvtNxg/QoE7GLSfklblbMMcviCj2XPt3i7QnweVCksslsQD2xrb+c59dW4dDrc2qQqkEG5kuJ1KdrqBlvwIzGlvM5pAZjW/i4nGMtcmve+V805fIQ3nLoFLXsUiKeRAGQ0Ob3S0ZqQn4y+njwxrrMW14Zod+P0REkcCgaZ4VsnrBEsThdM385sFFsNeVozpzNWqyV7Rq83BJHWZYq1ClJeJu+6W4N/afSDM0qJPKfPslqs3dsa8iw1DfoTbzY/+JFEMDypCEB+2zYTQ4cVPsf5BkbESVFo+/mqdhZd4m7+s4rGwQspvjcarpB8QabaiBGV+5xuHwmEKYDA7UIRZrtKHQDE6MNBUBBhfsElAZDLDDAJt+W7/Afd+h7kO1cT/m3m7zeV690ahqTcjaJT0p0yFZCCeSNBcSXJoKKhLlWpNrCTL0bS4kaZon8GjZlqBpMPkU/Qo2LmO0zYa3iveqLIQeaMi1f/DhzmhscvXHK053RuP+/s3YZDGFLCI2tBl4vrgGv/HpXpEHAme0pCbXA6YGb/eKet+JsZgxNg/5aXEYnNYPpx58MDMcRBSVwj1/dyoYefrpp/Hwww+jpKQEY8eOxYIFC3DMMfJNPrivvvoKt956KwoLC5Gfn48//elPuOaaa7r8zbTnu60VmPbaMDXd8vDBA2CVE4lc5FfguZapk0c0WVUbp8EAp+ck7pRZDJCTeRvXcrJHiOsePrHvrxiXhnhNQ5zmQpymwaJua7C4NL/7gdtatrtatVWPuzTcmZ2JrebYkIFCW785VxhBxPK4eFyXlx1yH8+UliG7IVMFGpK96G/Yp7bv0bLxttOT0TCU48zm+1VGAwY7kkY8CENM8IyGeg2OJMTvvgMOQzyqG92FzUJNpRXsXiGi3ijc83eHu2nefvttzJ07VwUkRx11FJ599lnMnDkT69evx6BBg1q13759O0477TRcddVVeP311/Htt9/iuuuuQ3Z2Ns477zx0J0lx32y7DgdnvQGr71RK/SRoMKjVTZckyTLr3UNOukYNcARJufezO5CiuRCraYjV4L6WJd01wKy2eS5yIL33fdv53zfrt32e13IfMGka5uVkBh0HMcrefmDQlsAAQb8vgcJmS+jujoL4ODWjJFSgoY+/UEEEfIIIxwm4yLQUB6EI/0jPdqdHgrx42e2jaRkoqvoT9iJbDQ7Njnc33NfkfpJsizM44C6uLgmOWCSW34Y5x2UjOd6MqoZm1YUyOi8FG0pq1X09oxFuoMHuFSLqyzqcGTn88MMxefJkLFy40Ltt9OjROPvss/HAAw+0an/77bfjgw8+wIYNG7zbJCvy888/47vv3OMaujMzcsHzBUgZ9RdoBnvQbggZTNjP6cTl1bUqUovxpPqDXmsSzbVcx8gJPdR9TVP7M/lcy34MYXQhGLrw5L+/GQQ9MGhrP4E/yztuwukOEEYZirBRG4jXnCejcNBn2GEJPoBTdp/SnIIn9zRggKG87WyFpwtE2NRv2N0tkmysg2v44zDG1Id8XwZnMi4b+DyGZ6d1aCYTB4cSEfVAZsRms2HlypW44447/LZPnz4dBQUFQZ8jAYc87uvUU0/Fiy++CLvdjtjY4CuEHghyAsno/w3sRjlpBT/Fy0mxNCYGA+0OHGO1hnUS3582Be1M7dQHQob3swJO/p7swCi4T/6vO07ZrwzCE2mpSGuw4I029hP4s2RJ+sBpqcmGJtRp8Z7Br9/DaAgeKMgklKoYF86234PsuNi2sxUhinpJELF00yTsqt6rto/Ld/9nWFfszmDItmOHDUF+cl672QpmL4iIDowOBSPl5eVwOp3o109OLC3kfmlpadDnyPZg7R0Oh9pfXl7rk0Bzc7O6+EZWXcEgZ9rUz901uttKN2jAP9JzMK24CJvaOYmHc6IP1Wa3loXb0xJh0GrVibfV69U0PJyejcW7TsbFps/a/VmhTv6JaEID4r33g3VDqAxCzONtBgaFMRk43X6H6qYIuR890PDcH5QSg3nnTWgVIOhZhhjzSxic7WoVIOhdHrbmeIw8fuB+ZStOHzNa8nd+28aFTgAREVE0TO01BKTUpacncFt77YNt10l3z/z589HVGh2NcGi2tgMR9cKAzTFJOMS+EHWu5LZP4gEn3w61sdqRGPOQyniEytJsjUnCL67j8ZbzpLB+lmyzJ8QC3kGTBiSnpON/2hk0GU4G4aihQ7B7n7nDhdnaHiPRP2SAMDGn9e+E2Qoioj4ejGRlZcFkMrXKgpSVlbXKfuhyc3ODto+JiUFmZvCaCPPmzVOzb3wzIwMHDsT+SjIn4fXTXseuul2oaqrGr2X70NjsQIIlBgPT3QWwGhuTEY98NQDx5EtHHfDKqZsqhsBiafIb/OgbDOzaZ4R9SkqHK7B2ZtBkOBmEgUG6/BggEBFRtw9gnTJlippNoxszZgzOOuuskANYP/zwQzXbRnfttddizZo13T6AlYiIiLpPuOfv1kuFtkMyFi+88AJeeuklNUPmlltuwa5du7x1QySrcckl7sJfQrbv3LlTPU/ay/Nk8Optt93W2fdGREREfXnMyPnnn4+Kigrce++9qujZuHHjsHjxYgwePFg9LtskONENHTpUPS5By1NPPaWKnj3++OPdXmOEiIiIIlOfKgdPREREvaCbhoiIiKgrMRghIiKiHsVghIiIiHoUgxEiIiLqUQxGiIiIqEcxGCEiIqLoW5umu+mzj7tqwTwiIiI68PTzdntVRKIiGKmrq1PXXbE+DREREXX/eVzqjUR10TOXy4Xi4mIkJye3uTpwR+kL8BUVFfXaYmq9/T3y/UU/HsPox2MY3WoP4HlCQgwJRKT6utFojO7MiLyBAQMGHLD9yy+/N56o+9J75PuLfjyG0Y/HMLqlHKDzRFsZER0HsBIREVGPYjBCREREPapPByMWiwV33323uu6tevt75PuLfjyG0Y/HMLpZIuA8ERUDWImIiKj36tOZESIiIup5DEaIiIioRzEYISIioh7FYISIiIh6VJ8ORp5++mkMHToUcXFxmDJlCpYvX45o9MADD+Cwww5TFWpzcnJw9tlnY+PGjX5tLrvsMlW91vdyxBFHIBrcc889rV57bm6u93EZgy1tpMJffHw8jj/+eBQWFiKaDBkypNV7lMv1118flcfv66+/xqxZs9Qxkdf63nvv+T0ezjFrbm7GjTfeiKysLCQmJuLMM8/E7t27Eenvz2634/bbb8f48ePV65Y2l1xyiaoi7Uvec+Ax/d3vfodoOYbh/E1G6zEUwf4/yuXhhx+OimP4QBjnhUj6f9hng5G3334bc+fOxV133YXVq1fjmGOOwcyZM7Fr1y5Em6+++kqdtL7//nssXboUDocD06dPR0NDg1+7GTNmoKSkxHtZvHgxosXYsWP9XvvatWu9j/3tb3/DI488gieffBIrVqxQgcopp5ziXdMoGsjr9n1/chzFb3/726g8fvK3N3HiRHVMggnnmMn/z3fffReLFi3CN998g/r6epxxxhlwOp2I5PfX2NiIVatW4c9//rO6fuedd7Bp0yb1IR7oqquu8jumzz77LKLlGIbzNxmtx1D4vi+5vPTSSyrYOO+886LiGH4Vxnkhov4fan3U1KlTtWuuucZv28EHH6zdcccdWrQrKyuT6draV1995d126aWXameddZYWje6++25t4sSJQR9zuVxabm6u9uCDD3q3Wa1WLTU1VXvmmWe0aHXzzTdrw4cPV+8v2o+f/C2+++67HTpm1dXVWmxsrLZo0SJvmz179mhGo1FbsmSJFsnvL5gff/xRtdu5c6d323HHHaeOczQI9h7b+5vsbcdQ3uuJJ57oty2ajmFZwHkh0v4f9snMiM1mw8qVK1WU6EvuFxQUINrV1NSo64yMDL/ty5YtU+m6gw46SEXzZWVliBabN29WqUTpVpM06LZt29T27du3o7S01O9YSuGe4447LmqPpfx9vv7667jiiiv8FoaM5uPnK5xjJv8/pbvDt40c/3HjxkXlcZX/k3Is09LS/La/8cYbKv0tmb/bbrstqrJ57f1N9qZjuHfvXnz88ce48sorWz0WLcewJuC8EGn/D6NiobyuVl5erlJM/fr189su9+XgRDMJ8m+99VYcffTR6g9GJ11QkvIfPHiw+iOUFPKJJ56o/tgivTrr4YcfjldffVV94MmHwv33348jjzxS9W3qxyvYsdy5cyeikfRdV1dXqz753nD8AoVzzKSN2WxGenp61P8ftVqtuOOOO3DhhRf6LUL2+9//XgXXkhpft24d5s2bh59//tnbRRfp2vub7E3H8J///Kcae3Huuef6bY+WY6gFOS9E2v/DPhmM6Hy/deoHLHBbtLnhhhvwyy+/qL49X+eff773tvwxHnrooepDRKL9wP9gkfihp5NBgdOmTcPw4cPVB4Q+YK43HcsXX3xRvWf5BtIbjl8onTlm0XZc5VulZPJcLpcaMO9LMgm+x3TkyJHquMo4k8mTJyPSdfZvMtqOoZDxIhJ4yGSHaDyGN4Q4L0TS/8M+2U0jKTWTydQqspMUY2CUGE1kxPMHH3yAL7/8EgMGDGizbV5envrgkO6PaCMjuiUokdeuz6rpLcdSvpF89tlnmDNnTq89fuEcM2kj3VVVVVUh20RDIDJ79myVNZBvyu0tzS4nr9jY2Kg8psH+JnvDMRQyy1JmobT3fzJSj+GNIc4Lkfb/sE8GI5J2kqm8gak0uS/p/2gjUapEvjJq/4svvlBpw/ZUVFSgqKhIfYBEG5lqtmHDBvXa9RSp77GU/zwykjwaj+XLL7+s+uBPP/30Xnv8wjlm8v9TPtR928hMBUmFR8Nx1QMROSlJcJmZmdnuc6TbUZ4Xjcc02N9ktB9D30ylvBeZeRNNx1Br57wQcf8PtT5KRgfLKOEXX3xRW79+vTZ37lwtMTFR27FjhxZtrr32WjUCetmyZVpJSYn30tjYqB6vq6vT/vCHP2gFBQXa9u3btS+//FKbNm2a1r9/f622tlaLdPLa5b1t27ZN+/7777UzzjhDS05O9h4rGQ0u7/+dd97R1q5dq11wwQVaXl5eVLw3X06nUxs0aJB2++23+22PxuMnr3n16tXqIh8zjzzyiLqtzyYJ55jJbLcBAwZon332mbZq1So1k0FmVTkcDi2S35/dbtfOPPNM9drXrFnj93+yublZPX/Lli3a/PnztRUrVqhj+vHHH6vZfIccckhEvL/23mO4f5PRegx1NTU1WkJCgrZw4cJWz4/0Y3htO+eFSPt/2GeDEfHUU09pgwcP1sxmszZ58mS/qbDRRP4jBbu8/PLL6nH545s+fbqWnZ2tAjA54cm0vF27dmnR4Pzzz1f/QeS15+fna+eee65WWFjofVymqMn0X5mmZrFYtGOPPVb9x4o2n3zyiTpuGzdu9NsejcdPTk7B/ibldYd7zJqamrQbbrhBy8jI0OLj41UQGinvua33JyemUP8n5XlC3oe8Z3lv8vkj07hvuukmraKiQosUbb3HcP8mo/UY6p599ln1umWKa6BIP4Zo57wQaf8PDZ4XTURERNQj+uSYESIiIoocDEaIiIioRzEYISIioh7FYISIiIh6FIMRIiIi6lEMRoiIiKhHMRghIiKiHsVghIiIiHoUgxEiIiLqUQxGiIiIqEcxGCEiIqIexWCEiIiI0JP+P/0HaaorT133AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "delta = 0.01\n", "\n", "for N in [100]:\n", " TEI = TestEulerIntegrator()\n", " TEI.delta = delta\n", " TEI.do_integration(N, 0.0, 2/delta)\n", " plt.plot(TEI.x, TEI.y, 'o-', label = 'Euler: ' + str(N))\n", "\n", "for N in [100]:\n", " THI = TestHeunIntegrator()\n", " THI.delta = delta\n", " THI.do_integration(N, 0.0, 2/delta)\n", " plt.plot(THI.x, THI.y, '*-', label = 'Heun: ' + str(N))\n", "\n", "TRKI = TestRK45Integrator()\n", "TRKI.delta = delta\n", "TRKI.do_integration(0.0, 2/delta, eps_allowed = 1e-5)\n", "plt.plot(TRKI.x, TRKI.y, '^-', label = 'RK45')\n", "lgnd = plt.legend(loc='best')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.12.13" } }, "nbformat": 4, "nbformat_minor": 1 }