{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Numerical Methods \n", "\n", "\n", "# Lecture 2: Numerical Differentiation\n", "\n", "\n", "## Exercise solutions" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# some imports we will make at the start of every notebook\n", "# later notebooks may add to this with specific SciPy modules\n", "\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2.1: Compute first derivative using forward differencing\n", "\n", "Use the forward difference scheme to compute an approximation to $f'(2.36)$ from the following data:\n", "\n", "$f(2.36) = 0.85866$\n", "\n", "$f(2.37) = 0.86289$\n", "\n", "You should get an answer of 0.423." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.42299999999999693\n" ] } ], "source": [ "dx = 2.37-2.36\n", "df = (0.86289-0.85866)/dx\n", "print(df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2.2: Compute first derivative using central differencing\n", "\n", "Use the data below to compute $f'(0.2)$ using central differencing:\n", "\n", "$$f(0.1) = 0.078348$$\n", "$$f(0.2) = 0.138910$$\n", "$$f(0.3) = 0.192916$$\n", "\n", "You should get 0.57284" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.57284\n" ] } ], "source": [ "dx=0.1\n", "df = (0.192916-0.078348)/(2*dx)\n", "print(df)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example: Write a function to perform numerical differentiation\n", "\n", "As covered above, the expression\n", "\n", "$$\\frac{f(x+\\Delta x) - f(x-\\Delta x)}{2\\Delta x},$$\n", "\n", "can be used to find an approximate derivative of the function $f(x)$ provided that $\\Delta x$ is appropriately small. \n", "\n", "Let's write a function `diff(f, x, dx = 1.0e-6)` that returns the approximation of the derivative of a mathematical function represented by a Python function `f(x)`.\n", "\n", "Let's apply the above formula to differentiate $\\,f(x) = e^x\\,$ at $\\,x = 0$, $\\,f(x) = e^{−2x}\\,$ at $\\,x = 0$, $\\,f(x) = \\cos(x)\\,$ at $\\,x = 2\\pi$, and $\\,f(x) = \\ln(x)\\,$ at $\\,x = 1\\,$, i.e. functions we know the exact derivative of.\n", "\n", "In each case, using $\\,\\Delta x = 0.01$, let's write out the error, i.e. the difference between the exact derivative and the result of the formula above." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The approximate derivative of exp(x) at x = 0 is: 1.000017. The error is 0.000017.\n", "The approximate derivative of exp(-2*x) at x = 0 is: -2.00013. The error is 0.00013.\n", "The approximate derivative of cos(x) at x = 2*pi is: 0.00000. The error is 0.00000.\n", "The approximate derivative of ln(x) at x = 0 is: 1.00003. The error is 0.00003.\n" ] } ], "source": [ "def diff(f, x, dx=1e-6):\n", " numerator = f(x + dx) - f(x - dx)\n", " derivative = numerator / ( 2.0 * dx )\n", " return derivative\n", "\n", "dx = 0.01\n", "x = 0\n", "f = np.exp\n", "derivative = diff(f, x, dx)\n", "print(\"The approximate derivative of exp(x) at x = 0 is: %f. The error is %f.\"\n", " % (derivative, abs(derivative - 1)))\n", "x = 0\n", "\n", "def g(x):\n", " return np.exp(-2*x)\n", "\n", "derivative = diff(g, x, dx)\n", "print('The approximate derivative of exp(-2*x) at x = 0 is: {0:.5f}. The error is {1:.5f}.'\n", " .format(derivative, abs(derivative - (-2.0))))\n", "\n", "x = 2*np.pi\n", "f = np.cos\n", "derivative = diff(f, x, dx)\n", "print('The approximate derivative of cos(x) at x = 2*pi is: {0:.5f}. The error is {1:.5f}.'\n", " .format(derivative, abs(derivative - 0)))\n", "\n", "x = 1\n", "f = np.log\n", "derivative = diff(f, x, dx)\n", "print('The approximate derivative of ln(x) at x = 0 is: {0:.5f}. The error is {1:.5f}.'\n", " .format(derivative, abs(derivative - 1)))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2.3: Compute the derivative of $\\sin(x)$\n", "\n", "Compute \n", "\n", "$$\\frac{d(\\sin x)}{dx}\\qquad\\textrm{at}\\qquad x = 0.8$$\n", "\n", "using (a) forward differencing and (b) central differencing. \n", "\n", "Write some code that evaluates these derivatives for decreasing values of $h$ (start with $h=1.0$ and keep halving) and compare the values against the exact solution.\n", "\n", "Plot the convergence of your two methods." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact derivative at sin(0.8) = 0.6967067093471654\n", "Forward differencing Central differencing\n", " 0.256492 (error= 0.44) 0.586258 (error= 0.11)\n", " 0.492404 (error= 0.2) 0.668038 (error= 0.029)\n", " 0.600269 (error= 0.096) 0.689472 (error= 0.0072)\n", " 0.650117 (error= 0.047) 0.694894 (error= 0.0018)\n", " 0.673843 (error= 0.023) 0.696253 (error= 0.00045)\n", " 0.685386 (error= 0.011) 0.696593 (error= 0.00011)\n", " 0.691074 (error= 0.0056) 0.696678 (error= 2.8e-05)\n", " 0.693897 (error= 0.0028) 0.6967 (error= 7.1e-06)\n", " 0.695304 (error= 0.0014) 0.696705 (error= 1.8e-06)\n", " 0.696006 (error= 0.0007) 0.696706 (error= 4.4e-07)\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAFYCAYAAAAmzDckAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd1zV1f/A8dcBxL0HGlZqWu5JKo5CrdTU0rShtqwEXGhqpllqWj+zcmeiado3zS3mnok5KFc5Sc1BriTci31+fxxAQFDAC5/L5f18PO6DPud+xvtyP/H2nM8ZSmuNEEIIIe5wsjoAIYQQwt5IchRCCCGSkeQohBBCJCPJUQghhEhGkqMQQgiRjCRHIYQQIhlJjsLhKaU8lVILlVLnlFKRSqmLSqkNSqm3lFLOVscnHoxS6pRSanYGjntbKfVOJoQkHIAkR+HQlFL9gO1AMeBD4BngHeAoMBVoa110wmJvY+4FIe7iYnUAQmQWpdRTwDjgG621X7K3f1ZKjQPyZ31kDy6uxqu01tFWxyKEI5Kao3Bkg4FLwKCU3tRaH9da74/fVkrVV0ptVErdUErdVEptUkrVT3yMUmq2UuqMUqqOUmqrUuqWUuqYUso32Xm0Uqpd8msqpaYqpf5TSuVKVNZdKbVPKRWulApTSs1UShVLdpxWSn2ulBqslDoJRAI14t6rGxdLuFLqtFLqI6XUp0opnewcLkqpIUqpv5RSEXHNzGOVUnkS7VMu7lo+SqmRSqnzSqkrSqkVSqmyKXye7kqpvUqp20qpy0qpLUqpRonez6eUGqOUOhnXpH1SKTVUKXXPvz2J4uiplBqnlAqN+12vVEqVu9exib6DVL9LpVQg8DTQOO46Oq5MCENrLS95OdwLcAZuAT+lcf+awG1gD9AJ6AjsiiurlWi/2cA1IBjwAZ4FfgI00CzRfn8BC5NdwxW4CExOVPYFEAWMBZ4DugFngd8B50T76bjyrXGxtQLcgBLAZeAQ8ArQHtgChJj/vZNcfz5wExiGaV7uA1wBliTap1zctU7Ffa7WwFtAGLAl2fm+jtt3BtAOaAOMAl6Le98lLt6LQD+gBTAUCAfG3uf7iI/jNLAi7tzdgPOYJvFcifY9BcxOz3cJVAX2AvuAhnGvqlbft/Kyn5flAchLXpnxikscGhidxv0XxyWKIonKCmFqnksTlc1OIRHmjkse0xOVDY37Y1w4UVn7uGPrx22XA2KAYcliaRy3X/tEZRo4B+RNtu//YWqRZROV5QUuJE6OQNO4c7yZ7PiuceW1E8WkU0iEA+PKH4rbrhgX+7h7/E7fiDvmqWTlQ+NiLnWPY+PjOAw4pfC7eTdRWfLkmNbvMhDYZvW9Ki/7fEmzqhDGU8BKrfWV+AKt9TVgOab5LbFbWuvNifaLAI4BjyTaZw4mab6cqOwN4IjWemfc9rOYRxtz45o8XZRSLpha47W4mBJbq7W+naysIRCktT6TKJ7bwKpk+7XCJKQlya61PtHnTyz58QfifsZ/xmfiYp9O6lpharA7UrhmrrjY72ex1jo2fkNrvR04A3je45j0fJdCpEiSo3BUFzE1t0fTuH8xTJNdcv8CRZOVXU5hvwgg4dmd1joE+BV4HUApVQTTNPhjomNKxf38G9O0mvhVCCie7BopxVcGCE2h/EKy7VKYZt0bya4Tf2zya11Kth0R9zP+M8bvf4bUlcL8/pN/tvh/HCS/ZkqSf474Mvd7HJOe71KIFElvVeGQtNbRcR0snlVK5Y6r3d3LJaB0CuWluTtRpNWPwHdKqUeBlpjkNDfR+xfjfj5Hygn3YrLtlNaXO8+dJJuYWwrnCsc0r6bkXCrlqQmL++kOHElln4vAScyz0JScSsN1kn+O+LI/73FMZnyXIoeRmqNwZF9gaidfpfSmUqq8Uqpm3OYWoI1SqmCi9wtiOppsyeD1F2ESUldMk+qvWutTid7fAMQCj2itd6fwOpmGa/wGeCbuSaqUyouppSa2FlPrK5zKtdKbHDfGxe59j33WAg8DN1K5Ztg9jo3XKXHPVqVUY6AsEHSPY9L6XUZgns8KcRepOQqHpbX+VSnVHxinlKqC6UzzD6ZprQXwHtAF2I/pZdkW2KSUGoOppX0I5ANGZvD615RSy4FemObP7snePx53rW+UUk9g/nCHYxLKs8CMxM82UzEO6AGsU0p9ivmD3z/uZ0JNU2sdqJSaByyOG9+5E5PcygHPAx9qrY+m47MdV0qNB/rHJZ7lmA469YG/tNYLMLXkbpjf6VhMz1BX4DHgBUyHo1v3uVRBYJlSahpQEhiNeb77v3sck9bv8jDQUyn1KnAcuK61Tq0WLHIaq3sEyUtemf0CGmFqcecxz7wuYTqFvE7SnpANMDWiG5ghD5uI61maaJ/ZwJkUrhEIBKZQ3gbzxzlJz9Vk+7yBqQHejLt2MPANSXugauCzVI6vC2zDJNazwCfAROBysv2cgL6YJBUOXI377y/jY+NOL9H3kh3rFVfulazcF/OPi4i432sg4Jno/TzACMzQlvh9dsWVudzjO4uPoyfmHwD/YYbmrALKJ9v3FIl6q6bjuywNrAaux13rru9PXjn3pbRO6TGGECK7ips9Zy8QprVuYXU8GRE30P8k0F1rPcPaaEROJM2qQmRzSqlRmB6vIZhnrO9hBsI/b2VcQmRnkhyFyP40Ztabh+L+ez/med4aS6MSIhuTZlUhhBAiGRnKIYQQQiSTY5pVS5QoocuVK2d1GCIDbt68Sf782XJlKWFjci+IeLa6F/bs2ROmtS6ZvDzHJMdy5cqxe/duq8MQGRAYGIiXl5fVYQg7IPeCiGere0EpFZJSucM3qyql2imlpl+9etXqUIQQQmQTDp8ctdYrtNbehQsXtjoUIYQQ2YTDJ0epOQohhEgvh0+OUnMUQgiRXg6fHIUQQoj0cvjeqkqpdkC7ihUr3nO/a9euERoaSlRUVNYEJtKscOHCBAcHWxpDrly5KFWqFIUKFbI0DiFE1nD45Ki1XgGs8PDw6J7aPteuXePChQu4u7uTN29elFJZGKG4n+vXr1OwYMH775hJtNbcvn2bs2fPAkiCFCIHkGZVIDQ0FHd3d/LlyyeJUdxFKUW+fPlwd3cnNDTU6nCEEFnA4ZNjWnqrRkVFkTevLAgu7i1v3rzS7C5EDuHwyTGtvVWlxijuR+4RIezD/PkwfXp5goIy7xoO/8xRCCFE9nf6tEmKM2bA0aMAj7BsGWzaBJ6etr+eJEchhBB26b//YPFi+Okn2LbNlJUtC0qB1orISAgMzJzk6PDNqjJDTvbSu3fvdE8mvHv3bpRSnDp1CjATEiulCAsLS9jn559/plKlSri4uPD222+nWiaEsNa1a/C//0Hr1lCmDPTsCZcvw2efwd9/w8KFkCcPODnF4uoKmTUPvcMnR0eeIeftt99GKXXX688//7Q6NEs1atSI8+fPU7x48YSy9957j44dOxISEsLEiRNTLRNCZL3wcFiyBDp1Ajc3eOst+OsvGDQI9u+Hgwdh6FB47DFTS9y0Cd5551SmNamCNKtme8888ww//vhjkrISJUpk+HxRUVHkypXrQcO6p9jYWLTWODs7Z8r5XV1dKV26dML2lStXCAsLo2XLlri7u6daJoTIOtHRJsnNmwdLl8L16yYxdu8OnTtDw4am+TQlnp4QEfEPnp4VMi0+h685ZrWgIBg9mkztRZVY7ty5KV26dJKXi4v5N09ERAT9+vXDzc2NPHny0LBhQ7bFN9xzp/lx9erV1K9fH1dXV9atW0euXLn4/fffE/YrW7YsVapUSdjesGED+fPnTxjWMG7cOGrWrEn+/Plxd3fnvffe48qVKwn7z549mwIFCrB69WqqV6+Oq6srwcHBxMTEMHDgQIoWLUrRokXp168fMTEx9/3Ma9eupXLlyuTJk4emTZty1Dydv+tzhYWFERgYSNGiRQFo3rw5SqlUy4QQmSs21jw77NULHnoIWrWCZcvg5ZdhwwY4cwYmTTLJz+rO4VJzTEW/fpDe1smrV00TQGwsODlBzZqQntbc2rVhwoT0XfNeBg0axMKFC/n++++pUKEC48aNo1WrVhw7dowyZcok7Pfhhx8yduxYKlasSMGCBalbty6bN2+mQYMGHDt2jKtXr3Lx4kXOnz9PmTJlCAwMpFGjRgk1TCcnJyZMmECFChUICQmhT58+9OnTJ0mNNjw8nM8++4xp06ZRsmRJypQpw9ixY/nuu+/47rvvqFmzJlOmTGHu3LnUrVs31c90+vRp2rdvT/fu3enVqxf79++nf//+qe7fqFEjDh06RLVq1ViyZAmNGjWiWLFiKZYJIWxPa/O3dN4809v09GnImxfatTM1xNatIXduq6O8m8PXHLOyQ87VqyYxgvmZFX2A1q5dS4ECBRJerVu3BuDmzZtMnTqVMWPG0KZNG6pUqYK/vz9ubm5MmTIlyTlGjBjBc889R4UKFShZsiReXl5s3rwZMLWwJk2aUL9+/YTaVfIVuPv160fz5s0pV64cTz/9NF9++SULFy4kNv6XAcTExDB58mQaN27M448/TsGCBZkwYQKDBg3ilVdeoXLlykycODFJc2hKpk6dyiOPPMKkSZOoXLkyr7zyCr6+vqnu7+rqSqlSpQAoVqwYpUuXTrVMCGE7R4/Cp59C1apQty6MH28qDHPmQGgoLFgA7dvbZ2KEHFBzTMvcqinJSA0uKAhatIDISHB1hblzM+9hcbynnnqK6dOnJ2zHz/Rz/PhxoqKiaNy4ccJ7zs7OeHp6cvjw4STn8PDwSLLt5eXFlClTiIqKIjAwkGbNmnHz5k0CAwN58cUX2bVrF19++WXC/r/88gujR48mODiYq1evEhMTQ2RkJP/++y8PPfQQAC4uLtSuXTvhmKtXr3L+/Hk8E/2CnJycaNCgAadPn0718wYHB9OwYcMkA/I9M/uXLIRIkzNnTNKbNw/27DFNo08/De+/Dx07QqI+cnbP4ZNjVorvRRUYaLoXZ8Xf7Hz58pHSiiNaayDlWV2Sl+XPnz/JdtOmTYmIiGDXrl1s2bKFfv36cePGDXx8fNi+fTu5cuWifv36AISEhNCmTRu6d+/OyJEjKV68OHv37qVz585ERkYmnDN37tw26YAT/7mEEPYhLMyMRZw3D7ZuNc2oHh4wdiy8+ipk1/5ukhxtzNMza5Li/VSsWBFXV1e2bdtGhQqmR1dMTAxBQUF06dLlnscWKFCAunXrMn36dK5fv07dunWJiorin3/+Ye7cuUmeN+7evZvIyEjGjx+fkPxWrlx53/gKFy5MmTJl+O2332jevDlgEt/OnTuTPA9NrmrVqixZsgStdUKS/+233+7/CxFC2Mz16/DzzyYhrl9vep5WqWKaUV97DSpVsjrCB+fwzxxzqvz589OjRw8GDx7M6tWrCQ4OpkePHly4cIGePXve93gvLy/mzJlD06ZNcXZ2Jk+ePDRo0IA5c+Yked5YqVIlYmNjmTBhAidPnmTevHlMSGObdN++ffnyyy9ZvHgxR44coV+/fpw/f/6ex/j6+nLq1Cn69evHkSNHWLx4Mf7+/mm63v38888/VK5cmRUrVtjkfEI4kvBwCAiAV16BUqXgjTfM+MMBA0yHm0OH4JNPHCMxgiRHhzZmzBheeeUVunXrRu3atdm/fz9r1669Z80sXrNmzYiJiUmSCFMqq1mzJhMnTmTcuHFUrVqVGTNm8PXXX6cpvgEDBtCtWzfee+89GjRoQGxsLF27dr3nMY888ghLly5l7dq11KpVi/Hjx/PFF1+k6Xr3ExkZyZEjR5DZlEROFz8kbetWUzPs1s2MQXzpJfPY6N13zZCMkyfhiy+gVi3rh17Ymsopz3A8PDz07t27U3wvODg4yTg+YV+sXuw4MblXrJW8p7Swve3b73QsjE8PhQqZxNi5MzRvDi528EDOVveCUmqP1tojebkdfEQhhBBW0hp27TI9TWfMgIiIO+917WrK8uSxLj4rOHxyVEq1A9ql1KNTCCFyqvjB+QsWmMm8T56EXLmgfn2TKGNizJC0Xr1yXmKEHJAcMzrOUQghHNGhQyYhLlhgBuo7O8Mzz5jONO3bQ9Gi5pljVg5Js0cOnxyFECKnO3r0TkI8dMhMb+nlZXqavvQSJF+rwF6GpFlJkqMQQjigkydNc+n8+XfmiW7SBCZPNktD3WemxhxPkqMQQjiIM2dMQlywAHbuNGUNGsC4cWbli7JlrY0vO5HkKIQQ2di//5rp2+bPN8MwAOrUgTFjzID9cuUsDS/bkuQohBDZTFgYLFliaohbtphVgKpXh1GjzHymjjJLjZUkOQohRDZw+bKZvm3BArPAQUwMPPEEfPyxSYhVq1odoWOR6eOEJUaMGEH16tXTdUxYWBhKqYR1JU+dOoVSisQzH23fvp2aNWvi6uqaMHtGSmVCZAfXrpn1D9u1M9O3vfsuHDsGH3wAf/wBwcF31kzMSYKCgpg7dy5BQUGZdo1smRyVUhWUUjOVUoutjsVqFy5coG/fvjz22GPkzp0bd3d3WrduzerVq216nYwks8z28MMPc/78+STrRPbt25datWpx/Phxli5dmmqZEPbq5k1TO3zppTsTfO/bB35+ppPN8eNm3tPatR1vPtO02Lp1K08//TQzZ86kRYsWmZYgs7xZVSn1PdAWCNVaV09U3gqYCDgDM7TWqc4mrbU+Abyb05PjqVOnaNy4MQULFmT06NHUqlWL2NhYNm3ahK+vL//880+WxxQZGYmrq2uWXMvZ2ZnSyfqj//333/Tq1YuHH374nmVC2JPwcFizxiTFFSvg1i0z1MLb2ywB1bChGZuY0/3++++8+uqrREVFAebvTWBgYOYseK61ztIX8BRQFziYqMwZOA5UAFyBfUBVoAawMtmrVKLjFqf1uvXq1dOpOXz4cKrvpdeOHTv0//3f/+kdO3bY7Jypad26tS5Tpoy+fv36Xe9dunQp4b+vXLmiu3fvrkuWLKkLFCign3rqKb1r166E92fNmqXz58+vN27cqKtVq6bz5cunvby89IkTJxLeB5K8Zs2apbXWGtDffPON7tChg86XL58eMGCAjo6O1u+8844uV66czpMnj65YsaIeM2aMjomJSbjm8OHDdbVq1e75+Xbu3Knr1q2rc+fOrWvXrq1XrlypAb1582attdYnT57UgN61a1fCfyePMbW4M8qW94pIv/jv3hFs2aL1m29q3bKl1gULag1alyihta+v1ps3ax0dbXWE9uPixYvax8dHK6V0iRIltKurq3ZyctJ58+Z94L+1wG6dQs7I8pqj1vpXpVS5ZMX1gb+1qRGilJoPvKi1Ho2pZWaIUsob8AZwc3NLeFaVXOHChbl+/XqSsg8//JADBw6k63rXrl3j4MGDxMbG4uTkRPXq1SlUqFCaj69RowZjxoxJ076XLl1i7dq1fPLJJ2it74rfxcWF69evo7WmVatWFCpUiAULFlC0aFF++uknmjdvzp49eyhdujTh4eFERETw2WefMXnyZPLkyYOvry/vvfcey5Yt4/nnn6dPnz6sXbs2obm2UKFCCdccMWIEw4cPZ8SIESiluHr1KiVKlGDWrFmUKFGCPXv20LdvX/Lnz8+bb74JQEREBLGxsXfFHe/mzZs8//zzNGnShMmTJ3PhwgX69+8PwK1bt7h+/To3btxI2LdixYocO3aMWrVqMWzYMDp27EiBAgXuKkscd0aEh4eneh+JzHfjxo1s/fuPjFTs2VOMn39+iN9/LwYoQNOw4UU6djxLnTpXcHY2S2Fs3WppqHZBa826deuYNm0a165do1OnTrz99tucPHmSnTt3Ur9+fSIiIjLlnrCX3qruwOlE22eABqntrJQqDnwO1FFKDYlLonfRWk8HpoNZsiq1zhjBwcF3LYnk6uqasLJ9Wl2/fp3Y2FiAhD/8RYsWTfPxrq6uaV6aKTg4GK01tWvXvucxv/zyCwcOHOC///4jb968ANSuXZv169cTEBDAoEGDyJMnD9HR0fj7+/PEE08AMGjQILp160b+/PkpWLAgxYoVw9XVlZQmcH/ttdfo3bt3krLESb569eoEBwcTEBBAr169AMidOzdOTk6pxj5v3jyioqL48ccf0VrTqFEjbt68yRtvvEG+fPkoWLAgBQoUAMzCzkWKFKFIkSIopXBzc0sSZ0plGZUnTx7q1KnzwOcRGZMdl6wKDzdrIi5aBMuXm042iSfydnZWvPBCCQYOLJH6SXKgQ4cO0aNHD7Zu3Yqnpyf+/v7UrFkz4f1q1apl6r1gL8kxpcfKqS40qbW+CPim6cQZXJUjravZJxYUFESLFi0SnrvNnTs3c9rCIb5Z+b727NnDrVu3KFmyZJLy8PBwjh8/nrCdO3fuhMQI8NBDDxEVFcWVK1coVqzYPa/h4XHXUmj4+/szY8YMQkJCuH37NlFRUTz66KNpihlM8q9ZsyYFChRIqOll1u9SCFu7fRvWrjWD81esgOvXzYTenTqZV/780KqVWTPR1dXMcyqMGzduMHLkSMaPH0+hQoWYMWMG3bp1wymLH7raS3I8AyTuLVEWOGeLE+ssXJXD09OTTZs2JfzrNjP/mFeqVAmlFMHBwXTo0CHV/WJjY3Fzc2NrCm00iZt8XZKtXqriusHF14TvJX/+/Em2FyxYQL9+/fj6669p1KgRhQoVYsqUKQQEBNz3XPHSmvyFsBe3bsHq1SYhrlxpep0WL27GIL78MjRrZpaEirdpk6x8kZjWmmXLltG3b19Onz7Nu+++yxdffEGJ5LOiZxF7SY67gEpKqfLAWeA1oIstTpzV6zl6enpmSQ2nWLFitGzZkm+++QY/P7+EJsZ4V65coUiRItStW5cLFy7g5OREhQoVMnw9V1dXYmJi0rTvtm3baNCgQZKm1sS11LSoWrUqP/zwAzdv3kwo++2339J1DiEy240bJiEuWmR+3roFJUvC66+bGqKXF7ik8ldWVr644+TJk/Tp04dVq1ZRo0YN5s2bR+PGjS2NKcs7Byul5gFBwBNKqTNKqXe11tFAb2AdEAws1FofssX1tNYrtNbehQsXtsXp7Mq3336L1hoPDw8WLVrEkSNH+Ouvv5g6dWpC2/wzzzxD48aNefHFF1mzZg0nT54kKCiI4cOHp1ibTE25cuUICQlh7969hIWFEZF4qfBkHn/8cfbu3cuaNWs4duwYo0aNYsuWLen6bF26dMHFxYV33nmH4OBgNmzYwOeff56uc6QmJiaGypUr4+/vb5PziZzl+nWYNw86djTjEF991XSeeest+OUXOHcO/P3NGompJUZhRERE8Pnnn1O1alW2bNnC2LFj2bt3r+WJESxIjlrrzlrrMlrrXFrrslrrmXHlq7XWj2utH9Na2+avIKbmqJSafvXqVVud0m6UL1+evXv38uyzz/Lhhx9Ss2ZNmjdvzvLly5k2bRpgmkdXr15N8+bN6d69O0888QSvvPIKR44c4aGHHkrztTp27Mjzzz9PixYtKFmyJPPmzUt1Xx8fH1555RW6dOnCk08+yalTpxgwYEC6PluBAgVYuXIlx44do2nTpgwcODDNPXnvR2vNkSNHCAsLs8n5hOO7ehXmzjWLAZcsCV26mAWB333XzG169ix8+61pOpWEmDa//PILtWrV4uOPP6Zt27YEBwfTv3//ux7xWCal8R2O+MqqcY7C9q5du2Z1CAnkXrFWVo5zvHxZ6x9+0LptW61dXc04RHd3rfv21XrrVq0TDdsV6XD+/HndpUsXDegKFSroNWvWZOg8troXsJdxjkIIYa8uXYKffzadajZsgKgoePhh6NXLdKpp0EBmqsmomJgYpk6dytChQwkPD2fYsGEMHjw4YYiZvXH45JjVHXKEENlLWJhJiIsWmR6k0dHw6KPQt6/pVFO/fs6cw9SWdu3aha+vb8JjoClTplDJztfVcvjkqLNwKIcQInv47z+z/NPixaYTTUwMVKgA/fubGmK9epIQbeHy5csMHToUf39/SpcuzYIFC3j55ZcThorZM4dPjkIIAXDhgkmIixaZ8YWxsVCxIgwaZGqIdepIQrQVrTVz5sxh4MCBhIWF4efnx8iRI9M1nabVHD45prVZNX4+VCFSk5YJEYT9CAoys9NERMDevfDrryYhPvEEfPSRSYg1a0pCtLXDhw/Ts2dPtmzZQsOGDVm3bl2SZeWyC4dPjmlpVs2fPz9nz57Fzc2NXLlyZYsqv8g6WmuioqK4cOHCXbMBCfsTEgLjxsE335hkCFCuHHz8sWkyrVZNEmJmuHnzJqNGjWLs2LEULFiQadOm8d5772XbSofDJ8e0KFu2LGFhYYSEhBAdHW11OCKZ8PBw8iSeqdkCLi4uFC5c2LKprMS9HT0KS5aY1549Sd9zdjbrIg4ZYk1sOcHy5cvx8/MjJCSEt99+my+//PKu+ZyzG0mOgJOTE6VKlaJUqVJWhyJSEBgYKCthiCS0hgMH7iTEQ3HzadWvD2PGQPnyZsYamdg7c506dQo/Pz9WrFhBtWrV+PXXX2natKnVYdmEwydHGcohhGPQGoKDC7JmDSxdCn//bZpHmzaFiROhQwczJjFe2bIysXdmiYyMZNy4cYwcORInJye++uor+vbtS67EM6tncw6fHGUohxDZV0wMbN9uaodLl8KZM/VwcYEWLeCDD+DFF8HNLeVjZWLvzBEYGEjPnj0TVgSaOHEiDyf+V4mDcPjkKITIXiIjYfNmkwyXLYPQUMid26x/+PrrwQwaVIV0rCEubOTChQsMHDiQOXPmUL58eVauXEmbNm2sDivTSHIUQlju9m1Yv97UEFesgCtXoEABaNMGXnoJnn/ebAcGXqBo0SpWh5ujxMTEMG3aND766CNu3brF0KFD+eijj8iXL5/VoWUqSY5CCEtcv27WQFyyxPy8eROKFjVNpS+9BM89BxZ3Us6xgoKCCAwMxM3NjalTp7J7926aN2/Ot99+yxNPPGF1eFlCkqMQIstcvgzLl5uEuH69GaBfqpRZHPill8ySTw7UpyNbCgoKokWLFoSHh6O1plixYvz000+89tprOWoMuMMnR+mtKoS1Llwwzw6XLDHPEqOjTa9SX1+zYHCjRmYsorCe1ppJkyZx+/ZtwKwH2xQSELcAACAASURBVLt3bzp37mxxZFnP4ZOj9FYVIuv984+Zx3TJEti2zQzDqFgRBgwwCdHDQ2apsTcnT56kV69erFmzBqUUTk5OuLq60qpVK6tDs4TDJ0chROYJCrozlrBECdPDdMkS2LXLvF+9OgwbZhJi9eqSEO1RVFQU48aN49NPP8XZ2ZkJEyZQt25dtm3bhpeXF545dDyMJEchRIbs2GHGG0ZEmG2tzU8PDxg92jxDfPxx6+IT9xcUFISPjw8HDhygffv2TJo0KWHMoqPMdJNRkhyFEGkWGws7d5oa4syZEB5+5702beDbb+GRR6yLT6TNlStXGDJkCNOmTcPd3Z2AgADat29vdVh2RZKjEOKeoqLMck/xg/LPnQMXF6hbF/74wyRMV1cYOlQSo73TWrNw4UL69u3Lf//9R9++fRk5ciQFCxa0OjS7I8lRCHGX27dhwwaTEJcvN0Mw8uaF1q3NHKZt2pgxiYmfOebQR1PZxsmTJ+nZsydr166lXr16rF69mrp161odlt1y+OQoQzmESJurV81g/KVLYc0aMyi/SBFo184kxJYtIfmkKDJ/qf1L3uFm4sSJ9OrVC2cZP3NPDp8cZSiHEKkLDTU1w6VLYeNG04RaujS88YZJiF5epslUZE9BQUF4e3tz8OBBOnTowKRJkyhbtqzVYWULDp8chRBJhYSYMYgBAWYMYmysWf/Qz8/0MG3YELLp4u0izuXLlxM63Dz88MP8/PPPvPDCC1aHla1IchQiBwgONrXDgADYs8eU1agBH39sEmLNmjIG0RForVmwYAH9+vXjv//+4/3332fkyJEUKFDA6tCyHUmOQjggrU0SXLrUvI4cMeUNG8KYMabJtFIla2MUtnXixAl69uzJunXrpMONDUhyFMJBxMSYZtL4GuLp02bOUi8v02T64ovg7m51lMLWoqKiGDt2LJ9++ikuLi5MmjSJnj17SoebByTJUYhsLCLCdKQJCICff4awMLMwcMuWMGoUtG0LxYtbHaXILNu3b8fHx4dDhw7x0ksvMXHiROlwYyOSHIXIZq5fN0MtAgJg1SqzXajQnYWBW7UyCwMLx3X58mUGDx7M9OnTpcNNJpHkKIQdix9kX7s2/PuvaTLdsMHUGEuWhFdfNQmxeXNTYxSOTWvN/Pnz6devH2FhYfTv359PP/1UOtxkgmybHJVS7YE2QClgitZ6vcUhCWFTy5bBK6+YsYfxHnkEevQwHWoaN5Z1EHOS48eP07NnT9avX8+TTz7J2rVrqVOnjtVhOSxLkqNS6nugLRCqta6eqLwVMBFwBmZorb9I7Rxa62XAMqVUUeBrQJKjyPaCg++MQdy9+065UtCzJ0yeLEMucprIyEi+/vprRo0aRa5cuZg8eTI9evSQDjeZzKqa42zgG+B/8QVKKWdgCvAscAbYpZRajkmUo5Md/47WOjTuvz+OO06IbEdrs/ZhfEKMH3JRvz74+sLs2abm6OoKXbtKYsxptm3bho+PD4cPH6Zjx45MnDgRd+lynCWUjl+ELasvrFQ5YGV8zVEp5QmM0Fq3jNseAqC1Tp4Y449XwBfABq31xlT28Qa8Adzc3OrNnz/fxp9CZIUbN2441DOVmBjFvn2F2bq1JNu2lSAsLDdOTprata/QpMl/NGlykZIlzSKJhw4V4s8/i1C79hWqVbtmceTWc7R7ITXXrl1j+vTprFq1Cjc3N/z8/GjUqJHVYdkVW90LzZo126O19khebk/PHN2B04m2zwAN7rF/H+AZoLBSqqLW2j/5Dlrr6cB0AA8PD+3l5WW7aEWWCQwMJLt/d7dvw/r1pna4YgVcumRWuWjZ0jw/bNtWUaxYUaBokuOy+ce2OUe4F+5Fa828efN4//33uXjxIgMHDmT48OE54h8E6ZXZ94I9JceUGoxSrdZqrScBk+57UlmVQ1jk8mUz1CIgANauhVu3kq5y8dxzkD+/1VEKe3H8+HF69OjBhg0bqF+/PuvWraN27dpWh5Vj2VNyPAM8nGi7LHDuQU8qq3KIrHTunBmMHxAAmzdDdDSUKQNvvXVnlYtcuayOUtiT5B1uvvnmG3x9faXDjcXsKTnuAioppcoDZ4HXgC4PelKpOYrMduzYnQ41v/1myipVgv79TUKsX19WuRB3CwoKYvbs2axfv55Tp07RqVMnJkyYIB1u7IRVQznmAV5ACaXUGWC41nqmUqo3sA7TQ/V7rfWhB72W1ByFrWkNf/xxJyEeirtL69Y1U7Z16ABVq0rPUpG6devW0aZNG2JiYlBK8dVXXzFw4ECrwxKJWJIctdadUylfDazO4nCEuK/4Sb0DAszg/JAQUxts2hQmTID27eHRR62OUti7+A433t7exMTEAODk5ERU4pkehF2wp2bVTCHNqiKjwsPvTOq9fPmdSb2ffRaGDTMda0qWtDpKkV0knuGmSpUqnDhxgujoaFxdXR26B2525fDJUZpVRXpcvQqrV5uEuGYN3LhxZ1LvDh3MpN4FC1odpchOUltSaufOnQnDETw9Pa0OUyTj8MlRao7iflatghkz4Px52LvXzEhTqhR07mwSokzqLTIqKCgIb29vDh48SIcOHZg0aVLCklKenp6SFO2YwydHqTmKlJw4YWqHP/wABw6YMqXMKhe9eoGnp0zqLTLu6tWrDBkyBH9/f9zd3Vm2bBkvvvii1WGJdHD45CgEmB6m+/bd6WEanxBLlzZJUWvTwaZmTWjSxNpYRfaltWbJkiX4+flx4cIF/Pz8GDVqFAWlLT7bkeQoHFZMDGzffqeH6alTJhE2bgxjx5oephcuQIsWEBlpJveWfhEio0JCQujVqxerVq2iTp06LF++HA+Pu6bsFNmEwydHeeaYs4SHw6ZNd3qY/vefSXrPPANDh8ILL5jnifEqVDD7BwaaxCiPgER6RUdHM2nSJD755BMAxo4di5+fHy4uDv/n1aE5/LcnzxwdX0o9TAsWTNrDtFCh1I/39JSkKDJm9+7deHt788cff9CmTRumTJnCozLg1SE4fHIUjunff80cpsuWmZpfVBS4uUGXLqa5VHqYisx0/fp1PvnkEyZPnoybmxuLFi2iY8eOKJkWyWFIchTZxt9/m2QYEABBQaYTTYUK4OdnaogNG0oPU5H5fv75Z3r37s3Zs2fx9fVl9OjRFC5c2OqwhI05fHKUZ47Zl9bw55/w/ffl6NMHDh405bVrw4gRJiFWry5zmIqscebMGfr06cOyZcuoXr06CxculHGKDszhk6M8c8xeUp7D9FGaNIHx402TablyVkcpcpKYmBi+/fZbhg4dSlRUFKNHj2bAgAHkkrXHHJrDJ0dh/8LDYcMGkwxTmsO0WLEdtG/f2OowRQ70559/4uPjw86dO3nuueeYOnUqFSpUsDoskQUkOQpLXL1qpm2L72F686bpUdq2rakdJp7DNDBQViwQWevmzZuMGDGC8ePHU7x4cebOnUvnzp2lw00OIslRZImgIFixAqKjzUw1mzebHqalS8Prr5vnh82amTGJQlhp9erV9OzZk5CQEN577z3GjBlDsWLFrA5LZDGHT47SIcdax47BxIkwdSrExpoyd3fo188kxAYNzLRtQljt/Pnz9OvXj4ULF1KlShV+/fVXmjZtanVYwiIOnxylQ07W0tqsbBHfoebQoaTvOztDz57w0UfWxCdEcrGxsUyfPp3BgwcTHh7OyJEjGTRoELlloGyO5vDJUWS+6GjYuvVOQjx92tQGn3rK1Boffhi6dr0zf2mzZlZHLIRx8OBBfHx82LFjB82aNcPf35/HH3/c6rCEHZDkKDLk1i3TwzQgwDxLvHQJ8uSB556DkSNNx5oSJe7sL/OXCnty+/ZtRo0axVdffUXhwoWZPXs2b775pnS4EQkkOYo0u3TpTg/TdetMgixSxCTCDh2gZUvInz/lY2X+UmEvNm7ciK+vL8ePH+ett97i66+/pkTif8kJgSRHcR9nzpim0mXLTM0vJsZ0qOnWzQy5ePppkLHQIjsIDQ1lwIABzJkzh0qVKrFp0yaaN29udVjCTklyFHcJDr7z/HDXLlNWuTIMGmQSooeH9DAV2YfWmlmzZvHBBx8kTBj+0UcfkSdPHqtDE3bM4ZOjDOW4v9hYkwTjE+KRI6a8fn0YPdokxMqVrY1RiPQKCgpi0aJFbN68mT///JMmTZowbdo0qlatanVoIhtw+OQoQzlSFhkJW7aYhPjzz3DuHLi4mA4zfn7w4oum+VSI7GjLli0888wzREdHAzB48GA+//xznKTJQ6SRwydHcceNG6YjTUAArFxppnDLlw9atza1wzZtoGhRq6MU4sFs2bKFTp06JSRGZ2dnChUqJIlRpIskRwcXFmaGWgQEmKEX4eFQvDi89JJJiM8+C3nzWh2lEA/u4sWLfPDBB8yaNYsyZcqQO3duoqOjcXV1xcvLy+rwRDYjydFBBAXdGUdYpsydHqZbt5pnio88Aj4+ZshF48amCVUIR6C1Zu7cufTv359Lly7x4YcfMmzYMPbt20dgYCBeXl6y7qJIN/kT6QB27IAWLSAiwmxrbX7WqAFDh5qEWLu2LAosHM/x48fp0aMHGzZsoEGDBmzcuJGaNWsC4OnpKUlRZJgkx2wqJsbUFpctg1mzTHNpvNatYdIkkA66wlFFRUXx9ddfM3LkSHLlysU333yDr68vzs7OVocmHIQkx2wkIsJMwxYQYBYFDg01c5XWrQt79pjmU1dX+OQTSYzCcQUFBeHt7c3Bgwfp2LEjEydOxF26Vgsbk+Ro565eNYsBBwTA6tWmx2nBgvD886a5tHVrs0hw4meO0pIkHNGVK1cYP348K1asoGzZsixfvpx27dpZHZZwUNkyOSqlqgB9gRLAJq31VItDsql//zVjD5ctMzXFqCgoVQo6dzYJsXlzSL6ajsxdKhyV1prFixfj5+dHaGgoffv2ZeTIkRQsWNDq0IQDS1dyVEoVBsK11hEZvaBS6nugLRCqta6eqLwVMBFwBmZorb9I7Rxa62DAVynlBHyX0VjsybFjJhkGBMBvv5lONY89Bn37miEXDRuatRCFyElCQkLo1asXq1atok6dOowYMQIfHx+rwxI5QJpHxSqlXICLwHMPeM3ZQKtk53YGpgCtgapAZ6VUVaVUDaXUymSvUnHHvABsAzY9YDyW0No8J/z4Y6heHR5/3MxdGhEBn34KBw6YhPnVV2bohSRGkZNER0czbtw4qlatyubNmxk7diw7d+7kiSeesDo0kUMoHd/vPy07K3UW6K61Xv1AF1WqHLAyvuaolPIERmitW8ZtDwHQWo9Ow7lWaa3bpPKeN+AN4ObmVm/+/PkPEvYDi4lR7NtXmG3bSrB9ewlCQ/Pg5KSpWfMKjRuH0aRJGKVLZ7hS7rBu3LhBgQIFrA5DZJEjR44wduxYjh07RsOGDenbty+lS5cG5F4Qd9jqXmjWrNkerbVH8vL0PnOcA7wHPFByTIE7cDrR9hmgQWo7K6W8gJeA3PeKRWs9HZgO4OHhoa2YJePWLVi//s6UbYkXBe7QAdq2VZQoURQoClTK8viyg/iB3MKxXb9+nWHDhjFp0iTc3NxYtGgRHTt2TLIAsdwLIl5m3wvpTY6ngC5KqV3Az8B5IEnVU2v9fQbiSGl4eqpVWq11IBCYphNbsCrHpUsmEcYvCnz7dtoXBRYiJ1q+fDm9e/fmzJkz+Pr6Mnr0aAoXLmx1WCIHS29ynBL30x2ol8L7GshIcjwDPJxouyxwLgPnuTugLFqV4/TpO1O2bdlyZ1Hgd96RRYGFSM3Zs2fx8/Nj6dKlVK9enQULFsisNsIupDc5ls+UKGAXUEkpVR44C7wGdLHFiTOr5qg1HD58p4fpnj2mPH5R4A4doF49WRRYiJTExMTg7+/PkCFDiIqK4v/+7/8YOHAgueRfkMJOpCs5aq1DHvSCSql5gBdQQil1BhiutZ6plOoNrMMM5fhea33oQa8Ftq05bt8Oc+bAtWtmceBjx0x5gwayKLAQabV//368vb35/fffefbZZ5k6dSqPPfaY1WEJkUSGJgFQSlUHngaKYYZ3/Kq1PpiWY7XWnVMpX43tO/rYzJIl0KnTne0nn4Rvv4UXXpBFgYVIi1u3bjFy5EjGjh1L0aJFmTNnDl26dEnS4UYIe5HeSQBcMOMUO5O0E41WSv0EvK21jrFdeA/OVs2qR46YVS20NmMOO3SAHj1sE6MQjm7dunX06NGDkydP8s477/Dll19SvHhxq8MSIlXpfSI2HHgFGIZ5/pg37ucw4NW4n3ZFa71Ca+39oD3fmjUzQzCcnc3k3tKbXIj7u3DhAl26dKFVq1a4uroSGBjIzJkzJTEKu5feZtXXgVFa688TlYUAn8fNctMNk0AdjqenmedUJvcW4v5iY2P5/vvvGTRoEDdv3mT48OEMGTKE3MknBRbCTqU3OT4EBKXy3g5g6IOFY3u27K0qk3sLcX/BwcH4+PiwdetWnn76afz9/aksPdVENpPeZtVzQONU3muEjcYm2pKtmlWFEPcWHh7O8OHDqVWrFgcPHmTmzJls3rxZEqPIltJbc5wLDFVKxcb993mgNGZc4lBgjG3DE0JkB5s3b8bX15ejR4/StWtXxo0bR6lSpawOS4gMS2/NcQSwGPgUOAbcAP4GPk9UbleUUu2UUtOvXr1qdShCOJy1a9dSr149mjdvTnR0NOvXr2fOnDmSGEW2l95JAKIxc6t+DjyFGed4CdiitT6cCfE9sKyaPk6InERrzfDhwxk1ahQALi4uzJgxg2bNmlkcmRC2kebkqJRyxTSb/qS13gXYZAYbIUT28vfff9OjRw82btyYUKa15rfffpPkKBxGmptVtdaRgA9mbKMQIoeJjIzk//7v/6hRowY7d+5k4MCB5M2bF2dnZ1xdXWUpKeFQ0tsh5w+gBvBrJsSSKaxYskoIR7Njxw68vb05dOgQnTp1YuLEiTz00EO89NJLCevqyWoawpGkNzkOAOYppUKAVVrrVNdctBfyzFGIjLty5QpDhgzB39+fhx9+mBUrVtC2bduE9z09PSUpCoeU3uS4CCiMWeg4WikVStJFibXW+lFbBSeEsIbWmsWLF+Pn50doaCjvv/8+I0eOpECBAlaHJkSWSG9y3ETSZCiEcDAhISH06tWLVatWUbduXVauXEm9eimtbS6E40rvUI63MykOIYTFoqOjmThxIsOGDUMpxbhx4+jTpw8uLhla2U6IbC3NvVWVUq5Kqb1KqecyMyBbk0kAhLi/3bt3U79+fQYOHEjz5s05fPgw77//viRGkWOldyhHeSA688KxPZlbVYjUXb9+nX79+tGgQQP+/fdfFi1axPLly3nkkUesDk0IS6V3+rgNQLaqOQohUrZ8+XKqVavGpEmT8PHxITg4mE6dOqGUuv/BQji49LaZTAbmKKVcgGWYiceTdNDRWp+wUWxCiExw9uxZ/Pz8WLp0KdWrV2fBggUyHEOIZNKbHLfE/ewPvJ/KPs4ZD0cIkVliYmLw9/dnyJAhREVFMXr0aAYMGECuXLmsDk0Iu5Pe5NgtU6IQQmSq/fv34+3tze+//86zzz7L1KlTeeyxx6wOSwi7ld6hHD+k9p5SyhkzQYAQwk7cunWLkSNHMnbsWIoWLcqcOXPo0qWLPFcU4j7u2yFHKXVJKVU30bZSSi1XSlVItqsH8J+tA3xQMpRD5FTr16+nevXqjBkzhjfffJPg4GC6du0qiVGINEhLb9UiJK1hOgFt48rtngzlEDlNaGgoXbt2pWXLlri6uhIYGMjMmTMpXry41aEJkW2kdyiHEMJOxcbGMnPmTCpXrszixYsZPnw4+/bt4+mnn7Y6NCGyHZn+QggHEBwcjI+PD1u3buWpp55i2rRpVK5c2eqwhMi2pOYoRDYWHh7O8OHDqVWrFgcPHmTGjBls3rxZEqMQDyitNUf3RB1wnBOVXUm0T1nbhSWEuJ/AwEB8fHw4evQoXbp0Yfz48ZQqVcrqsIRwCGlNjotTKFuWbFshy1kJkekuXrzIBx98wKxZsyhfvjxr166lZcuWVoclhENJS3KUgf9C2AGtNXPnzqV///5cvnyZwYMH88knn5AvXz6rQxPC4dw3Od5r4L8QInMFBQURGBhIxYoV+e6779iwYQMNGzZk+vTp1KhRw+rwhHBY0ltVCDsVFBREixYtCA8PR2tNvnz5mDJlCj4+Pjg7yxTGQmSmbNtbVSmVXym1RynV1upYhMgM//vf/7h9+zZam0f5fn5+9OzZUxKjEFkgy5OjUup7pVSoUupgsvJWSqkjSqm/lVKD03CqD4GFmROlENa5evUqPXv2xN/fH6UUTk5O5M2blxdeeMHq0ITIMaxoVp0NfAP8L74gbtLyKcCzwBlgl1JqOWbYyOhkx78D1AQOA3myIF4hsoTWmiVLluDn58eFCxfo168fbdq0YdeuXXh5ecmai0JkIRXfZJOlF1WqHLBSa109btsTGKG1bhm3PQRAa508McYf/zmQH6gK3AY6aK1jU9jPG/AGcHNzqzd//nybfxaR+W7cuEGBAgWsDiNT/fvvv0ycOJHffvuNSpUqMWDAAJ544gmrw7I7OeFeEGljq3uhWbNme7TWHsnL7aVDjjtwOtH2GaBBajtrrYcCKKXeBsJSSoxx+00HpgN4eHhoLy8vG4UrslJgYCCO+t1FR0czadIkhg0bhtaasWPH4ufnh4uLvfyvaV8c+V4Q6ZPZ94K9/B+Y0ho6963Saq1n3/fESrUD2lWsWDEDYQmRefbs2YO3tzd79+6lTZs2TJkyhUcffdTqsIQQ2E9v1TPAw4m2ywLnbHFiWbJK2JsbN27w/vvvU79+fc6dO8fChQtZsWKFJEYh7Ii9JMddQCWlVHmllCvwGrDcFieWxY6FPVmxYgVVq1ZlwoQJeHt7ExwczMsvvywLEAthZ6wYyjEPCAKeUEqdUUq9q7WOBnoD64BgYKHW+pAtric1R2EPzp07R6dOnXjhhRcoVKgQ27dvZ+rUqRQpki3WDBcix8nyZ45a686plK8GVtv6evLMUVgpJiaGadOmMWTIECIjI/n8888ZOHAgrq6uVocmhLgHe2lWzTRScxRW2b9/P40bN6ZXr17Ur1+fAwcO8NFHH0liFCIbcPjkKERWu3XrFoMHD6ZevXocP36cH3/8kfXr1yOtF0JkH/YylCPTSLOqyErr16+nR48enDhxgm7duvHVV19RvHhxq8MSQqSTw9ccpVlVZIXQ0FC6du1Ky5YtcXZ25pdffuH777+XxChENuXwyVGIzKS1ZubMmVSuXJlFixYxbNgw9u/fT7NmzawOTQjxAKRZVYgM+uuvv/Dx8eHXX3+ladOmTJs2jSpVqlgdlhDCBhy+5ijNqsLWwsPDGTFiBLVq1WL//v189913BAYGSmIUwoE4fM1RCFsKDAzEx8eHo0eP0rlzZ8aPH4+bm5vVYQkhbMzha45C2MLFixd55513aNasGVFRUaxZs4affvpJEqMQDsrhk6PMrSoehNaaOXPmUKVKFf73v/8xaNAgDh48SKtWrawOTQiRiRw+OcozR5FRx48fp2XLlrzxxhuUL1+ePXv2MGbMGPLly2d1aEKITObwyVGI9IqKimL06NFUr16d3377jcmTJ7Njxw5q1apldWhCiCwiHXKESCQoKAhvb28OHjxIhw4dmDx5Mu7u7laHJYTIYpIcRY4XFBTEmjVrOHjwIMuWLcPd3Z1ly5bx4osvWh2aEMIiDp8cZRIAcS87duygWbNmREZGAvDyyy8zc+ZMChYsaHFkQggrOfwzR+mQI1ITEhLCO++8k5AYnZ2dqVOnjiRGIYTjJ0chkouOjmbcuHFUq1aNU6dOkStXLpydnXF1dcXLy8vq8IQQdsDhm1WFSGzPnj14e3uzd+9enn/+eaZMmcL58+cJDAzEy8sLT09Pq0MUQtgBSY4iR7hx4waffPIJkyZNomTJkixYsICXX34ZpRTlypWTpCiESEKSo3B4K1eupFevXvzzzz/4+PjwxRdfUKRIEavDEkLYMXnmKBzWuXPnePnll2nXrh0FChRg27Zt+Pv7S2IUQtyXwydHmVs154mNjeXbb7+lSpUqrFixgs8++4w//viDxo0bWx2aECKbcPjkKEM5cpYDBw7QuHFjevXqhYeHBwcOHGDo0KG4urpaHZoQIhtx+OQocobbt28zZMgQ6taty7Fjx/jhhx/YuHEjlSpVsjo0IUQ2JB1yRLa3YcMGfH19OXHiBG+//TZfffUVJUqUsDosIUQ2JjVHkW2Fhoby+uuv89xzz+Hs7Mwvv/zCrFmzJDEKIR6YJEeR7Wit+f7776lSpQoLFy7kk08+Yf/+/TRr1szq0IQQDkKaVUW2cuTIEXx8fNiyZQtNmjRh2rRpVK1a1eqwhBAORmqOIluIiIjg008/pWbNmuzbt4/p06ezZcsWSYxCiEwhNUdh9/bt20ePHj3466+/eO211xg/fjylS5e2OiwhhAOT5Cjs1qVLlxg0aBAzZ86kXLlyrFmzhlatWlkdlhAiB8iWzapKKS+l1FallL9SysvqeIRtaa356aefqFKlCrNnz+bVV1/l4MGDkhiFEFkmy5OjUup7pVSoUupgsvJWSqkjSqm/lVKD73MaDdwA8gBnMitWkfVOnDhBq1at6Nq1K48++ii7d+/G19eX/PnzWx2aECIHsaLmOBtIUgVQSjkDU4DWQFWgs1KqqlKqhlJqZbJXKWCr1ro18CHwaRbHLzJBVFQUY8aMoXr16gQFBTF58mSCgoKoXbu21aEJIXKgLH/mqLX+VSlVLllxfeBvrfUJAKXUfOBFrfVooO09TncZyJ3am0opb8AbwM3NjcDAwIwHLjLN4cOHGTt2LCdOnKBp06b06dOHkiVLsnXrVsCsxSjfnQC5F8QdmX0v2EuHHHfgdKLtM0CD1HZWSr0EtASKAN+ktp/WejowHcDDw0N7AoqiWgAADpJJREFUeXnZIlZhI9euXeOjjz7i22+/5aGHHmLZsmW8+OKLd+0XGBiIfHcC5F4Qd2T2vWAvyVGlUKZT21lrvRRYmqYTK9UOaFexYsUMhiZsTWtNQEAAffr04fz58/Tp04fPPvuMggULWh2aEEIA9tNb9QzwcKLtssA5W5xYlqyyL6dPn6Z9+/Z07NiRkiVL8vvvvzNx4kRJjEIIu2IvyXEXUEkpVV4p5Qq8Biy3xYllsWP7sG3bNtq2bcvjjz/Oxo0b+eqrr9i9ezdPPvmk1aEJIcRdrBjKMQ8IAp5QSp1RSr2rtY4GegPrgGBgodb6kC2uJzVH682ePZunnnqKVatWERkZyY8//sjAgQNxcbGXVn0hhEjKit6qnVMpXw2stvX15JmjdW7cuMHw4cMZP348WptHyEopjhw5YnFkQghxb/bSrJpppOZojVWrVlGtWjXGjRvHCy+8QN68eXF2dsbV1VV6Gwoh7J60awmbOn/+PH379mXRokVUrVqVrVu30qRJE4KCghK6Xnt6elodphBC3JPDJ0dpVs0asbGxTJ8+ncGDBxMeHs5nn33GBx98gKurKwCenp6SFIUQ2YY0q4oHdvDgQZo2bUqPHj2oV68eBw4cYOjQoQmJUQghshuHT44i89y+fZuhQ4dSp04djhw5wg8//MDGjRupVKmS1aEJIcQDkWZVkSGbNm3C19eXv//+m7feeouvv/6aEiVKWB2WEELYhMPXHKVZ1bb+++8/3nzzTZ555hnAJMnZs2dLYhRCOBSHT47CNrTWzJo1i8qVKzN//nw+/vhjDhw4QPPmza0OTQghbM7hm1XFgzt69Cg+Pj4EBgbSuHFjpk+fTtWqVa0OSwghMo3D1xxlbtWMi4iIYOTIkdSoUYM//viDadOm8euvv0piFEI4PIdPjvLMMWO2bt1KnTp1GD58OB06dOCvv/7C29sbJyeHv2WEEMLxk6NIn8uXL9O9e3eeeuopbt26xerVq5k/fz6lS5e2OjQhhMgykhwFYDrczJs3j8qVKzNr1iwGDhzIoUOHaN26tdWhCSFElpMOOYKTJ0/So0cP1q1bx5NPPsm6deuoXbu21WEJIYRlHL7mKB1yUhcVFcWXX35JtWrV2L59O5MmTSIoKEgSoxAix3P45CgdclL2+++/4+HhwYcffkjLli0JDg6mT58+ODs7Wx2aEEJYzuGTo0jq2rVr9O7dG09PTy5evMjSpUsJCAigbNmyVocmxP+3d/+xVdVnHMffDzWlZBkEFY1DMmnmkKKMJhekLF3IggOdkx8jCiEMGKOtAhI3dS4yCGHEQXRBixWKU0BTwChUQRBcBasLBhokKlRZA04aQ4QFUcEMvDz7o1WasxZoufeecw+fV3L/uN9zzvc+KQ/99Ht6eo5IZCgcLyHr16+nb9++VFRUMH36dPbt28fo0aPDLktEJHIUjpeAQ4cOMWrUKMaMGUOPHj145513KC8vp2vXrmGXJiISSQrHGEsmkzz++OMUFBSwdetWFi1axK5duxg0aFDYpYmIRJr+lCOm3n33XUpKSqirq2PEiBFUVFTQu3fvsMsSEckKWjnGzIkTJ7j//vsZOHAgn3zyCatXr2bTpk0KRhGRdoj9yvFSedjxjh07WLZsGVu2bOHw4cNMmzaNhQsX0r1797BLExHJOrEPR3ffAGxIJBLTwq4lXTZu3MioUaNIJpOYGRUVFdx9991hlyUikrV0WjWLnTlzhmXLljF27FiSySQAnTp14vPPPw+5MhGR7KZwzFJ79+6luLiYsrIy+vXrR15eHjk5OeTm5jJ06NCwyxMRyWoKxyzz9ddfM3v2bAoLC/nwww959tlnqaur44033mD+/PnU1NRQVFQUdpkiIlkt9r9zjJOamhrKyspoaGhg4sSJPPbYY/To0QOAoqIihaKISIpo5ZgFjh49yqRJkxg2bBjuzuuvv86qVau+C0YREUkthWOEuTsrV67khhtuoKqqiocffpj333+fYcOGhV2aiEis6bRqRO3fv5+ysjK2bdvGkCFDqKyspF+/fmGXJSJyScjKlaOZdTKzBWZWbmaTwq4nlU6dOsX8+fPp378/u3fvZunSpbz11lsKRhGRDMp4OJrZM2b2mZl9EBgfYWYfmVmDmT10nmlGAj2B00BjumrNtLfffpsBAwYwZ84cRo4cSX19PaWlpXTqlJU/w4iIZK0wvuuuAEa0HDCzHOBJ4FagABhvZgVmdpOZbQy8rgL6ADvc/fdA1t8K5tixY5SUlFBcXMzJkyd59dVXWbt2Lddcc03YpYmIXJIy/jtHd681s+sCw4OABnc/AGBma4CR7v4IcHtwDjNrBE41v0229VlmVgKUAFx99dVs3779YstPKXdn27ZtLFmyhOPHj3PnnXcyefJkunTpErlaw/TVV1/p6yGAekHOSncvROWCnJ7AoRbvG4Gbz7H/OqDczIqB2rZ2cvdKoBIgkUh4lO4cc/DgQe655x5ee+01EokElZWVFBYWhl1WJG3fvl13/RFAvSBnpbsXohKO1sqYt7Wzu58Epl7QxBF7Ksfp06dZvHgxc+fOJScnh8WLFzNjxgxycnLCLk1ERJpF5UqPRqBXi/fXAp+mYmJ33+DuJd26dUvFdBdl586dDBw4kAcffJBbbrmFffv2MWvWLAWjiEjERCUcdwHXm1lvM8sFxgGvpGJiM/uVmVUeP348FdN1yBdffMG9997L4MGDOXLkCOvWraO6uppevXqd/2AREcm4MP6UYzWwA+hjZo1mNtXdvwFmAFuAeuAFd9+bis8Le+VYXV1NQUEBS5YsYfr06dTX1zN69GjMWjuTLCIiURDG1arj2xjfBGxK9eeF9TvHxsZGZs6cSXV1Nf379+ell17i5pvPdY2RiIhERVROq6ZNpleOyWSS8vJyCgoK2LJlCwsXLqSurk7BKCKSRaJytWos7Nmzh9LSUnbu3Mnw4cOpqKggPz8/7LJERKSdYr9yzMQFOSdOnOCBBx4gkUjw8ccfU1VVxebNmxWMIiJZKvbhmO7Tqps3b+bGG2/k0UcfZcqUKdTX1zN+/HhdcCMiksViH47pcvjwYcaNG8dtt91GXl4etbW1LF++nMsvvzzs0kRE5CLFPhxTfVr1zJkzVFZW0rdvX9avX8+8efPYs2cPxcXFKZlfRETCF/twTOVp1aqqKvLz8yktLWXAgAG89957zJkzh86dO6egUhERiQpdrXqB1qxZw4QJEwDIzc1lwYIF9OnTJ+SqREQkHWK/ckyVgwcPfneRTTKZ5M033wy5IhERSZfYh2Oqfuc4dOhQ8vLyyMnJITc3V4/NERGJsdifVnX3DcCGRCIx7WLmKSoqoqam5rtniBUVFaWoQhERiZrYh2MqFRUVKRRFRC4BsT+tKiIi0l4KRxERkQCFo4iISEDswzETNx4XEZF4iX04Zvp5jiIikv1iH44iIiLtpXAUEREJUDiKiIgEmLuHXUNGmNkR4N8pmq4bkK4rfC527o4c355jzrfvxWxva9uVwNELqi7zotwLHZlDvdBx6oX0bE93L/zQ3Xv836i769XOF1AZ1bk7cnx7jjnfvhezva1tQF3Y/+bZ2AsdmUO9oF5QLzS9dFq1YzZEeO6OHN+eY86378VsT+fXNV2i3AsdmUO90HHqhfRsD6UXLpnTqpK9zKzO3RNh1yHhUy/It9LdC1o5SjaoDLsAiQz1gnwrrb2glaOIiEiAVo4iIiIBCkcREZEAhaOIiEiAwlFERCRA4ShZy8xGmdlyM3vZzH4Rdj0SHjPLN7O/m9mLYdcimWdm3zOzlc3fDyakYk6Fo4TCzJ4xs8/M7IPA+Agz+8jMGszsoXPN4e7V7j4NmAzclcZyJY1S1AsH3H1qeiuVTGpnX4wBXmz+fnBHKj5f4ShhWQGMaDlgZjnAk8CtQAEw3swKzOwmM9sYeF3V4tDZzcdJdlpB6npB4mMFF9gXwLXAoebdkqn48MtSMYlIe7l7rZldFxgeBDS4+wEAM1sDjHT3R4Dbg3OYmQF/BTa7++70VizpkopekPhpT18AjTQF5B5StOjTylGipCdnf/qDpobveY79ZwLDgLFmVpbOwiTj2tULZnaFmS0FCs3sT+kuTkLTVl+sA35tZk+RonuxauUoUWKtjLV5Cyd3fwJ4In3lSIja2wv/AfQDUvy12hfufgKYksoP0spRoqQR6NXi/bXApyHVIuFSL0hrMtYXCkeJkl3A9WbW28xygXHAKyHXJOFQL0hrMtYXCkcJhZmtBnYAfcys0cymuvs3wAxgC1APvODue8OsU9JPvSCtCbsv9FQOERGRAK0cRUREAhSOIiIiAQpHERGRAIWjiIhIgMJRREQkQOEoIiISoHAUEREJUDiKxJSZPW1mbmZ/C7sWkWyjmwCIxJCZdQEOA98HjgA9m+8uIiIXQCtHkXgaDXQFFgFXEXhorIicm8JRJJ4mAQeBP9O0cvxNy41m9iMzO21m8wLjT5nZl2aWyFypItGjcBSJGTP7AU0PgX7e3U8Da4A7zKz7t/u4ewPwNHCfmV3ZfNwc4LfAaHevy3zlItGhcBSJn4k0/d9+vvn9KqAzcFdgv3lADvBHM5sKzAUmuvs/MlWoSFTpghyRmDGzvcCX7j64xVg9cMzdhwT2XQD8AbgMmOXuT2a0WJGI0spRJEbMbCBQADwX2PQcUGRmPw6M/4umVeUOBaPIWQpHkXiZBJwG1gbGnwecFhfmmNnPgWU0PVD2p2b2k0wVKRJ1Oq0qEhNmlgt8CvzT3Ue2sn0bkA9cBxQC22laUd4H7Af2uvsvM1WvSJRdFnYBIpIytwNXAIfMbFQr2w8AQ4HfAX8BtgIz3f1M8590PGNmP3P32kwVLBJVWjmKxISZvQzccQG7OlALDHf3/zYfmwN8QCsX7YhcihSOIiIiAbogR0REJEDhKCIiEqBwFBERCVA4ioiIBCgcRUREAhSOIiIiAQpHERGRAIWjiIhIwP8ALfQUjTGxAH4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import math\n", "\n", "def forward_diff(f, x, dx):\n", " fx = f(x)\n", " fxph = f(x + dx)\n", " return (fxph - fx) / dx\n", "\n", "\n", "def central_diff(f, x, dx):\n", " fxph = f(x + dx)\n", " fxnh = f(x - dx)\n", " return (fxph - fxnh) / (2 * dx)\n", "\n", "\n", "# for this example we know trivially what the exact solution should be\n", "exact = np.cos(0.8)\n", "\n", "print('Exact derivative at sin(0.8) = ', exact)\n", "# headers for the following errors outputs\n", "print('%20s%40s' % ('Forward differencing', 'Central differencing'))\n", "\n", "# we're going to store all the values for plotting, initialise variable for these\n", "fd_errors = []\n", "cd_errors = []\n", "dx_all = []\n", "dx = 1.0 # an initial mesh spacing\n", "for i in range(10):\n", " fd = forward_diff(np.sin, 0.8, dx)\n", " cd = central_diff(np.sin, 0.8, dx)\n", " print('%10g (error=%10.2g) %10g (error=%10.2g)' %\n", " (fd, abs(fd - exact), cd, abs(cd - exact)))\n", " # store the h and the errors\n", " dx_all.append(dx)\n", " fd_errors.append(abs(fd - exact))\n", " cd_errors.append(abs(cd - exact))\n", " dx = dx / 2 # halve h for the next iteration\n", "\n", "# as we expect a polynomial relationship between dx and the errors,\n", "# a log-log plot will demonstrate this if we get straight lines\n", "# the slopes of these lines indicating the order of the relationship:\n", "# slope 1 for forward diff and slope 2 for central diff\n", "\n", "# set up figure\n", "fig = plt.figure(figsize=(7, 5))\n", "ax1 = plt.subplot(111)\n", "\n", "ax1.loglog(dx_all, fd_errors, 'b.-', label='Forward diff.')\n", "ax1.loglog(dx_all, cd_errors, 'k.-', label='Central diff.')\n", "ax1.set_xlabel('$\\Delta x$', fontsize=16)\n", "ax1.set_ylabel('Error', fontsize=16)\n", "ax1.set_title('Convergence plot', fontsize=16)\n", "ax1.grid(True)\n", "ax1.legend(loc='best', fontsize=14)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2.4: Compute second derivative\n", "\n", "Calculate the second derivative $f''$ at $x = 1$ using the data below:\n", "\n", "$f(0.84) = 0.431711$\n", "\n", "$f(0.92) = 0.398519$\n", "\n", "$f(1.00) = 0.367879$\n", "\n", "$f(1.08) = 0.339596$\n", "\n", "$f(1.16) = 0.313486$\n", "\n", "You should get 0.36828" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.36828124999999967\n" ] } ], "source": [ "dx = 0.08\n", "ddf = (0.339596 - 2*0.367879 + 0.398519)/(dx*dx)\n", "print(ddf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2.5: Implementing Forward Euler's method\n", "\n", "Write a function *euler*( *f*, *u0*, *t0*, *t_max*, *h*) that takes as arguments the function $f(u,t)$ on the RHS of our ODE,\n", "an initial value for $u$, the start and end time of the integration, and the time step.\n", "\n", "Use it to integrate the following ODE problems up to time $t=10$\n", "\n", "$$u'(t)=u(t),\\quad u(0)=1$$\n", "\n", "and \n", "\n", "$$u'(t)=\\cos(t),\\quad u(0)=0$$\n", "\n", "and plot the results. A template to get you started is below." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAFDCAYAAABlbTTTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhc9X3v8fdXm2XZknfLtmzwCsZ2WIwKhgSqQMoScjFpkqfOzeInoXWa0pClvSlpe0ubNL1pmzYJSZt7HSCQlIZQQouTkBAHoxAabDDGeF/kXZYsyZYlS7LWme/9Y47MWB7Z1nrOjD6v55lnZn7nd85854fxx+ec35xj7o6IiIj0XVbYBYiIiKQrhaiIiEg/KURFRET6SSEqIiLSTwpRERGRflKIioiI9FNO2AVEzeTJk3327NkD2kZLSwtjxowZnIIyjMYmNY1LahqX3mlsUhuKcXn99dePu/uUVMsUoj3Mnj2bjRs3Dmgb5eXllJWVDU5BGUZjk5rGJTWNS+80NqkNxbiY2aHelulwroiISD8pREVERPpJISoiItJPClEREZF+UoiKiIj0k0JURESknxSiIiIi/aQQFRER6SeFqIiISD8pREVEJOM0tnby/VcOcrShdUg/RyEqIiIZZ/exJv73s9vZW9M0pJ+jEBURkYyzr64ZgHlTxg7p5yhERUQk4+yrbSY/N4uS8aOH9HMUoiIiknH21TUzd/JYsrJsSD9HISoiIhlnX10L86YO7aFcUIiKiEiGaeuMceTkaeZNGfqblitERUQkoxw80YL70E8qAoWoiIhkmH21LYBCVEREpM/21TVjBnMm63CuiIhIn1TUNlMyfjSj87KH/LMUoiIiklH21TUPy6FcUIiKiEgGibuzv65FISoiItJXJ9uc1s4Y86YO/flQUIiKiEgGqW6JA8MzMxcUoiIikkGqmx0YASFqZo+aWa2ZbUux7E/NzM1scvDezOwhM6swsy1mtjSp70oz2xs8Via1X2tmW4N1HjKzob2AooiIhK66JU5Rfg6Tx+YNy+eFuSf6GHBHz0YzmwX8DnA4qflOYEHwWAV8O+g7EXgQuB64DnjQzCYE63w76Nu93jmfJSIimaW6Jc68qWMZrv2m0ELU3V8C6lMs+hrwecCT2pYD3/OE9cB4M5sO3A6sdfd6dz8JrAXuCJYVufsr7u7A94B7hvL7iIhI+KpbfNgO5ULEzoma2d3AUXd/s8eiEuBI0vvKoO187ZUp2kVEJEOdauukod2ZPwx3b+mWM2yfdAFmVgD8BXBbqsUp2rwf7b199ioSh34pLi6mvLz8QuWeV3Nz84C3kak0NqlpXFLTuPROY3Ou/Q0xAFprDlBefuQCvQdHZEIUmAfMAd4MjmXPBDaZ2XUk9iRnJfWdCVQF7WU92suD9pkp+qfk7quB1QClpaVeVlbWW9eLUl5ezkC3kak0NqlpXFLTuPROY3Ou469Xwvo3Wf7OZczN9Nm5Pbn7Vnef6u6z3X02iSBc6u7HgDXAR4NZusuARnevBp4HbjOzCcGEotuA54NlTWa2LJiV+1Hg2VC+mIiIDIt9dc1kG8yaWDBsnxnmT1x+ALwCXG5mlWZ273m6PwfsByqA7wB/BODu9cCXgNeCxxeDNoBPAg8H6+wDfjYU30NERKJhX20zUwuM3Ozhi7bQDue6+wcvsHx20msH7uul36PAoynaNwJLBlaliIiki721zZSMHd59w8gczhUREemvts4YB0+0KERFRET6qqK2GXeYWagQFRER6ZPdx5oAtCcqIiLSV3tqmsjLzqK4YHgvk64QFRGRtLe7pol5U8eSnaUQFRER6ZM9x5q4vHj4LvfXTSEqIiJp7VRbJ1WNbVw2rXDYP1shKiIiaW1vTWJS0eXFClEREZE+2X2sGYDLFKIiIiJ9s6emiTF52cycMHrYP1shKiIiaW33sSYum1ZIcAewYaUQFRGRtLanpimU86GgEBURkTR2vLmdEy0doZwPBYWoiIiksT3B5f4uD+HnLaAQFRGRNLYrCFHtiYqIiPTRnpomJo7JY/LYvFA+XyEqIiJpa3dNE5cVjw1lZi4oREVEJE25e3DN3HAO5YJCVERE0tTRhlZaOmKhXDO3m0JURETS0q7q8K6Z200hKiIiaWlH9SnMYOH0otBqUIiKiEha2l7VyOxJYxg7Kie0GhSiIiKSlnZUn2LRjPD2QkEhKiIiaaixtZMj9a0sCvFQLoQYomb2qJnVmtm2pLZ/NLNdZrbFzP7TzMYnLfuCmVWY2W4zuz2p/Y6grcLMHkhqn2NmG8xsr5n90MzC+SWuiIgMup3VpwBYPIL3RB8D7ujRthZY4u5XAnuALwCY2SJgBbA4WOdfzSzbzLKBfwHuBBYBHwz6Avw98DV3XwCcBO4d2q8jIiLDZXtVIkRH7OFcd38JqO/R9gt37wrergdmBq+XA0+6e7u7HwAqgOuCR4W773f3DuBJYLklLl1xC/B0sP7jwD1D+oVERGTY7Kg6xZTCUUwtzA+1jiifE/048LPgdQlwJGlZZdDWW/skoCEpkLvbRUQkA2yvagz9fChAePOCz8PM/gLoAp7obkrRzUn9jwA/T//ePm8VsAqguLiY8vLyvpR7jubm5gFvI1NpbFLTuKSmcendSB6bzrizt+Y08wvazhmD4R6XyIWoma0E3gPc6u7dwVcJzErqNhOoCl6naj8OjDeznGBvNLn/Odx9NbAaoLS01MvKygb0HcrLyxnoNjKVxiY1jUtqGpfejeSx2Xa0kdgvXuaOZUsou3LGWcuGe1widTjXzO4A/gy4291PJy1aA6wws1FmNgdYALwKvAYsCGbi5pGYfLQmCN8XgfcH668Enh2u7yEiIkNnR1X3zNxxIVcS7k9cfgC8AlxuZpVmdi/wLaAQWGtmm83s/wK4+3bgKWAH8HPgPnePBXuZfww8D+wEngr6QiKMP2dmFSTOkT4yjF9PRESGyPaqRsbkZXPpxIKwSwnvcK67fzBFc69B5+5fBr6cov054LkU7ftJzN4VEZEMsqP6FFdMLyIrK5x7iCaL1OFcERGR84nHnR1V4V/ur5tCVERE0sbh+tO0dMRCv1JRN4WoiIikjTNXKpoe/qQiUIiKiEga2VHdSE6WsaB4bNilAApRERFJI9urTjF/6ljyc7PDLgVQiIqISJpwd7ZUNrKkJBqHckEhKiIiaaLyZCv1LR1cNWv8hTsPE4WoiIikhS2VjQBcNVN7oiIiIn3yZmUDedlZLJwWjZ+3gEJURETSxJtHGrhiRhF5OdGJruhUIiIi0otY3Nl2tDFSh3JBISoiImlgf10zLR0xrpwZnUlFoBAVEZE0sPlIAwBXz9KeqIiISJ9sqWxk7Kgc5k6OxpWKuilERUQk8rZUNrCkJBq3P0umEBURkUhr74qxo/pUpC6y0E0hKiIikbaruonOmHNVxCYVgUJUREQibktlYlLRlRH7eQsoREVEJOI2H2lk8tg8SsaPDruUcyhERUQk0rZUNnDlzPGYRWtSEShERUQkwprbu6ioa47koVxQiIqISIRtrWzEnUhOKgKFqIiIRNimwycBuOYShaiIiEifbDxYz/ypYxlfkBd2KSmFFqJm9qiZ1ZrZtqS2iWa21sz2Bs8TgnYzs4fMrMLMtpjZ0qR1Vgb995rZyqT2a81sa7DOQxbFM9IiItKreNzZdLiB0ksnhF1Kr8LcE30MuKNH2wPAC+6+AHgheA9wJ7AgeKwCvg2J0AUeBK4HrgMe7A7eoM+qpPV6fpaIiETYvrpmGls7uVYhei53fwmo79G8HHg8eP04cE9S+/c8YT0w3symA7cDa9293t1PAmuBO4JlRe7+irs78L2kbYmISBrYeChxPjTKIZoTdgE9FLt7NYC7V5vZ1KC9BDiS1K8yaDtfe2WK9pTMbBWJvVaKi4spLy8f0Jdobm4e8DYylcYmNY1LahqX3o2Esfnp1nYK8+DQttc4fJFn5IZ7XKIWor1JNXrej/aU3H01sBqgtLTUy8rK+lHiW8rLyxnoNjKVxiY1jUtqGpfejYSx+ZuN5SybP4l3vrP0otcZ7nGJ2uzcmuBQLMFzbdBeCcxK6jcTqLpA+8wU7SIikgaON7dz4HhLpA/lQvRCdA3QPcN2JfBsUvtHg1m6y4DG4LDv88BtZjYhmFB0G/B8sKzJzJYFs3I/mrQtERGJuE3B+dAoz8yFEA/nmtkPgDJgsplVkphl+xXgKTO7FzgMfCDo/hzwbqACOA18DMDd683sS8BrQb8vunv3ZKVPkpgBPBr4WfAQEZE08Pqhk+RlZ7GkJJqX++sWWoi6+wd7WXRrir4O3NfLdh4FHk3RvhFYMpAaRUQkHBsPnWRJSRH5udlhl3JeUTucKyIiI1x7V4ytlY2Uzp4YdikXpBAVEZFI2Xa0kY5YPPKTikAhKiIiEbPxYGJS0dJLFKIiIiJ98vqhk8yeVMCUwlFhl3JBClEREYmMeNx59WB9WpwPBYWoiIhEyO6aJhpOd3LD3Elhl3JRFKIiIhIZr+w7AcCyeQpRERGRPlm//wSXTCygZPzosEu5KApRERGJhHjc2XCgnmVz0+N8KChERUQkInYeO0Vjayc3pMmhXFCIiohIRJw5H5omk4pAISoiIhGxfn89sycVMH1cepwPBYWoiIhEQCzubDhwIq0O5YJCVEREImBn9Sma2rrS6lAuKERFRCQC0vF8KChERUQkAtbvP8HcyWMoLsoPu5Q+UYiKiEioumJxXj1QnzZXKUqmEBURkVBtrzpFU3v6nQ8FhaiIiITs5YrjAGlz0flkClEREQnVr/bUsXhGUVrcP7QnhaiIiISmqa2TTYdOctOCKWGX0i8KURERCc0r+07QFXduvmxy2KX0S05fOpvZMuAOYBkwAxgNHAd2A78C/svdTw52kSIikple2ltHQV42pZemz51bkl3UnqiZrTSzrcBvgM8ABcBeYANwErgeeBg4amaPmdmcgRRlZp81s+1mts3MfmBm+WY2x8w2mNleM/uhmeUFfUcF7yuC5bOTtvOFoH23md0+kJpERGTwvbTnODfMnUReTnoeGL1g1Wb2JvAV4DngWmCCu9/s7u9z9w+7+7vd/QpgIvAHwFRgu5n9Xn8KMrMS4H6g1N2XANnACuDvga+5+wISwX1vsMq9wEl3nw98LeiHmS0K1ltMYu/5X80suz81iYjI4Dt4vIXD9ae5+bL0PB8KF7cn+l1gjrv/mbu/4e6eqpO7N7r7E+7+buAGoGEAdeUAo80sh8RebzVwC/B0sPxx4J7g9fLgPcHyW83MgvYn3b3d3Q8AFcB1A6hJREQG0Ut76wAyO0Td/evu3gZgZk+Z2ayLWOdNd3++PwW5+1Hgq8BhEuHZCLwONLh7V9CtEigJXpcAR4J1u4L+k5LbU6wjIiIhe2lPHbMmjmb2pIKwS+m3Pk0sAt5PIuCO9FxgZhOBhe7+m4EUZGYTSOxFziGxN/sfwJ0punbvEVsvy3prT/WZq4BVAMXFxZSXl/et6B6am5sHvI1MpbFJTeOSmsald+k+Nl1x5+U9p7lhRg6/+tWvBm27wz0uFwxRM7s86LfzAl0XAL8mcQ5zIN4FHHD3uuDznwFuBMabWU6wtzkTqAr6VwKzgMrg8O84oD6pvVvyOmdx99XAaoDS0lIvKysb0BcoLy9noNvIVBqb1DQuqWlcepfuY/PKvhO0xdazouwqyhZPG7TtDve4XMw50RXAVqCZxJ7cn5vZ/Wb2DjMbm9RvHNA2CDUdBpaZWUFwbvNWYAfwIok9YYCVwLPB6zXBe4Ll64LztmuAFcHs3TkkQv7VQahPREQG6KW9deRkGTem4UXnk13M4dx/Bl4ClgL/CCwhMds1D4ib2T4SIXcVsGWgBbn7BjN7GtgEdAFvkNhL/CnwpJn9bdD2SLDKI8D3zayCxB7oimA7283sqaC2LuA+d48NtD4RERm4X+2uY+klEyjMzw27lAG5YIi6exOJvcAXzexe4CMkwnIxiWC9BlhEItj+ajCKcvcHgQd7NO8nxezaYNLTB3rZzpeBLw9GTSIiMjiqGlrZUX2KB+5cGHYpA9aniUXuvijp7ebgISIictFe2FkDwLuuKA65koFLz0tEiIhI2lq7s5Y5k8cwb8qYsEsZMIWoiIgMm6a2Tl7Zd5x3XTGVxNzR9HYxl/171syuudgNBte5/ZyZ/eHAShMRkUzz673H6Yx5RhzKhYvbEz0MrA8u7n6/mS0Nfo95hpnNMLN7zOwRElcZ+jiJ2bUiIiJn/HJHDeMLcrn20glhlzIoLmZ27qfM7Osk7t7y1yR+D+pmdgpoByYAuSSuEPRq0O/77h4fqqJFRCT9dMXirNtdyy0Lp5KTnRlnEy9qdq677wM+ZWZ/QuLi8teTuJ9oPnAC2AW85O6HhqpQERFJb68fOknD6U5+J0MO5ULfr527xN1/ReIG3CIiIhftlztryMvO4qY0vmtLT33dn37RzN45JJWIiEjGcnfW7qjhhnmTGDuqr/tv0dXXEP134Dkze1/PBcG1dF8enLJERCST7Ktr4eCJ07xrUeYcyoU+hqi7fxL4PySuYfuHAGb2NjP7MYnr62bGdCsRERlUz28/BsC7rpgaciWDq8/71O7+RTM7CnzbzD4IvJ3E/UU/DnxvkOsTEZEM8JMt1ZReOoHp40aHXcqg6vMc4+Dm25cBMeAmYD2wwN0f089aRESkp311zeysPsVdV04Pu5RB16cQNbMHSdxN5T7gn0jsfZaSuF2aiIjIOZ7bUo0Z3Lkk80K0r4dz/wJ4GPgbd68BMLPDwH+aWTHwYXfvHOQaRUQkjf10a+JQ7rRx+WGXMuj6ejj3Cnf/o+4ABXD3dcA7gd8Gfj6YxYmISHqrqG1i17Em7npb5u2FQt9n5+7rpX0T8A5g9iDUJCIiGeKnW44lDuUqRM/P3SuAGwdreyIikv5+urWK35o9keKizDuUC4N8P9Hkw7wiIjKy7alpYk9NM+/JwFm53TLjMvoiIhI5Pw1m5d6xZFrYpQwZhaiIiAw6d+fHW6q4bvZEphZm5qFcUIiKiMgQ2FLZyP66Ft57TUnYpQwphaiIiAy6H22qZFROFu/O4POhoBAVEZFB1tEVZ82bVdy2eBpF+blhlzOkIhmiZjbezJ42s11mttPMbjCziWa21sz2Bs8Tgr5mZg+ZWYWZbTGzpUnbWRn032tmK8P7RiIiI8e6XbU0nO7kd5dm9qFciGiIAt8Afu7uC4GrgJ3AA8AL7r4AeCF4D3AnsCB4rAK+DWculP8gcD1wHfBgd/CKiMjQeWZTJVMKR3HT/MlhlzLkIheiZlYE3Aw8AuDuHe7eACwHHg+6PQ7cE7xeDnzPE9YD481sOnA7sNbd6939JLAWuGMYv4qIyIhT39LBi7truefqGeRkRy5iBl0Uv+FcoA74rpm9YWYPm9kYoNjdqwGC5+47u5aQuJ9pt8qgrbd2EREZIj9+s4rOmPO+a2eGXcqw6PNNuYdBDrAU+JS7bzCzb/DWodtULEWbn6f93A2YrSJxKJji4mLKy8v7VHBPzc3NA95GptLYpKZxSU3j0ruojs1jr7RySWEWx3Zt4tiu4f/84R6XKIZoJVDp7huC90+TCNEaM5vu7tXB4drapP6zktafCVQF7WU92stTfaC7rwZWA5SWlnpZWVmqbhetvLycgW4jU2lsUtO4pKZx6V0Ux6aitokDP3+Jv7zrCspumhtKDcM9LpE7nOvux4AjZnZ50HQrsANYA3TPsF0JPBu8XgN8NJiluwxoDA73Pg/cZmYTgglFtwVtIiIyBJ589Qg5Wcbyq0fOmbMo7okCfAp4wszygP3Ax0gE/lNmdi9wGPhA0Pc54N1ABXA66Iu715vZl4DXgn5fdPf64fsKIiIjR1tnjKc3VXL74mlMKRwVdjnDJpIh6u6bgdIUi25N0deB+3rZzqPAo4NbnYiI9PSzbdU0nO7kQ9dfEnYpwypyh3NFRCT9PLH+MHMnj+GGeZPCLmVYKURFRGRAdh07xcZDJ/mf11+CWaofRmQuhaiIiAzIv284TF5OFu9bOjJ+G5pMISoiIv3W0t7FM5uO8p63TWfCmLywyxl2ClEREem3H79ZRXN7Fx9aNrImFHVTiIqISL+4O09sOMzlxYUsvWRk3t9DISoiIv3y2sGTbD3ayEduuHTETSjqphAVEZF++c6v9zOhIHdETijqphAVEZE+O3C8hV/urOEjyy5ldF522OWERiEqIiJ99ujLB8jNyuLDN1wadimhUoiKiEifNJzu4D9eP8Lyq2cwtTA/7HJCpRAVEZE+eWLDYdo64/x+SLc7ixKFqIiIXLT2rhiP/eYgNy2YzOXTCsMuJ3QKURERuWhrNldR19SuvdCAQlRERC5KVyzOv5bvY9H0Im5eMDnsciJBISoiIhflJ1uqOXC8hftvXTBiL67Qk0JUREQuKBZ3vvViBZcXF3LbouKwy4kMhaiIiFzQz7ZVU1HbzKdunU9WlvZCuylERUTkvOJx55svVDB/6ljuXDI97HIiRSEqIiLn9YsdNeyuaeKP3zmfbO2FnkUhKiIivYrHnYde2MucyWN4z5XaC+1JISoiIr368ZYqdlSf4v5b55OTrcjoSSMiIiIpdXTF+eovdnPF9CKWX1USdjmRFNkQNbNsM3vDzH4SvJ9jZhvMbK+Z/dDM8oL2UcH7imD57KRtfCFo321mt4fzTURE0tMTGw5xpL6VB+5cqBm5vYhsiAKfBnYmvf974GvuvgA4CdwbtN8LnHT3+cDXgn6Y2SJgBbAYuAP4VzMbuTe9ExHpg6a2Tr65roIb503S1YnOI5IhamYzgbuAh4P3BtwCPB10eRy4J3i9PHhPsPzWoP9y4El3b3f3A0AFcN3wfAMRkfT2nZf2U9/SwQN3LtTVic4jkiEKfB34PBAP3k8CGty9K3hfCXQfoC8BjgAEyxuD/mfaU6wjIiK9qD3Vxnd+fYD3XDmdK2eOD7ucSMsJu4CezOw9QK27v25mZd3NKbr6BZadb52en7kKWAVQXFxMeXl5X0o+R3Nz84C3kak0NqlpXFLTuPRuKMfmka3tdHTFuGlcQ9qN/3D/mYlciAJvB+42s3cD+UARiT3T8WaWE+xtzgSqgv6VwCyg0sxygHFAfVJ7t+R1zuLuq4HVAKWlpV5WVjagL1BeXs5At5GpNDapaVxS07j0bqjGZtPhk/z657/hE789l9+784pB3/5QG+4/M5E7nOvuX3D3me4+m8TEoHXu/iHgReD9QbeVwLPB6zXBe4Ll69zdg/YVwezdOcAC4NVh+hoiImknFncefHY7xUWj+NQtC8IuJy1EcU+0N38GPGlmfwu8ATwStD8CfN/MKkjsga4AcPftZvYUsAPoAu5z99jwly0ikh5++NoRth5t5BsrrmbsqHSKh/BEepTcvRwoD17vJ8XsWndvAz7Qy/pfBr48dBWKiGSGky0d/MPzu7huzkTuvmpG2OWkjcgdzhURkeH3T2t309TWxReXL9ZPWvpAISoiMsJtPFjPExsO89EbLmXhtKKwy0krClERkRGsrTPG53+0hRnjRvMnt10edjlpJ9LnREVEZGg99MJe9te18L2PX6fJRP2gPVERkRFq29FG/t9L+/nAtTO5+bIpYZeTlhSiIiIjUGcszv96eguTxuTxl3ctCructKV9dxGREeib6yrYWX2K1R+5lnEFuWGXk7a0JyoiMsK8drCeb63by+8uLeG2xdPCLietKURFREaQxtZOPvPkZmZNLOCLy5eEXU7a0+FcEZERwt358//cSs2pNp7+5I2ajTsItCcqIjJCPP16JT/dUs3nbruMq2fpPqGDQSEqIjIC7D7WxINrtnPD3El84uZ5YZeTMRSiIiIZrrG1k098fyNjRuXwjRVXk52la+MOFoWoiEgGi8edz/1wM5UnW/n2h5YytSg/7JIyikJURCSDPbRuLy/squWv/sciSmdPDLucjKMQFRHJUGt31PD1X+7lfUtn8pFll4ZdTkZSiIqIZKAtlQ3c/4M3uHLmOL783iW6R+gQUYiKiGSYI/Wn+fhjG5k4Jo+HV5aSn5sddkkZS7+0FRHJII2nO/nYY6/R0RXjyVXXM7VQE4mGkkJURCRDtHXG+MS/beTQiRYe//h1zJ9aGHZJGU8hKiKSATq64vzRE5tYv7+er//e1dw4b3LYJY0IOicqIpLmumJxPvPDN1i3q5a/vWcJ91xTEnZJI4ZCVEQkjcXjzuef3sJzW4/xl3ddwYf1U5ZhpRAVEUlTsbjzwDNbeOaNo3z2XZfx+zfNDbukESdyIWpms8zsRTPbaWbbzezTQftEM1trZnuD5wlBu5nZQ2ZWYWZbzGxp0rZWBv33mtnKsL6TiMhg64o7n37yDZ7aWMn9t8zn/lvnh13SiBS5EAW6gD9x9yuAZcB9ZrYIeAB4wd0XAC8E7wHuBBYEj1XAtyERusCDwPXAdcCD3cErIpLO2jpjfOuNdn6ypZov3LmQz912uS6mEJLIhai7V7v7puB1E7ATKAGWA48H3R4H7gleLwe+5wnrgfFmNh24HVjr7vXufhJYC9wxjF9FRGTQNbV1cu/jr7G5LsaX7lnCJ35btzULU6R/4mJms4FrgA1AsbtXQyJozWxq0K0EOJK0WmXQ1lu7iEhaqmpo5eOPvUZFbTN/8LY8XQ83AiIbomY2FvgR8Bl3P3WeQxWpFvh52lN91ioSh4IpLi6mvLy8z/Uma25uHvA2MpXGJjWNS2oal7ccbIzx9U3ttMeczy7N59L8Vo1NCsP9ZyaSIWpmuSQC9Al3fyZorjGz6cFe6HSgNmivBGYlrT4TqAray3q0l6f6PHdfDawGKC0t9bKyslTdLlp5eTkD3Uam0tikpnFJTeOS8Ivtx/iHdZuZUJDPUx/7LS4rLtTY9GK4xyVy50Qtscv5CLDT3f85adEaoHuG7Urg2aT2jwazdJcBjcFh3+eB28xsQjCh6LagTUQkLcTizj8+v4tV33+dBVPH8p/33chlxbqUX5REcU/07cBHgK1mtjlo+3PgK8BTZnYvcBj4QLDsOeDdQAVwGvgYgLvXm9mXgNeCfl909/rh+QoiIgNzsqWD+598g1/vPRiDcUgAAAyBSURBVM6K35rFX9+9WHdjiaDIhai7v0zq85kAt6bo78B9vWzrUeDRwatORGTobTxYz6ef3ExdUztf+d23seK6S8IuSXoRuRAVERmpumJxvrmugm+u20vJhNH8xx/ewFWzxoddlpyHQlREJAIOnWjhsz/czKbDDfzu0hL+5u7FFObnhl2WXIBCVEQkRLG4893/PsBXf7Gb3OwsHvrgNdx91Yywy5KLpBAVEQnJ7mNN/NmPtrD5SAO3LJzKl9+7hOnjRoddlvSBQlREZJg1tXXyzXUVPPryAYpG5/KNFVdz91UzdP3bNKQQFREZJu7Of20+yt89t4u6pnY+cO1MHrhzIZPGjgq7NOknhaiIyDB4Zd8J/v7nu9h8pIGrZo5j9Ueu5ZpLdGOpdKcQFREZQturGvmHn+/mV3vqmFaUzz+8/0rev3QmWVk6dJsJFKIiIkNg29FGvrluL89vr2Hc6Fy+cOdCVt44W1cdyjAKURGRQeLuvH7oJP/yYgUv7q6jMD+H+2+Zz703zWXcaP3mMxMpREVEBqgzFudn247xyMsHePNIA+MLcvnT2y7jozfOpkgXTMhoClERkX6qOdXGD187wpOvHqaqsY05k8fwpeWLed+1MynI01+vI4H+K4uI9EFXLM6v9x7nB68e5oVdtcTizk0LJvOle5bwzsunasLQCKMQFRG5AHdne9Upntl0lDVvHuV4cweTxuTxBzfN5YPXzeLSSWPCLlFCohAVEUnB3dlRfYqfbqnmua3VHDxxmrzsLG69YirvvaaEssunkpeTFXaZEjKFqIhIoDMW59UD9fxyZw2/3FnDkfpWsrOMG+ZOYtXN83j326YxviAv7DIlQhSiIjKiHak/zUt76/j1nuP8977jNLV1kZeTxTvmT+aPyuZz++JpTByj4JTUFKIiMqJUNbSy4cAJ1u+rZ/2BExw6cRqAGePyuett07ll4VTesWCyZtfKRdGfEhHJWB1dcXYdO8UbhxvYeOgkmw6d5GhDKwBF+TlcP3cSK2+Yzc2XTWHelDG6i4r0mUJURDJCW2eMPTVN7Kg6xbaqRrZWNrKzuomOWByAaUX5XDt7Ah9/xxyWzZ3IwmlFZOvnKDJAClERSSsdXXEO17dQUdvM7mPN7KltYs+xJvYfbyEWdwDGjsphSUkRH3v7bK6cOZ6rLxlPyXjd7FoGn0JURCKnvSvG0ZOtbK3r4vArBzl4/DSHTrRw4HgLh+pPnwlLM5g1oYDLigu5ffE0Fs8oYtGMImZNKNBFD2RYKERFZFjF487xlnZqGts5dqqNY42tHG1oo6qhlaqGVipPtlLT1IZ7sMLr2xmdm82lkwpYOL2Qu66czrwpY5k7ZQzzp47VBCAJlf70iciAtXXGOHm6g/qWtx4nmjs43tx+5rm2qZ26pnaON7fTFfez1s/NNqaNy2f6uNHcOH8SsyYUMGtiAScO7ea973o7UwpHadKPRFLGh6iZ3QF8A8gGHnb3r4RckkikdMbitLR30dIRo6W9i+b2Lprb3no+1dZJU9JzY2snja2dnGrtpOF0Jw2tHbR1xlNuOzvLmDQmj0ljRzG1cBQLpxUypXAU08blU1yUz7SifKaNy2fK2FEpD7+WN1UwtSh/qIdApN8yOkTNLBv4F+B3gErgNTNb4+47wq1MJMHd6Yo77TGnsbWTrlicrrjT0RWnM+l1RyxOZ/Dc0RU/09beGac9Fqe9M0Z7V+K5rfu5M05rZ4y2zthZz6c7YrR2JJ5Pd3TRGfMLFwoUjsqhMD+HotG5jBudy6yJBbytJJcJY/IYX5DLhII8JhTkMWlsHhMKcpk0ZhTjRufq3KRktIwOUeA6oMLd9wOY2ZPAcmDIQrS2qY03arvo3FEDJP6SPJ/kpefvevbC5L49Vzt7WV/W85T9urfT3XbmOWk9f6vjWZ/rnmhKPDu7D3dSuf5QYmmwnvvZ23hrnURL3P2c7ZxZxyGe3EZyu5+1LB5PvE9s763X8WBb3a8T/ZxY9/Lu93GCZz/z3P26K57oF3OnK/bWsuT3XfFEKMZiTmc8TlfMzz6sufYXDIb83Czyc7PJz8l+63VuNqNzsykuzCU/L5uC3GwK8rIpGJVDQW42Y0blMGZU93MOhaNyGJufw5i8RGiOHZWjn4OIpJDpIVoCHEl6Xwlc37OTma0CVgEUFxdTXl7e7w/cXNvFNza1w6aN/d5GxtuxbUg3byRmbVry6+B9lr3VlnVmmZ27LLn/mb5GVrAsK1h29sPINsjrfp8D2d3tWYm27O5HVnbSa4h1dlCQP+pMW05W4lBoTlbifW4W5ATvc7IgN8vICfrlZVuwPNHv3HOH8eDR2fugdQWPlsTb5uARtubm5gH9/5jJNDapDfe4ZHqIpvqn8zn7e+6+GlgNUFpa6mVlZf3+wKVtnYwf9RKlpaUXvU7y33mWsuRz+11ovbOXnW87565nZ173WJa0bvfnnfU5dvZ6ltRuGGbwym9+w41vv/HMNs6EXRBkZ22/O8ySlp1ZJ+l1VhBu6TzxpLy8nIH8uctUGpfeaWxSG+5xyfQQrQRmJb2fCVQN5QcW5ecye1w2S0rGDeXHpK3x+VlMLdREERHJDJl+M7zXgAVmNsfM8oAVwJqQaxIRkQyR0Xui7t5lZn8MPE/iJy6Puvv2kMsSEZEMkdEhCuDuzwHPhV2HiIhknkw/nCsiIjJkFKIiIiL9pBAVERHpJ4WoiIhIPylERURE+kkhKiIi0k8KURERkX6yC91lZKQxszrg0AA3Mxk4PgjlZCKNTWoal9Q0Lr3T2KQ2FONyqbtPSbVAIToEzGyju1/8FehHEI1NahqX1DQuvdPYpDbc46LDuSIiIv2kEBUREeknhejQWB12ARGmsUlN45KaxqV3GpvUhnVcdE5URESkn7QnKiIi0k8K0UFmZneY2W4zqzCzB8KuJwrMbJaZvWhmO81su5l9OuyaosTMss3sDTP7Sdi1RImZjTezp81sV/Bn54awa4oCM/ts8P/RNjP7gZnlh11TWMzsUTOrNbNtSW0TzWytme0NnicMZQ0K0UFkZtnAvwB3AouAD5rZonCrioQu4E/c/QpgGXCfxuUsnwZ2hl1EBH0D+Lm7LwSuQmOEmZUA9wOl7r4EyAZWhFtVqB4D7ujR9gDwgrsvAF4I3g8Zhejgug6ocPf97t4BPAksD7mm0Ll7tbtvCl43kfjLsCTcqqLBzGYCdwEPh11LlJhZEXAz8AiAu3e4e0O4VUVGDjDazHKAAqAq5HpC4+4vAfU9mpcDjwevHwfuGcoaFKKDqwQ4kvS+EoXFWcxsNnANsCHcSiLj68DngXjYhUTMXKAO+G5wqPthMxsTdlFhc/ejwFeBw0A10Ojuvwi3qsgpdvdqSPwDHpg6lB+mEB1clqJN058DZjYW+BHwGXc/FXY9YTOz9wC17v562LVEUA6wFPi2u18DtDDEh+XSQXB+bzkwB5gBjDGzD4db1cimEB1clcCspPczGcGHWpKZWS6JAH3C3Z8Ju56IeDtwt5kdJHHo/xYz+7dwS4qMSqDS3buPWDxNIlRHuncBB9y9zt07gWeAG0OuKWpqzGw6QPBcO5QfphAdXK8BC8xsjpnlkTjhvybkmkJnZkbi3NZOd//nsOuJCnf/grvPdPfZJP6srHN37VUA7n4MOGJmlwdNtwI7QiwpKg4Dy8ysIPj/6lY04aqnNcDK4PVK4Nmh/LCcodz4SOPuXWb2x8DzJGbNPeru20MuKwreDnwE2Gpmm4O2P3f350KsSaLvU8ATwT9I9wMfC7me0Ln7BjN7GthEYtb7G4zgKxeZ2Q+AMmCymVUCDwJfAZ4ys3tJ/KPjA0Nag65YJCIi0j86nCsiItJPClEREZF+UoiKiIj0k0JURESknxSiIiIi/aQQFRER6SeFqIicw8yKzOyvzeyKsGsRiTKFqIikUkrih+u5YRciEmUKURFJ5RqgHV1qT+S8dMUiETmLme0EFvZo/pG7vz+MekSiTCEqImcxs98icVeZ7cDfBc3V7n4ovKpEokkXoBeRnt4kcRu/b7r7+rCLEYkynRMVkZ4WA3kk7hQiIuehEBWRnpYCDmy+UEeRkU4hKiI9XQPsc/dTYRciEnUKURHpaRH6aYvIRdHEIhHpqQFYama3A43AXnc/EXJNIpGkn7iIyFnMbAnwCHAlkA/c5O4vh1uVSDQpREVERPpJ50RFRET6SSEqIiLSTwpRERGRflKIioiI9JNCVEREpJ8UoiIiIv2kEBUREeknhaiIiEg/KURFRET66f8DRGdwjszIQksAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def euler(f,u0,t0,t_max,dt):\n", " u=u0; t=t0\n", " # these lists will store all solution values \n", " # and associated time levels for later plotting\n", " u_all=[u0]; t_all=[t0]\n", " while tExercise 2.6: Implementing Heun's method\n", "\n", "Repeat the previous exercise for this method.\n", "\n", "For some ODEs you know the exact solution to compare the errors between Euler's and Heun's method, and how they vary with time step." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdEAAAFDCAYAAABlbTTTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU1f3/8dcnGwk7yFpAAUURUJBdRUVRxFYFq1YrKq6otbVg7Rftt1b7rQpttWrVXy0CrrhQCeK+kBAFFJStJBGURYQAsoUlAZKQzPn9cSeQZRKTkMydCe/n4zGdmXPu8pkr5cM599xzzDmHiIiIVF+M3wGIiIhEKyVRERGRGlISFRERqSElURERkRpSEhUREakhJVEREZEaivM7gEjTqlUr17lz5yM6xr59+2jUqFHtBFTP6NpUTNemcro+FdO1qVhtXJslS5bscM61DlWnJFpG586dWbx48REdIy0tjaFDh9ZOQPWMrk3FdG0qp+tTMV2bitXGtTGz7yuqU3euiIhIDSmJioiI1JCSqIiISA3pnmgVHDx4kKysLPLy8qq0fbNmzVi5cmUdR1X3EhMT6dixI/Hx8X6HIiISkZREqyArK4smTZrQuXNnzOxHt8/JyaFJkyZhiKzuOOfYuXMnWVlZdOnSxe9wREQikrpzqyAvL49jjjmmSgm0vjAzjjnmmCq3vkVEjkZKolV0NCXQYkfjbxYRqQ4l0XrowQcf5NFHH/U7DBGRek9JVCgsLPQ7BBGRqKQkGiVeeeUVBg4cSJ8+fbjtttsoKiqicePGh+rffPNNbrjhhnL7rV27lhEjRtCvXz/OOussVq1aBcANN9zA3XffzbnnnsuECRPC9TNERMJiz549PPPMM2zdurVOz6MkGgVWrlzJG2+8wYIFC1i+fDmxsbFMnz69SvuOHTuWp556iiVLlvDoo4/yq1/96lDdt99+y5w5c3jsscfqKnQREV8sX76cX//613z/fYUz9tUKPeJSTePGwfLllW9TVJREbGzVj9mnDzzxRMX1KSkpLFmyhAEDBgBw4MAB2rRp86PHzc3N5fPPP+fKK688VJafn3/o85VXXklsdQIVEYkS6enpAHX+iJ6SaBRwzjFmzBgmTpxYqrxkCzLUoyiBQIDmzZuzvIKsr1UfRKS+ysjIoHnz5rRq1apOz6MkWk2VtRiL5eQcqNXJFoYNG8bIkSMZP348bdq0ITs7m5ycHNq2bcvKlSs56aSTmDVrVrlzNm3alC5duvCf//yHK6+8EuccK1asoHfv3rUWm4hIJEpPT+eUU06p80f1dE80CvTo0YOHHnqI4cOHc+qpp3LBBRewZcsWJk2axMUXX8x5551H+/btQ+47ffp0pk6dSu/evenZsyezZ88Oc/QiIuHlnCMjI4NevXrV+bnUEo0SV111FVdddVW58iuuuKJc2YMPPnjoc5cuXfjwww/LbfPCCy/UZngiIhFj48aN7N27l1NOOaXOz6WWqIiI1CvFg4qUREVERKopIyMDgJ49e9b5uZRERUSkXklPT6djx460aNGizs+lJCoiIvVK8cjccFASFRGReuPgwYOsWrUqLCNzQUlURETqkdWrV1NQUKCWqJS2fv36sP3LSkQkWhUPKlISFRERqab09HRiY2Pp3r17WM6nJBpFioqKuPXWW+nZsyfDhw/nwIEDlS519uabbx7at3jZtLS0NIYOHcoVV1xB9+7dGT16NM45X36PiEhtS09Pp1u3biQmJoblfEqiUWT16tXceeedZGZm0rx5c2bOnFnpUmcVWbZsGU888QRff/0169atY8GCBWGIXkSk7oVrur9imvavuqqwFlpSURG1uhZaUJcuXejTpw8A/fr1Y/369ZUudVaRgQMH0rFjx+Cp+7B+/XqGDBlS9XhFRCLQvn37WLduHddff33YzqkkGkUaNGhw6HNsbCxbt26tcKmzuLg4AoEA4E3GXFBQUOFxCgsL6zBqEZHwyMzMxDkXtkFFoCRafVVoMR7IyanVpdAqUtlSZ507d2bJkiX84he/YPbs2Rw8eLDO4xER8VPxyNxwduf6dk/UzDqZ2VwzW2lmmWb222B5SzP7xMxWB99bBMvNzP5pZmvMbIWZ9S1xrDHB7Veb2ZgS5f3MLD24zz+trheW80FFS53deuutfPrppwwcOJBFixZpAW4RqffS09NJSkqia9euYTunny3RQuB3zrmlZtYEWGJmnwA3ACnOuUlmdi9wLzABuAjoFnwNAv4FDDKzlsADQH/ABY/ztnNuV3CbscBC4H1gBPBBGH9jrencufOhf2UB3HPPPYc+h1rqrG3btixcuPDQ94kTJwIwdOhQhg4deqj86aefroNoRUTCLyMjg549exJbnTEpR8i3lqhzbotzbmnwcw6wEugAjAReDG72IjAq+Hkk8JLzLASam1l74ELgE+dcdjBxfgKMCNY1dc594bxnOF4qcSwREaln0tPTwz4pTUQ84mJmnYHTgEVAW+fcFvASLdAmuFkHYGOJ3bKCZZWVZ4UoFxGRemb79u1s3bo1rIOKIAIGFplZY2AmMM45t7eS25ahKlwNykPFMBav25e2bduSlpZWqr5Zs2bk5ORUFFc5RUVF1do+kuXl5ZW7HkciNze3Vo9Xn+jaVE7Xp2K6Nt7z7wCBQKDUtajra+NrEjWzeLwEOt05lxws3mpm7Z1zW4JdstuC5VlApxK7dwQ2B8uHlilPC5Z3DLF9Oc65ycBkgP79+7uS9wwBVq5cWa3RtjlhGp0bDomJiZx22mm1drziGZOkPF2byun6VEzXBlasWAHA6NGjad++/aHyur42fo7ONWAqsNI5948SVW8DxSNsxwCzS5RfHxylOxjYE+zu/QgYbmYtgiN5hwMfBetyzGxw8FzXlziWiIjUI+np6RxzzDG0a9curOf1syV6JnAdkG5mxbMF/AGYBMwws5uBDUDxdDzvAz8F1gD7gRsBnHPZZvYX4Kvgdv/nnMsOfr4DeAFIwhuVG5Ujc0VEpHIZGRmccsophPtJRt+SqHNuPqHvWwIMC7G9A+6s4FjTgGkhyhcDWj+sEmlpaSQkJHDGGWf4HYqISI0EAgEyMjK44YYbwn7uiBidK/5JS0vj888/9zsMEZEa27BhA7m5uWEfmQtKolHjlVdeYeDAgfTp04fbbruN77//nm7durFjxw4CgQBnnXUWH3/8MQCjRo2iX79+9OzZk8mTJx86xocffkjfvn3p3bs3w4YNY/369Tz77LM8/vjj9OnTh3nz5vn180REaiw9PR0I73R/xXx/xEV+3MqVK3njjTdYsGAB8fHx/OpXv+LTTz9lwoQJ3H777QwaNIgePXowfPhwAKZNm0bLli05cOAAAwYM4PLLLycQCHDrrbfy2Wef0aVLF7Kzs2nZsiW33347jRs3LjUDkohINFESjSLjxo0LuWpKSUVFRdWadqpPnz48UcnE9ikpKSxZsoQBAwYAcODAAdq0acODDz7If/7zH5599tlSMf3zn/9k1qxZAGzcuJHVq1ezfft2zj77bLp06QJAy5YtqxyfiEgky8jI4LjjjqNp06ZhP7eSaBRwzjFmzJhD898W279/P1lZ3qRMubm5NGnShLS0NObMmcMXX3xBw4YNGTp0KHl5eTjnwj5qTUQkHPyY7q+Ykmg1VdZiLFbbky0MGzaMkSNHMn78eNq0aUN2djY5OTk8+uijjB49muOOO45bb72Vd999lz179tCiRQsaNmzIqlWrDk1Cf/rpp3PnnXfy3XfflerObdKkCXv37q21WEVEwqmgoIBVq1Zx8cUX+3J+DSyKAj169OChhx5i+PDhnHrqqVxwwQWsX7+er776igkTJjB69GgSEhJ4/vnnGTFiBIWFhZx66qncf//9DB48GIDWrVszefJkfv7zn9O7d2+uuuoqAC655BJmzZqlgUUiEpW+/fZbCgsLfRmZC2qJRo2rrrrqUOIrVnKps+Tk5EOfP/gg9JwSF110ERdddFGpshNPPPHQdFkiItHGz0FFoJaoiIhEsYyMDOLi4ujevbsv51cSFRGRqJWens5JJ51EQkKCL+dXEhURkajl58hcUBKtMm/q3qPL0fibRSR65OTksH79et8GFYGSaJUkJiayc+fOoyqpOOfYuXMniYmJfociIhJSZmYmgK9JVKNzq6Bjx45kZWWxffv2Km2fl5dXL5JPYmIiHTt2/PENRUR84PfIXFASrZL4+PhD0+VVRVpaGqeddlodRiQiIhkZGTRq1IjOnTv7FoO6c0VEJCqlp6fTs2dPYmL8S2VKoiIiEnWcc6Snp1d8P/S77+Ckk2i+dGmdxqEkKiIiUWfbtm3s2LGj4iS6aBF8+y2FtTiPeShKoiIiEnV+dFDRwoWQlMS+rl3rNA4lURERiToZGRlAJY+3LFoE/fvjqrG2c00oiYqISNRJT0+nTZs2tGnTpnxlfj4sXQqDBtV5HEqiIiISdSqd7u+//4WCAgguBVmXlERFRCSqBAIBMjMzK+/KBbVERUREyvruu+/Yv39/xUl04ULo0AHCMOOakqiIiESVHx2Zu2hRWFqhoCQqIiJRpnhkbs+ePctX7tgBa9cqiYqIiISSnp5O165dady4cfnK4vuhYRhUBEqiIiISZTIyMirvyo2NhX79whKLkqiIiESN/Px8vvnmm8pH5vbqBY0ahSUeJVEREYkaq1atoqioKHRLNBDwkmiYunJBSVRERKJIpdP9ffst7NkTtkFFoCQqIiJRJD09nfj4eE488cTylQsXeu9KoiIiIuWlp6fTvXt34uPjy1cuWgRNm0L37mGLR0lURESiRkZGRuWDigYOhJjwpTYlURERiQp79uxhw4YNoZPo/v2wYkVYBxWBkqiIiESJzMxMoILp/pYsgaKisN4PBSVRERGJEsVz5oZsiYZx5ZaSlERFRCQqpKen06RJE4499tjylQsXQteu0Lp1WGNSEhURkahQPN2fmZWvDOPKLSUpiYqISMRzzpGenh66K3fTJsjKUhIVEREJZcuWLWRnZ1d+PzTMI3NBSVRERKJA8XR/IUfmLloECQnQp0+Yo1ISFRGRKFA8MjdkEl240EugDRqEOSolURERiQIZGRm0a9eOVq1ala4oLITFi33pygUlURERiQIVDirKzPRmK/JhUBEoiYqISIQrKioiMzMzdBItXrlFLVEREZHy1q1bR15eXsWDilq1gi5dwh8YPiZRM5tmZtvMLKNE2YNmtsnMlgdfPy1Rd5+ZrTGzb8zswhLlI4Jla8zs3hLlXcxskZmtNrM3zCwhfL9ORERqy49O9zdoEISagCEM/GyJvgCMCFH+uHOuT/D1PoCZ9QCuBnoG9/l/ZhZrZrHAM8BFQA/gl8FtAf4aPFY3YBdwc53+GhERqRPp6emYGT169ChdsWcPrFzpW1cu+JhEnXOfAdlV3Hwk8LpzLt859x2wBhgYfK1xzq1zzhUArwMjzZsT6jzgzeD+LwKjavUHiIhIWGRkZHD88cfTsGHD0hVffQXO+TaoCCLznuivzWxFsLu3RbCsA7CxxDZZwbKKyo8BdjvnCsuUi4hIlKlwZG7xTEUDBoQ3oBLifDtzaP8C/gK44PtjwE1AqM5uR+h/BLhKtg/JzMYCYwHatm1LWlpatYIuKzc394iPUV/p2lRM16Zyuj4Vq8/XpqCggNWrVzNo0KByv7HXe++RdNxxfLV8eYX71/W1iagk6pzbWvzZzJ4D3g1+zQI6ldi0I7A5+DlU+Q6guZnFBVujJbcPdd7JwGSA/v37u6FDhx7R70hLS+NIj1Ff6dpUTNemcro+FavP12bZsmUEAgEuvvji0r/ROVizBn72s0p/e11fm4jqzjWz9iW+XgYUj9x9G7jazBqYWRegG/Al8BXQLTgSNwFv8NHbzjkHzAWuCO4/Bpgdjt8gIiK1Z8WKFUCIkbnffQfbt/t6PxR8bIma2WvAUKCVmWUBDwBDzawPXtfreuA2AOdcppnNAL4GCoE7nXNFweP8GvgIiAWmOecyg6eYALxuZg8By4CpYfppIiJSSz777DNatGhBt27dSlf4uHJLSb4lUefcL0MUV5jonHMPAw+HKH8feD9E+Tq80bsiIhKlUlNTOffcc4mNjS1dsWgRNGwIoSZgCKOI6s4VEREptm7dOtavX895551XvnLhQujfH+L8HdqjJCoiIhEpJSUFgGHDhpWuyM+HZct8vx8KSqIiIhKhUlNTad++PSeddFLpiv/+FwoKlERFRERCcc6RmprKsGHDsLLz4vq8cktJSqIiIhJxMjMz2bZtW/muXPAGFXXo4L18piQqIiIRp/h+aMhBRcUrt0QAJVEREYk4KSkpnHDCCRx77LGlK7Zvh7VrI6IrF5RERUQkwhQWFvLpp5+GboV++aX3rpaoiIhIeUuWLGHv3r2h74cuXAixsdCvX/gDC0FJVEREIkpqaioA5557bvnKRYvglFOgUaMwRxWakqiIiESUlJQUTj31VFq3bl26IhDwunMjpCsXlERFRCSC5OXlsWDBgtBdud98A3v2RMygIlASFRGRCPLFF1+Ql5dX8aMtoJaoiIhIKCkpKcTGxnL22WeXr1y0CJo1g7LTAPpISVRERCJGamoqAwYMoGnTpuUrFy6EgQMhJnJSV+REIiIiR7W9e/fy5Zdfhr4fum8fpKdHVFcuKImKiEiE+OyzzygqKgqdRJcsgaIiJVEREZFQUlNTSUxM5PTTTy9fGYGDikBJVEREIkRKSgpnnnkmiYmJ5SsXLYKuXaHss6M+UxIVERHfbd++nRUrVoR+tAW8QUUR1goFJVEREYkAc+fOBQh9P3TTJu8VQZMsFFMSFRER36WkpNC0aVP6hZpYPkLvh4KSqIiIRIDU1FTOOecc4uLiylcuXAgJCdCnT/gD+xFKoiIi4qsNGzawZs2a0F254LVETzsNGjQIb2BVoCQqIiK+Kl76LOSgosJCWLw4IrtyQUlURER8lpKSQuvWrenVq1f5yowM2L8/IgcVgZKoiIj4yDlHSkoK5513HmZWfoP58713tURFRERK++abb9iyZUvF90Nnz/ZWbenSJbyBVVGIYVAVM7PBwAhgMPATIAnYAXwDfAq85ZzbVdtBiohI/ZSSkgJUcD80OxvmzoX/+R8I1UqNAFVqiZrZGDNLBz4HxgENgdXAImAXMAiYAmwysxfMLDL/ySAiIhElNTWV4447jq5du5avfOcdb9L5n/88/IFV0Y+2RM3sv0Ab4CXgemC5c86F2K4ZcDEwGsg0sxudc2/UcrwiIlJPFBUVMXfuXC677LLQ90OTk6FTJwg1AUOEqEp37vPAs865vMo2cs7tAaYD082sN9CuFuITEZF6avny5ezatSt0V25uLnz0Edx+e8R25UIVunOdc08UJ1Azm2Fmnaqwz3+dcx/VRoAiIlI/Vfp86AcfQH5+RHflQvVH514BtA9VYWYtzeyMIw9JRESOBikpKZx88sm0bx8irSQne8uenXlm+AOrhh9NomZ2kpn1NLMf27YbMK92whIRkfqsoKCAefPmhX60JT8f3nsPRo2C2NjwB1cNVWmJXg2kA7mAA/5gZneZ2RAza1xiu2ZApfdNRUREABYtWsT+/ftDJ9GUFMjJifiuXKjawKJ/AJ8BfYG/A73wnhVNAAJmthb4GugNrKijOEVEpB5JSUkhJiaGc845p3xlcjI0bQoVLdAdQaoysCjHOTfXOfcYsAq4CmiCl1THAh8DTYFlwC11GKuIiNQTqamp9O3blxYtWpSuKCz0Zim6+GJv+bMIV60Zi5xzPUp8XR58iYiIVNm+fftYuHAh48ePL185fz7s2BEVXbmguXNFRCTM5s+fz8GDB0PfD01OhsREGDEi/IHVQFVG5842s9OqekAzSzSzu83s9iMLTURE6qOUlBTi4+M5s+zjK4GAl0RHjIBGjfwJrpqq0hLdACw0s0XBUbl9zaxUN7CZ/cTMRpnZVGALcBOwtA7iFRGRKJeamsrpp59Oo7KJcvFi2LQparpyoWoDi34D9AC+BB4EvgLyzCzbzLaYWR6wEUgGeuJNUH+qc+7LOotaRESiUnZ2NkuXLq24KzcuzhtUFCWqNLDIObcW+I2Z/Q44HW/Vlp8AicBOvFG7nznnvq+rQEVEJPqlpaXhnCs/1Z9zMHOm91hL2RG7Eay6o3ML8NYN/bRuwhERkfosNTWVRo0aMXDgwNIVmZmwZg3cc48/gdVQtUbnmtllWitURERqKiUlhbPOOouEss+AJid7q7WMHOlPYDVUrZYoMBNwZpYD/JfDz4ouBzKccwdrOT4REaknNm/ezKpVq7j55pvLVyYne5PNt4uuVTSr+5zoscBI4DG8e6GXAlOBxUCumVV58gUzm2Zm28wso0RZSzP7xMxWB99bBMvNzP5pZmvMbIWZ9S2xz5jg9qvNbEyJ8n5mlh7c558WcsVXEREJl+Klz8oNKlq7Fv7736galVusWknUOZflnHvXOfcX59zPnXNd8FZveQ7IwRvBW1Uv4M3BW9K9QIpzrhuQEvwOcFHwPN3wphr8F3hJF3gAb6DTQOCB4sQb3GZsif2i48ldEZF6KiUlhZYtW9K7d+/SFbNmee+XXRb+oI7QEc9Y5Jxb65y7HXiRakwD6Jz7DMguUzwyeByC76NKlL/kPAuB5mbWHrgQ+MQ5l+2c2wV8AowI1jV1zn3hnHPASyWOJSIiYeacIzU1lXPPPZeYmDKpJzkZ+vaFzp19ie1I1Oa0f08DRzpLUVvn3BaA4HubYHkHvGdRi2UFyyorzwpRLiIiPli7di0bNmwo/2jLli3wxRdR2QqFag4sMrPxeKu1LHfO7S5TnQecUFuBlT11iDJXg/LQBzcbi9f1S9u2bUlLS6tBiIfl5uYe8THqK12biunaVE7Xp2LRcG3eeecdABo3blwq1p/Mns2JwJcdO7K/Dn5DXV+b6o7OnVS8j5l9TzCh4k319wtg/RHGs9XM2jvntgS7ZLcFy7OATiW26whsDpYPLVOeFizvGGL7kJxzk4HJAP3793dDhw6taNMqSUtL40iPUV/p2lRM16Zyuj4Vi4Zr89hjj9GxY0euu+46So3zfPhhOOkkBo4Z4z3iUsvq+tpUtzu3MdAfuBV4D2gH/B4vAfUEfn2E8bwNFI+wHQPMLlF+fXCU7mBgT7C79yNguJm1CA4oGg58FKzLMbPBwVG515c4loiIhNGmTZt4//33yyfQ7GyYO9cblRulD1BUd8aig3itz2Uly82siXMupzrHMrPX8FqRrcwsC2+U7SRghpndjDfx/ZXBzd8HfgqsAfYDNwbjyTazv+DN5wvwf8654sFKd+CNAE4CPgi+REQkzF588UUCgQA33XRT6Yp33oGioqh8tKVYdbtzQ6puAg3u88sKqsrNShwcYXtnBceZBkwLUb4Y6FXduEREpPYEAgGmTZvGOeecwwknlBk2k5wMnTpBv37+BFcLtCi3iIjUmc8++4y1a9dyyy23lK7IzYWPPorqrlxQEhURkTo0ZcoUmjVrxuWXX1664oMPID8/qrtyQUlURETqyO7du5k5cybXXHMNSUlJpSuTk6F1a2++3CimJCoiInXi1VdfJS8vr3xXbl4evPsujBoFsbH+BFdLlERFRKROTJkyhT59+tC3b9/SFSkp3j3RKO/KBSVRERGpA8uWLWPZsmUVL3vWtCmUnQIwCimJiohIrZs6dSoNGjRg9OjRpSsKC2H2bLjkEii7MHcUUhIVEZFadeDAAV555RUuv/xyWrRoUbpy3jzYuTNqJ5wvS0lURERqVXJyMnv27Km4KzcxEUbUjyWelURFRKRWTZ06la5du5af+D0Q8BbgHjECGjXyJbbapiQqIiK1Zu3atcydO5ebbrqp/OLbX30FmzbVi1G5xZRERUSk1kybNo2YmBhuuOGG8pWzZkFcHFx8cdjjqitKoiIiUisKCwt54YUXuOiii+jQoUPpSudg5kzvsZayg42imJKoiIjUio8++ojNmzeHHlCUmQlr1tSrrlxQEhURkVoyZcoU2rRpw8WhumvffNNbrWXkyPAHVoeUREVE5Iht3bqVd999lzFjxhAfH1+6Mi8P/v1vOP98aNfOnwDrSK0syi0iIke3l156icLCQm666abylS+8AD/8AK++Gva46ppaoiIickScc0yZMoUhQ4bQvXv30pWFhfC3v8GgQVD2udF6QElURESOyIIFC/j2229DDyh64w347ju47z7vnmg9oyQqIiJHZOrUqTRp0oQrr7yydEUgAJMmQY8e3oTz9ZDuiYqISI3t3buXGTNmcO2119Ko7FR+770HGRnw8stQdvaieqJ+/ioREQmL119/nf3795fvynUOHnkEOneGq6/2JbZwUEtURERqbOrUqfTq1YsBAwaUrvj0U1i4EJ55xpvqr55SS1RERGokPT2dL7/8kltuuQUrO2ho4kRo0wZuvNGf4MJESVRERGpk6tSpJCQkcO2115auWLIEPv4Y7r4bkpL8CS5MlERFRKTa8vPzefnllxk1ahTHHHNM6cqJE6FZM7jjDn+CCyMlURERqbbZs2eTnZ3NLbfcUrpi1SpIToY774SmTf0JLoyUREVEpNqmTJnCcccdx7Bhw0pX/O1v0KAB/Pa3/gQWZkqiIiJSLd9//z1z5szhxhtvJKbk858bN3rPhN56qzeo6CigJCoiItXy/PPPA3Bj2ZG3jz3mvd9zT5gj8o+SqIiIVFlRURHTpk1j+PDhHHvssYcrtm+HyZNh9GgoWV7PKYmKiEiVzZkzh40bN5afoeif//TWDZ0wwZ/AfKIkKiIiVfbcc8/RqlUrLr300sOFe/fC00/DZZfBySf7F5wPlERFRKRKvvzyS2bOnMltt91GgwYNDlf8+9+we7e33NlRRklURER+VCAQ4K677qJ9+/ZMKNllm5cH//gHnH8+9O/vX4A+qb+zAouISK2ZPn06ixYt4sUXX6RJkyaHK154AX74AaZP9y02P6klKiIilcrNzWXChAkMHDiw9Dy5hYXe5AqDBsG55/oXoI/UEhURkUpNnDiRLVu2kJycXHpyhRkz4Lvv4PHHoewqLkcJtURFRKRC69at47HHHuO6665j8ODBhysCAY7dkyYAABr6SURBVG+i+R494JJL/AvQZ2qJiohIhe655x7i4uKYNGlS6Yr33oOMDHjpJYg5ettjSqIiIhJSSkoKs2bN4pFHHuEnP/nJ4QrnvFZo585w9dW+xRcJlERFRKScwsJCxo0bR5cuXRg/fnzpys8+gy++gGeegfh4fwKMEEqiIiJSzuTJk8nIyCA5OZnExMTSlY884q3SUnYC+qPQ0duRLSIiIWVnZ3P//fdz3nnnMWrUqNKVS5bAxx/D+PGQlORPgBFESVREREp54IEH2L17N0888QRW9tGVhx6CZs3gjjv8CS7CKImKiMghGRkZ/Otf/+KOO+7glFNOKV352mvw1lveSi3NmvkTYIRREhUREQCcc4wbN46mTZvy5z//uXTlxo1e6/P00+H3v/cnwAikgUUiIgLA22+/TUpKCk899RTHHHPM4YpAAMaMgaIiePlliFPqKBaRLVEzW29m6Wa23MwWB8tamtknZrY6+N4iWG5m9k8zW2NmK8ysb4njjAluv9rMxvj1e0REIl1eXh533303PXv25Pbbby9d+fjjMHcuPPkkHH+8PwFGqIhMokHnOuf6OOeK19a5F0hxznUDUoLfAS4CugVfY4F/gZd0gQeAQcBA4IHixCsiIqU98cQTrFu3jieeeIK4ki3NFSvgD3+AUaP0SEsIkZxEyxoJvBj8/CIwqkT5S86zEGhuZu2BC4FPnHPZzrldwCfAiHAHLSIS6TZv3sxDDz3EqFGjOP/88w9X5OXB6NHQsiU899xRO8l8ZSK1Y9sBH5uZA/7tnJsMtHXObQFwzm0xszbBbTsAG0vsmxUsq6i8HDMbi9eKpW3btqSlpR1R8Lm5uUd8jPpK16ZiujaV0/Wp2JFem0mTJlFQUMAVV1xR6jjH/7//R6eMDFZMmkR2RsaRB+qDOv9z45yLuBfwk+B7G+C/wNnA7jLb7Aq+vwcMKVGeAvQDfg/8sUT5/cDvfuzc/fr1c0dq7ty5R3yM+krXpmK6NpXT9anYkVybRYsWOcDde++9pSvmzHEOnLvzziMLzme18ecGWOwqyBkR2Z3rnNscfN8GzMK7p7k12E1L8H1bcPMsoFOJ3TsCmyspFxERIBAIcNddd9GuXTv+8Ic/HK7Ytcsbjdu9u7fotlQo4pKomTUysybFn4HhQAbwNlA8wnYMMDv4+W3g+uAo3cHAHud1+34EDDezFsEBRcODZSIiAkyfPp1FixYxadIkmjRp4hU65z0PunUrvPIKNGzob5ARLhLvibYFZgWnmooDXnXOfWhmXwEzzOxmYANwZXD794GfAmuA/cCNAM65bDP7C/BVcLv/c85lh+9niIhErtzcXCZMmMDAgQO57rrrDle8+iq88QY8/DD06+dfgLUgEKj7c0RcEnXOrQN6hyjfCQwLUe6AOys41jRgWm3HKCIS7SZOnMiWLVtITk4mpnhR7e+/hzvvhDPP9Kb2i2KrV8Nll8FvftOEoUPr7jwRl0RFRKRuzZs3j0cffZRrr72WwYMHe4VFRd590EDAm5UoNtbfII/A2rVw7rlQUAAJCUV1ei4lURGRo8i3337LqFGj6NKlC08++eThin/8Az79FF54Abp08S2+I/Xdd14CzcvzJlnauXN/nZ4v4gYWiYhI3dixYwc/+9nPiImJ4b333qNly5ZexfLl8L//C5dfDtdf72+QR2D9ei+B7tsHKSlQdhGauqCWqIjIUSAvL4+RI0eyceNG5s6dy/HFc+AeOADXXgutWsG//x21sxJt2ADnnQd79ngJtHe5kTV1Q0lURKSeCwQC3HDDDXz++efMmDGD008//XDlffdBZiZ8+CGUXLklimRleQk0OxvmzIG+fX98n9qiJCoiUs/98Y9/5I033uCvf/0rV1555eGKTz7xVmb5zW/gwgv9C/AIbN7sJdBt27yf07//j+9Tm5RERUTqsSlTpjBx4kTGjh3L70supp2dDTfcACefDH/9q2/xHYkffvAS6JYt8NFHMGhQ+GNQEhURqac++eQTbr/9di688EKeeeYZrPh+ZyAAY8fC9u3w7ruQlORvoDWwbZuXQLOyvJ7oM87wJw4lURGReigjI4MrrriCHj16MGPGjMNrhBYVwW23wcyZ3ry4p53mb6A1sH07DBvmjcb94AMYMsS/WJRERUTqmc2bN/PTn/6Uxo0b895779G0aVOvorDQ68KdPh3+9Ce45x5f46yJnTvh/PNhzRp47z045xx/41ESFRGpR3Jzc7nkkkvIzs5m3rx5dOoUXMyqoACuucZrgT7yiDcqN8pkZ3sJ9Ntv4Z13vO5cvymJiojUE0VFRVxzzTUsX76cd955h9OKu2rz8uCKK7ym2+OPw7hx/gZaA7t2wfDh8PXX8PbbXjKNBEqiIiL1xPjx43nnnXd45pln+OlPf+oV7tsHo0Z5D1A++6x3PzTK7NnjPYGzYgW89VZkPY2jJCoiUg88+eSTPPXUU9x999386le/8gpzcuBnP4MFC7w5cceMqfQYkWjPHhgxwpuZcOZMKP63QaTQ3LkiIlFu/vz5jB8/nssuu4y///3vXuHu3XDBBfD5594aoVGYQBct8gYPL17sLXF6ySV+R1SekqiISBRbvHgxDz/8MAMGDOCVV17x1gbdscMbdbNsmdd8u+oqv8OslkAAJk3yHl0JBLzFZS67zO+oQlN3rohIlFqxYgWXXHIJzZs35+2336Zhw4beND4XXOA9AzJ7ttcXGkU2b4brroPUVC/3P/ssNG/ud1QVU0tURCQKvfTSSwwePBgzY+LEibRt29abvuecc7xFNd97L+oS6LvvwqmnwsKFMHUqvPZaZCdQUBIVEYkq+fn53HHHHYwZM4aBAweydOlSOnfu7E3fc/bZhyeSjYSHKKsoLw/uusu759mpEyxdCjfdFB2rsimJiohEiQ0bNnD22Wfz7LPP8vvf/545c+bQrl07krKy4KyzvMFEKSlw5pl+h1plK1d6E8c/9ZT3+OrChXDSSX5HVXW6JyoiEgU+/vhjrrnmGgoKCpg5cyY///nPvYqvv6bPb38LsbEwd274VqM+Qs7BlCnw299Co0ZeV+7PfuZ3VNWnlqiISAQLBAI89NBDjBgxgvbt27N48eLDCbR48lgzSEuLmgS6axf84hfeQjJnnulNohCNCRSUREVEItauXbu49NJLuf/++/nlL3/JwoULOfHEE71HWK69Fi6+GNq1Y/mTT0KPHn6HWyULFkCfPt7MQ3/9q3f7tn17v6OqOSVREZEItGzZMvr168fHH3/M008/zSuvvEKjhg29WQd69IAZM+DBB2HJEg506OB3uD+qqAj+7/+8sU/x8V4y/Z//gZgoz0JRHr6ISP3z/PPPc8YZZ1BQUMCnn37KnXfeiW3Z4s04cPXV0LmzN4T1gQcgIcHvcCvlnHe/s18/L9xrrvFCHzjQ78hqh5KoiEiEyMvLY+zYsdx0002cccYZLF26lNMHD/YemuzRw+v7fPRR+OIL6NXL73B/VEoKnHGG9+hKbq7XiH75ZShe3rQ+UBIVEYkA69evZ8iQITz33HPcd999fPzxx7TZt8+bfeiWW7wbienp8LvfeSNxI9jnn3uPqZ5/vjf/w+TJ3qMsv/iF35HVPiVREREfFRYWMm3aNPr168eaNWt46623eOQvfyH26ae91uaXX3pz36Wmwgkn+B1upZYu9VZZOfNMyMyEJ5+E1avh1lu9+6D1kZKoiIgPAoEAb7zxBj179uTmm2/mhBNOYPHixYw88URv4oRx42DoUC8b3XZbRI/Aycz01vzu18+bLGHiRFi3zpuFKDHR7+jqVuT+VxERqYecc7zzzjv07duXq6++moSEBGbNmsXCefM44Y03vG7bb77xbh6++643D16EWrPGmyz+lFO827V/+pM3be+993oTKBwNlERFRMIkNTWVM844g0svvZTc3FymT5/O8mXLGJWUhA0aBH/8I4wa5d1AvPbaiJ08duNGb6KE7t29ldbuucdLnn/+MzRr5nd04aUkKiJSx7744guGDRvGsGHDyMrKYvLkyaxcupRrcnOJ7d3bW21l2zaYNcsbwtqmjd8hl+Oct0j27bd7t2ZfeAHuuAPWroW//Q1atfI7Qn9o7lwRkTqyfPly7r//ft59913atGnDE088wW2XXkri1KnQtSvs3AmnnQYvveQNXW3QwO+Qy9m4EV55BV580etlTkz0unDvvx+OO87v6PynJCoiUsu++eYb/vSnPzFjxgyaN2/Oww8/zF1nn03jf/8bfv97KCyEkSO9wUNnnx1x3bb79kFyspc4U1O9VuhZZ3mhX3HF0ddlWxklURGRWrJq1Sr+9re/8eKLL5KUlMT/3ncf9/ToQfPnnoP//V9o3NjrA73rLjj+eL/DLSUQgE8/9RrFb77pTY7QpYs3WOj6672Gs5SnJCoicgQ2btzI66+/zmuvvcayZcto0KABd91+O/e1a0eb55/3Rtwcdxw89hjcfHPENeNWr/YS58svw/ffQ5MmcNVVMGaM97xnBD9ZExGUREVEqmnHjh385z//4bXXXmPevHkADBgwgMcnTODqXbto9/LLkJPjZaG//93ruo2LjL9unfOe4fzoI+9e5xdfeInyggvgkUe8wcENG/odZfSIjP+qIiIRLicnh7feeovXXnuNTz75hMLCQk7u3p2/jBnD1YEAJ8yf763tFRfnDRIaNw4GDPA7bMAb+Jua6s1lO2cOrF/vlffo4YU8ejREwUIwEUlJVESkAvn5+XzwwQe8+uqrvPPOO+Tl5XFshw787oIL+GV+PqcuXIitWuUNWT3/fJgwAS691PcFMnNzYd48L2HOmeMteg3QvDmce643QGjYMDjxxIgb0xR1lERFRErYtWsXCxYsIDk5meTkZPbs2UPrFi24+dRT+eXevZy+ahUxmzZBx47eiJuLL/Yyk499oIWFxuefH06aCxfCwYPeEzNDhnjT8A0bBn37Rvzc9VFHSVREjlrOOTZs2MD8+fMPvTIzM3HO0SQpiZ936sQvY2MZlp1N3FdfeYtg/uUv3tpep57qSzPOOdiwAZYtg+XLYfFiSE09kwMHvHD69fMWejn/fG8ZsqSksId4VFESFZGjRlFREenp6cyfP58FCxYwf/58srKyAGiamMgZxxzD1a1bc+b27Qw+cIDEzZvhwgu91uZFF0HbtmGN9+BBbwbA5csPJ83ly2H3bq8+JgZOOgkuuGArY8Z0YOhQaNkyrCEe9ZRERaTe2r9/P19++eWhVubnn39OTk4OAB0bNmRIfDxDzBjiHL3y8og9eBD69/eac2ed5U2EEKZZhPbu9e5dlkyWGRlQUODVJyV5jd+rr/bmqO/Tx5v4vWFDSEtbzdChGhnkByVREYlqzjl27NjBqlWrDr1WrlzJqsxM1m/ciHMOA3olJnJtfj5DgCHAsY0aHU6Yxe8dOtRZF61zsH27N9fsunXl3zdvPrxt69bebIDjxnnJ8rTToFs33c+MREqiIhIVioqKWLNmjZcgV61i1ddfs2rFClatXk12sHUJkGTGSWYMCgQYAwwATm/ZkhYDBpROmh071nrCzMvz5poNlSTXrfNGzZbUsaM3E9CFF3qTuvfu7SXM9u01ajZaKImKiO8CgQDbt29n06ZNh1/r17Np7Vo2b9zIxk2bWL15MwcDgUP7tAW6A1cG37vHx3PyscfS6cQTiTnhBG9aveOP9/o8jz22xlnJOe8e5JYt8MMP3ntFn4vvVRZLTPSSZNeu3gDerl29kLp29abUq+8LVh8NlERFpE4EAgH27t1Ldnb2odf2TZvY9O23bFq3jk0bN7Lphx/YnJ3NlpycUgkSvHUa2wIdgBOBi4HujRpx8rHHctLJJ9Oie/fDibJrV68r9kfmqCsqgj17IDsbdu0q/V62bOvWw8kxP7/8sRITvRZj+/Zw8slw3nne5w4dDofUvr2mzavv6n0SNbMRwJNALDDFOTfJ55BEosLBgwfJzc09/Nq9m9wdO8jduZM927aR/cMPZG/bRvbOnV6S3LOH7JwcsvftI/vAAXYVFBCo4NiN8ZJjB+AcoENiIh2aNaND69Z0aN+eDscdR9vOnYnr0IFAm3bkt2jHnDVbOfn0C8nJgYxcb1a93FzIWQO5y0t8zyn9effuw0lyzx6vZVmRRo2gRQvv1aaN94xlu3aHk2XJz02bqstV6nkSNbNY4BngAiAL+MrM3nbOfe1vZCI/zjlHUVERBw8epCA/n4N5eRTs28fBAwcoOHCAgwcOeGX795O/bx95ubnk7dvHgZwc8vbvL/3KyyPvwAHy8vO9z/n5HMjPZ39eHrl5eeTm55NbUEDuwYPkFhaSW1REfmXZpoTmQMvgqwXGsXEJNI1PomnTY2iU0JiGDZqR1KAF8YmtiGvQhphGJ5DboAvbY9ryA+34IdCG7/ITOHAA77UGDqRz6Hvx6NSqaNDAWyilSZPD761be4+BtGzpJcdQ78WfExJq8l9Kjmb1OokCA4E1zrl1AGb2OjASqLMkmpG6gvf+PJ2vWn9c7X1dJX9pBQIV/Zuecv+0LnmcQInPhiv3r3Dngsd13v84d+gLDjDncKFiOxSPt4/DQcl357zzBcMrPs/u3bv4uNlLh7aBQLDegXM43KHPJY/tAoHg98P1zpX8HvA+l9g/4AJltgscKgu4AC4QIMDhbZ1zh8pc8HsgEDi0n3MBiorLD5UVfz58/IBzFAW3DQTfiyhZ7igi+Bl36HuhcxQCB3EcxFGN3FElcUBiiVcDIIkYkoglkTjak0ACjYkjkTgSgzUNgUZAEwI0wdGUQppzkJYcoCP76MAejmEPzcikOQdIgkKDQuBAmfPHeUkuKSn0q2nTiuuSkmDz5m/p1+/EQ8mxZKJs3Nh7KQlKuNX3JNoB2FjiexYwqOxGZjYWGAvQtm1b0tLSanzCDx9+nUc/m1Lj/aXuGd79tpgf+WxALBb8bsSWeS+utxLbxJR5xRJDDDHEHfrk1cRYyc9eHRZDrMURY3GYxRFLHBYTj1k8RhwxMfFgCZjFExOTADHBd0vAYpMgLgmLbQhxjSCuIcQ38l4JjXHxjSGhIUVx8RTGJeDi4yE2lri4ALGxjthYR1xc6ffY2EC5spJ1CQmO+PgAcXEHiI/fT3z8JuLji8sChz7Hxzvi4gJHfG8wNzeXxo0PPwdS3FLdvv3Ijlsf5ObmHtHfW/VZXV+b+p5EQ92xKNfcc85NBiYD9O/f3w0dOrTGJ+zaqiutnuxE1xquYBtTyd80VskNGIstvZ9hh359qf2s+K//Eucs3teCexq44DZmFtyn+HOIWC3Guzdkhpkd+hwTY8HjBY8VE8PXX39Nz149sRhvHzMvecTEGBYT3DYm+N1iDh8zxoiJifH2iw2+m/ceE2vExHrfY2JjvP3j4g59j4n1vsfGxRITF+udJ8aC5z90Wcq9vN9Yuqzs99q8J5aWlsaR/Nmr73R9KqZrU7G6vjb1PYlmAZ1KfO8IbK5g21pxbK9j6T/6TP2BrkBRWhFDhg7xOwwRkVpR3wdffwV0M7MuZpYAXA287XNMIiJST9TrlqhzrtDMfg18hPeIyzTnXKbPYYmISD1Rr5MogHPufeB9v+MQEZH6p75354qIiNQZJVEREZEaUhIVERGpISVRERGRGlISFRERqSElURERkRpSEhUREakhq2zlkKORmW0Hvj/Cw7QCdtRCOPWRrk3FdG0qp+tTMV2bitXGtTnOOdc6VIWSaB0ws8XOuf5+xxGJdG0qpmtTOV2fiunaVKyur426c0VERGpISVRERKSGlETrxmS/A4hgujYV07WpnK5PxXRtKlan10b3REVERGpILVEREZEaUhKtZWY2wsy+MbM1Znav3/FECjPrZGZzzWylmWWa2W/9jinSmFmsmS0zs3f9jiWSmFlzM3vTzFYF//yc7ndMkcLMxgf//5RhZq+ZWaLfMfnFzKaZ2TYzyyhR1tLMPjGz1cH3FrV9XiXRWmRmscAzwEVAD+CXZtbD36giRiHwO+fcycBg4E5dm3J+C6z0O4gI9CTwoXOuO9AbXSMAzKwDcBfQ3znXC4gFrvY3Kl+9AIwoU3YvkOKc6wakBL/XKiXR2jUQWOOcW+ecKwBeB0b6HFNEcM5tcc4tDX7OwfuLsIO/UUUOM+sI/AyY4ncskcTMmgJnA1MBnHMFzrnd/kYVUeKAJDOLAxoCm32OxzfOuc+A7DLFI4EXg59fBEbV9nmVRGtXB2Bjie9ZKFGUY2adgdOARf5GElGeAP4HCPgdSITpCmwHng92dU8xs0Z+BxUJnHObgEeBDcAWYI9z7mN/o4o4bZ1zW8D7hzzQprZPoCRauyxEmYY/l2BmjYGZwDjn3F6/44kEZnYxsM05t8TvWCJQHNAX+Jdz7jRgH3XQJReNgvf3RgJdgJ8AjczsWn+jOvooidauLKBTie8dOYq7V8oys3i8BDrdOZfsdzwR5EzgUjNbj3cL4Dwze8XfkCJGFpDlnCvutXgTL6kKnA9855zb7pw7CCQDZ/gcU6TZambtAYLv22r7BEqitesroJuZdTGzBLyb/G/7HFNEMDPDu6+10jn3D7/jiSTOufuccx2dc53x/sykOufUogCccz8AG83spGDRMOBrH0OKJBuAwWbWMPj/r2Fo0FVZbwNjgp/HALNr+wRxtX3Ao5lzrtDMfg18hDdSbppzLtPnsCLFmcB1QLqZLQ+W/cE5976PMUl0+A0wPfgP03XAjT7HExGcc4vM7E1gKd7o92UcxTMXmdlrwFCglZllAQ8Ak4AZZnYz3j86rqz182rGIhERkZpRd66IiEgNKYmKiIjUkJKoiIhIDSmJioiI1JCSqIiISA0piYqIiNSQkqiIlGNmTc3sQTM72e9YRCKZkqiIhNIf72H1eL8DEYlkSqIiEsppQD6aYk+kUpqxSERKMbOVQPcyxTOdc1f4EY9IJFMSFZFSzGwA3moymcAjweItzrnv/YtKJDJpAnoRKeu/eMv4PeWcW+h3MCKRTPdERaSsnkAC3uogIlIJJVERKasv4IDlP7ahyNFOSVREyjoNWOuc2+t3ICKRTklURMrqgR5tEakSDSwSkbJ2A33N7EJgD7DaObfT55hEIpIecRGRUsysFzAVOBVIBM5yzs33NyqRyKQkKiIiUkO6JyoiIlJDSqIiIiI1pCQqIiJSQ0qiIiIiNaQkKiIiUkNKoiIiIjWkJCoiIlJDSqIiIiI1pCQqIiJSQ/8fSs0IMBk8e5sAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "def euler(f,u0,t0,t_max,dt):\n", " u=u0; t=t0; u_all=[u0]; t_all=[t0];\n", " while t