{ "nbformat": 4, "nbformat_minor": 0, "metadata": { "colab": { "provenance": [] }, "language_info": { "name": "python" }, "kernelspec": { "name": "python3", "display_name": "Python 3" } }, "cells": [ { "cell_type": "markdown", "source": [ "#Solving Differential Equations with the Euler-Cromer Method" ], "metadata": { "id": "3xOVbovRVSRj" } }, { "cell_type": "markdown", "source": [ "You will encounter differential equations that cannot be solved analytically. However, there are numerical methods that can be used to solve them with a computer. For many second-order, ordinary differential equations, a fairly simple approach called the Euler-Cromer Method is sufficient. This method was used frequently in Physics 231, but the name may not have been mentioned. (Note that it is referred to as the \"last point approximation\" or LPA in the refrence below.)" ], "metadata": { "id": "bqcQAsODVSRl" } }, { "cell_type": "markdown", "source": [ "##Notation" ], "metadata": { "id": "5N1H1j33VSRm" } }, { "cell_type": "markdown", "source": [ "We will use shorthand notation for time derivatives where \n", "\n", "$$\n", "\\dot{x} = \\frac{dx}{dt} \n", "$$\n", "and \n", "\n", "$$\n", "\\ddot{x} = \\frac{d^2x}{dt^2}.\n", "$$ \n", "\n", "In other words, the number of dots indicates the number of time derivatives. \n", "We’ll be finding the value of $x(t)$ at discrete values of\n", "$t$. Label these times with integer subscripts as \n", "\n", "$$\n", "t_1 = t_0 + \\Delta t \\\\\n", "t_2 = t_0 + 2\\Delta t \\\\\n", "\\vdots \\\\\n", "t_i = t_0 + i\\;\\Delta t\n", "$$ \n", "\n", "where $t_0$ is the initial time and $\\Delta t$ is the time interval between the subsequent solutions. We will also label the values of the function at these time with integers as \n", "\n", "$$\n", "x_0 = x(t_0) \\\\\n", "x_1 = x(t_1) \\\\\n", "\\vdots \\\\\n", "x_i = x(t_i).\n", "$$ \n", "\n", "Similarly, the derivatives at step $i$ are \n", "\n", "$$\n", "\\dot{x}_i = \\dot{x}(t_i)\n", "$$ \n", "\n", "and \n", "\n", "$$\\ddot{x}_i = \\ddot{x}(t_i).$$" ], "metadata": { "id": "Tbbe5pgZVSRp" } }, { "cell_type": "markdown", "source": [ "##The Problem" ], "metadata": { "id": "gSbnGQpmVSRr" } }, { "cell_type": "markdown", "source": [ "The form of the second-order differential equation to be solved is \n", "\n", "$$\n", "\\frac{d^2x}{dt^2} = f\\left(x, \\frac{dx}{dt}, t\\right),\n", "$$ \n", "\n", "or, in the more compact notation, it is \n", "\n", "$$\n", "\\ddot{x} = f\\left(x, \\dot{x}, t\\right),\n", "$$ \n", "\n", "where the initial values of the variable ($x_0$) and its first derivative ($\\dot{x}_0$) are known at the initial time ($t_0$).We want to find $x$ and $\\dot{x}$ at later times." ], "metadata": { "id": "9enCbmy1VSRs" } }, { "cell_type": "markdown", "source": [ "##Euler Method" ], "metadata": { "id": "mxj70-LOVSRt" } }, { "cell_type": "markdown", "source": [ "The simplest numerical method to solve differential equations is the Euler Method. The derivatives at step $i$ can be approximated as \n", "\n", "$$\n", "\\dot{x}_i = \\frac{dx}{dt}\\bigg|_{t_i} \\approx \\frac{x(t_{i+1}) - x(t_{i})}{\\Delta t} = \\frac{x_{i+1} - x_i}{\\Delta t}\n", "$$ \n", "\n", "and \n", "\n", "$$\n", "\\ddot{x}_i = \\frac{d^2x}{dt^2}\\bigg|_{t_i} \\approx \\frac{\\dot{x}(t_{i+1}) - \\dot{x}(t_{i})}{\\Delta t} = \\frac{\\dot{x}_{i+1} - \\dot{x}_i}{\\Delta t} .\n", "$$ \n", "\n", "Suppose that $x$ and $\\dot{x}$ are known at some time. Rearranging the equation for the first derivative, the approximate the solution after moving a step forward in time is \n", "\n", "$$\n", "x_{i+1} \\approx x_i + \\dot{x}_i\\Delta t.\n", "$$ \n", "\n", "To find the solution at another step forward in time, we will need to know the first derivative at that later time. Rearranging the equation for the second derivative, the approximate first derivative after moving a step forward in time is \n", "\n", "$$\n", "\\dot{x}_{i+1} \\approx \\dot{x}_i + \\ddot{x}_i\\Delta t \n", "= \\dot{x}_i + f(x_i,\\dot{x}_i,t_i)\\Delta t,\n", "$$ \n", "\n", "where we insert the function for the second derivative. These final two equations can be used iteratively to find the approximate solution at later times. Unfortunately, this method is not very accurate after a while. (For classical mechanics problems, this method often doesn't conserve energy.)" ], "metadata": { "id": "TWEoZPnSVSRu" } }, { "cell_type": "markdown", "source": [ "##Euler-Cromer Method" ], "metadata": { "id": "5PEtZEsBVSRu" } }, { "cell_type": "markdown", "source": [ "A simple modification to the method described above greatly improves the accuracy of the numerical solution. First, calculate the value of the second derivative at the current time, \n", "\n", "$$\n", "\\ddot{x}_i = f(x_i,\\dot{x}_i,t_i),\n", "$$ \n", "\n", "and use it to approximate first derivative after a step forward in time, \n", "\n", "$$\n", "\\dot{x}_{i+1} \\approx \\dot{x}_i + \\ddot{x}_i\\Delta t.\n", "$$ \n", "\n", "Next, the value of the varible at a later time is approximated by \n", "\n", "$$\n", "x_{i+1} \\approx x_i + \\dot{x}_{i+1}\\Delta t,\n", "$$ \n", "\n", "where the first derivative at the later time is used. \n", "Note the slight difference between this approach and the Euler Method. The current value of $\\ddot{x}$ (which may depend on the current values of $x$, $\\dot{x}$, and $t$) is used to calculate a new value of $\\dot{x}$. However, the new value of $\\dot{x}$ is used to calculate the new value of $x$." ], "metadata": { "id": "qk70KD0hVSRw" } }, { "cell_type": "markdown", "source": [ "##Implementation in Python" ], "metadata": { "id": "GMjacrFvVSRx" } }, { "cell_type": "markdown", "source": [ "As an example, we can numerically solve the differential equation \n", "\n", "$$\n", "\\frac{d^2x}{dt^2} = -kx,\n", "$$ \n", "\n", "for $k=5$ and with initial conditions $x=10$ and $\\dot{x}=0$ at $t=0$. In the program, the first and second derivatives of $x$\n", "will be called $xd$ and $xdd$, respectively. (Thinks of each \"$d$\" as a dot in the shorthand notation for derivatives.) The first part of the program prepares for the calculation. The pylab library is imported so that the results can be plotted. Constant(s), initial values of $x$, $xd$, and $xdd$, time step $dt$, and the final time $tf$ are set. Three lists ($t\\_list$, $x\\_list$, and $xd\\_list$) are started with the inital values of $t$, $x$, and $xd$." ], "metadata": { "id": "IWtckzAfVSRy" } }, { "cell_type": "code", "source": [ "import pylab as pl\n", "\n", "# constants\n", "k = 5\n", "\n", "# set inital values, time step, & final time\n", "t = 0\n", "x = 10\n", "xd = 0\n", "dt = 0.001\n", "tf = 10\n", "\n", "# start lists with the initial values\n", "t_list = [t] \n", "x_list = [x]\n", "xd_list = [xd]" ], "outputs": [], "metadata": { "id": "fUgaH7itVSRy" }, "execution_count": 1 }, { "cell_type": "markdown", "source": [ "The main part of the program is the loop that calculates $x$ and its first derivative ($xd$) at time steps of $dt$ using the Euler-Cromer method. The results for each step are appended to the lists. For this example, the calculation continues until the time reaches $tf$." ], "metadata": { "id": "G2l8hGNXVSR0" } }, { "cell_type": "code", "source": [ "while (t < tf):\n", " # calculate new values\n", " xdd = -k*x # Calculate d2x/dt2 using the current x & dx/dt\n", " xd = xd + xdd*dt # Use d2x/dt2 to update dx/dt\n", " x = x + xd*dt # Use the updated dx/dt to update x\n", " t = t + dt # Advance t by a step dt\n", " # append new values to lists\n", " t_list.append(t)\n", " x_list.append(x)\n", " xd_list.append(xd)" ], "outputs": [], "metadata": { "id": "DlXvSR9gVSR0" }, "execution_count": 2 }, { "cell_type": "markdown", "source": [ "The results were stored in lists so that they can be used. After the calculations are complete, $x$ and $\\dot{x}=dx/dt$ can plotted versus $t$ as shown below." ], "metadata": { "id": "9RILCwGBVSR1" } }, { "cell_type": "code", "source": [ "# plot the x and dx/dt vs. time (from the lists)\n", "pl.figure()\n", "pl.plot(t_list,x_list, ls='-', c='b')\n", "pl.xlabel('t')\n", "pl.ylabel('x')\n", "\n", "pl.figure()\n", "pl.plot(t_list,xd_list, ls='-', c='b')\n", "pl.xlabel('t')\n", "pl.ylabel('dx/dt')\n", "\n", "pl.show()" ], "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deZRV1ZX/v5sqCmSGohhkhqpiHoQSBIOiYtSo0I5tlm2rMU38NemkO/mtxPx6/TpZ3atX59dD7E5MR01ixwwmMUYDaY3EMVFxoCgpKGrgPUpkVAoQUEAoqP37Y7/TVRSvpvfuveece/dnrbfuq/vuu3e/W/vc7zn77HMOMTMURVEUpaf0sm2AoiiK4icqIIqiKEpOqIAoiqIoOaECoiiKouSECoiiKIqSE4W2DYiS4cOH88SJE22boSiK4hUbN248wMwl7fcnSkAmTpyIyspK22YoiqJ4BRG9m22/hrAURVGUnFABURRFUXJCBURRFEXJCRUQRVEUJSdUQBRFUZScsCogRPQIEe0nopo2+4YR0XNElMpsh3bw3Tszx6SI6M7orFYURVEA+y2QHwG4ut2++wC8wMxlAF7I/H0WRDQMwNcBLAKwEMDXOxIaRVEUJRysjgNh5j8S0cR2u1cCWJZ5/yiAlwF8td0xVwF4jpkPAQARPQcRop+HYed3vgP07g386Z8CQ1Wm8qalBVi7FqitBZYsAS69FCCybZX/7NwJPPkkUFQE3HorMHy4bYv8hxl4+mlg82Zg0SLg8svVV9vi4kDCkcy8L/P+PQAjsxwzBsCuNn/vzuw7ByJaBWAVAIwfPz4ngx56CNi6Ffj614FnngEWLMjpNAqAY8eAlSuBF15o3feZzwAPPwwUFNizy3fWrAFuv13uLwB84xsi0hddZNUsrzlxArjxRuDZZ1v33X478KMfAYUuPjktYDuE1Sksq13lteIVMz/MzBXMXFFScs5I/G6xZQuwYQNw3nnAtdcCe/fmY1FyYQb+7M+Al14CHnwQOHwY+NrXgEceAf7v/7Vtnb9UVkrreOZMIJ0GqqqAQYPEV3ft6vr7SnY+8xlg3Trg298WX/3GN4Cf/Qz4yldsW+YQzGz1BWAigJo2fzcAGJ15PxpAQ5bvfBrAQ23+fgjAp7u61oIFCzgftm5l7teP+cYb8zpNYvnJT5gB5n/917P333MPc69ezJWVduzymVOnmOfMYR4zhvnAgdb9qZT46jXXMLe02LPPV371K/HVf/iHs/evXi37X3vNjl22AFDJ2Z7f2XZG+coiIP8C4L7M+/sA/HOW7wwD8A6AoZnXOwCGdXWtfAWEWRwKYP7DH/I+VaI4eZJ5/HjmCy9kPn367M8OH2YeOZL5kkvs2OYzDz0k/vjkk+d+9m//Jp+98EL0dvlMczNzaakIc3Pz2Z99+CHz2LHMCxcmS5idFBBIp/c+AM2Qfox7ABRDsq9SAJ43wgCgAsAP2nz3MwDSmdfd3bleEAJy/DjziBFSs1O6z4MPiretW5f983//92TW7PKhuZl58uSOH2YnTkjLZOnS6G3zmR/9SHxxzZrsnxtfTpIwOykgUb+CEBDm1lZITU0gp4s9LS3M06dL66OjWttHHzEXFzPfcEO0tvnM44933PowfOtbcsymTdHZ5Tvz5zPPnt2xr544wTxqFPPVV0drl006EhCnO9Fd5d57Ja33hz+0bYkfvP46UFcHfO5zHadA9u8P3H038NvfAvv3R2ufrzzyCDBunGS1dcSddwJ9+kiWm9I1VVXy6sxX+/YFPvtZ6WDfvTta+1xDBSQHhg8Hrr8eeOwx4PRp29a4zw9+AAwYIJlCnXH33XI/f/rTaOzymb17gd//HrjjDqBXJ6V42DDgllvknn78cXT2+cr3vy/Zlrff3vlxd90lWYU//nEkZjmLCkiO3HEH8P77wHPP2bbEbU6dAn79a+Dmm0VEOmPGDGDhQhWQ7vDYYzIg88//vOtjb78dOHpUBEfpmNOngSeekBbdkCGdHztlCnDJJeqrKiA58qlPAYMHy8NR6ZgXX5SH1003de/4W24B3n4b2LEjVLO858knZUDr1KldH3vFFTKDwhNPhG+Xz7z6KnDgQM98ta4OaGgI1y6XUQHJkaIi4JprJGbf0mLbGnd56ilpeSxf3r3j/+RPZLtmTXg2+c7+/cAbbwArVnTv+N695b6uWQOcPBmubT7z5JPSv3F1+9n5OsD0Pf3mN+HZ5DoqIHmwYoUU5rfesm2Jm5w5I4Xr2mulYHaH0lJg1iwRHiU7Tz8t8ffuCggA3HCDtARffTU8u3yGWXzuqqu6DrUaxo2TVqAKiJITV18t8zetXWvbEjfZtEkE9vrre/a9FSvkQXfkSDh2+c7atfLwmju3+9+57DJpiaxbF55dPlNbKxlVPfXVlSulNXjgQDh2uY4KSB4MHQpcfLF2TnbE88/LtrvhK8OVV0rr5Q9/CN4m32lulsSNa6/t2aywAwaIr6qAZMf46pVX9ux75vgXXwzWHl9QAcmTyy+XvPEPPrBtiXs8/zwwezYwMtt8yp2weLGkUppCrbRSWSkz7vZUlAEJz2zeDOzb1/WxSeP554GyMqCnE3ZXVMjElW1nl04SKiB5cvnlEj/V2vLZnDgBvPJKbg+6Pn0kRVIF5Fxeekm2l17a8+9edZVsNfX8bJqbgZdfzs1XCwslPJhUX1UByZNFi6S2nNQmbEesXy8ZP1dckdv3ly+XFMk9e4K1y3deekladbksFjV3roRd//jH4O3ymbfeAj76KDcBAeR7jY3yShoqIHlSVAQsXaoC0p6XXpIEg0suye37l18uW33YtXLyJPDaa1LjzYVevaQf5JVXgrXLd0yrbtmy3L5vfDWJUQgVkABYtkxWLExqJkY21q8H5s0DBg7M7ftz5sj8WOvXB2uXz2zYIKHBXB90gFR2tm3T+cbasn69LMY1bFhu3582TVp2SfRVFZAAWLJEtm+8YdcOVzh9WsICixfnfo7CQgkPvvZacHb5jmk55NqqA0RAAB0PYmhpkck+TRnOhV69xNdVQJScqKiQcM3rr9u2xA22bJFMoXwKJSDhlupqiU8rUkEpLweKi3M/x4IFMqhTw1hCQ4MsV5tPZQcQX6+tBQ4dCsYuX1ABCYD+/aWDUlsggqmJ5SsgS5ZIDfHNN/O3yXeY5T4sWpTfeYqK5BzaAhGC8tWLL5Zt0p4BTgoIEU0lok1tXkeJ6K/bHbOMiI60OebvbNkLSA3mrbdkAFzSWb8eGD265zn17bnoIhksl8TQQHt27pTZn/MVEEDua3W1zosFiG8NGyYtu3y48EKJQiQt5OqkgDBzAzPPY+Z5ABYAOA4g2+xIr5jjmPnvo7XybBYvllBLTY1NK9zAxJR7MlI6G0OGyBTvKiCtrbCLLsr/XBdeKGMfqqvzP5fvvP66lN18fbV/f0kaSZqvOikg7bgCwHZmfte2IZ1hYqhJ7wd5/33gnXfyjykbLrwQ2LhRQjhJ5s03pe9izpz8z3XhhbLdsCH/c/nM4cMy1ihIX62qStbs3D4IyG0Aft7BZ4uJqJqIfkdEM7MdQESriKiSiCqbmppCM3LSJBncVVkZ2iW8oKpKtuYhlS8LFgBNTbp06JtvAvPny4SI+TJuHDBihArI22/LNkhfPXoU2L49mPP5gNMCQkRFAFYA+FWWj6sATGDmuQC+AyDrpMrM/DAzVzBzRUlJSYi2SgE3TplUjIDMmxfM+RYskO3GjcGcz0eam+X3B9H/AYivXnihVnaMr15wQTDnS6KvOi0gAK4BUMXM77f/gJmPMvNHmffPAOhNRDlM8BAc8+dLCuupUzatsEtVlUxKN2hQMOebO1fy7JP8sKurk/XMg6opA5J6XleX7BTpqippjQVVr5w5U7LcVEDc4dPoIHxFRKOIpOuLiBZCfsvBCG07hwsukNri1q02rbBLVZUIaVD06ycFM0mFsj2bNsk2qJoyIGLU0tJaC08iQftqUZFUeJLkq84KCBH1B3AlgCfb7LuXiO7N/HkzgBoiqgbwbQC3MdvtajXOmNRCeeiQrGUeZKEEJDSQ5I70TZtkws6ysuDOaVozSW3ZffSRDCIMy1eT0pHurIAw8zFmLmbmI232PcjMD2beP8DMM5l5LjNfxMzWE+gmT5bQTVIFxPT/hFEok9yRvmmTZF8VFAR3zhEjZKxOUlN5q6ulQhKGryapI91ZAfGRXr0kzJBUAQm6U9KQxM5JA7MISFBJCW2ZOze5AmJ8NQwBAZLjqyogATN/vhTK06dtWxI9VVXAhAn5zdWUjTlzJHNo8+Zgz+sDu3bJapdhCMicOTJ/U3Nz8Od2nY0bZaXM0aODPe/MmTIRaFJ8VQUkYC64QKbcbmiwbUn0BN0paejfH5gyJTmFsi0mLBhWC6S5GaivD/7crmN8Nd8R6O0pKgKmTpVszCSgAhIwZqRwUhzIcPw4kEoFM1I6G7NnJ++eAhK+IpLfHzRz58o2acJ86pSkMIflq3PmJOeeqoAEzLRp0tmZtDmx6uokXh/Ggw6Q86bT0rpLEps2yUR//fsHf+7ycqkxJ60fJJWSEHOYvrpzJ3DkSNfH+o4KSMD06SMFM2m1ZSOYs2aFc/45cyQ1srY2nPO7Slgd6IBMizJzZvIEJApfbXudOKMCEgKzZiXDedpSUyPiOWVKOOc3tcWkhAYASQfdsSO8UAsgYawk3VNAfLWgQKIFYZAkX1UBCYHZs4HGRlmVLynU1ADTp0sGShhMmSKD6ZLUsjOtrbBqyoCI03vvJWuN9JoaiRL06RPO+ceNAwYPToavqoCEgCnwSZrSpKYm3AddQYGEW5JQKA1GQGZmnWc6GJLYkb5lS7i+apIekuCrKiAhYJqwSQljHT4so8TDLJRAsrJbAKmA9O0LTJwY3jWMONXVhXcNlzh2TKIDYfuqEZC4T7+jAhICkyYlK9xiWlphZbUYZs+WUEtSwi1bt0pYMMgpTNozYoQs6ZqU5ASTLRhFZefIERkIGmdUQELAhFuS0gIJO6vFkLSWXW2tLOkbJkQiUkkREPXVYFEBCYkkZWJt2QIMHCidh2EyfbpskxBuOXpUaq9h9n8YZsxIloCEmS1oSIqvqoCExOzZkt1y4IBtS8LHdKAHPS1Ee0aPltmO414ogdYHetgtEHONAwdkxuO4U1MjvzfMsCAgYcERI+LvqyogIWFqIHGfZ4g5/Awsgwm3xL1QAtFkYBmMSCWhFRKVrwIyziTuvqoCEhJmkFLcBaSpCTh4MJqaMpCMQgm0ZmBNmhT+tZIiIEePAnv2ROerprIT50wsZwWEiHYQ0RYi2kRE56ybRsK3iShNRJuJKIR5YHNn/Hh5AMRdQMzvMy2usJk+Hdi3L/7zDNXWts6rFjZjxkgfVtwFxMyQHdYI9PZMny5T8cc5a9BZAclwGTPPY+aKLJ9dA6As81oF4HuRWtYFBQUyrXPca8umUE6dGs31ktI5uXVrNOErQEKDSehIN5WdKAUEiLevui4gnbESwI9ZeAPAECIKeHmY/Jg2LRktkL59pcUVBUkolCYDK6pQC5AcASksDD8Dy5AEX3VZQBjA74loIxGtyvL5GABth+nszuw7CyJaRUSVRFTZFHGaybRpwDvvAB9/HOllI6WhQeYV6hWRJ02aJFOQx7lQmt8WtYC89x5w6FB014yahgZg8mSZhTgKxo4FBgyIt6+6LCCfYOb5kFDVaiK6JJeTMPPDzFzBzBUlJSXBWtgF06dLB1oqFellI6W+PrqQACA1yPLyeBdK4y9RhQWBVrGK832N2leJ4p/04ayAMPOezHY/gKcALGx3yB4AbYeujc3scwbjrHF1oJMnpYUV5YMOiH8qbzotD58oMrAMcQ+3nDkjwhylgADx91UnBYSI+hPRQPMewCcBtB/XvRbAn2eysS4CcISZ90VsaqeUl8uDIK79IOm0LPJko1DGOTSYSrVm8UXF+PEyQnvbtuiuGSU7dshStjYqO3v2SL9WHHFSQACMBPAqEVUDeAvA08z8LBHdS0T3Zo55BkAjgDSA7wP4Szumdsx558lMqnEVEPO7bBTKlpb4PuzSaaC0NNprFhTINeN6T6NO4TXEfUBxSMv/5AczNwKYm2X/g23eM4DVUdqVC3HOxIo6hdfQNtwS5mp9tkingVtvjf665eXx9VVblZ22YeyF7YPwMcDVFkhsMALS0mLbkuCpr2/NNImS8nLZxrG2fOiQvKJugQByX9Np4PTp6K8dNvX1wPDhQHFxtNedMkVad3H0VUAFJHSmTQNOnIjnugANDdHX6AAJDY4bF89CmU7Ltqws+mtPnQo0NwPvvhv9tcOmoSH68BUgKcOTJ8fTVwEVkNCJa3YLc/RpkW0pL49noTQpvLZaIEA872t9vZ3KDiCVgbim8quAhExcC+X770tmic1CuW1b/CaqMym8kydHf23jq6ZvKy6Y+ahsVnZSqXiGsVVAQmbECJmoLm41kKjnFWpPebmsxX7woJ3rh0UqJeG5KFN4DcOHA0OHxq+yYysDy1BeDhw/Duzda+f6YaICEjJErTWQOGErA8tgastxu682UngNxlfjKiA2W8tA/HwVUAGJhDjGQOvrgX79JAvLBqZQxu1hl0rZ6UA3lJfHL4RVXy+d2VGO7G9LXMPYgApIJJSVtY6EjQuplNSUo5pEsT2TJkl6ZJyE2WYKr2HqVGD3buDYMXs2BE0qJem0hZZGvY0dKyFJFRAlJ8rKpAPtnXdsWxIc6bTdmrKpUcapUG7fLlvbLRCgNZ04DtgMCwJSySotjVdlx6ACEgFxi4GeOQM0NtotlED84vU2U3gNccvEYrYvIED8fNWgAhIBcYvX79olA85cKJSpVHxSeU0Kb1QLHmUjbr763nsSjrPtq2Vl0sKM2yh/FZAIKC6W9Mi4tEBMeMOFQhmn9MhUqjVebot+/eI1yt/myP62lJeLeMRtlL8KSETEKRPLFQGJWyqvC6EWIF6ZWK75alyE2aACEhFxE5C+fYHzz7drR9zCLbZTeA1xiten05J9NX68XTvi1g9qUAGJiPJy6TuIwyJI6bTE6W2l8BrGjZNFkOJQKD/4QEbV264pA2LD4cPxWB89lZJsPVspvIYRI4BBg+IjzAbnBISIxhHRS0RUS0RbieiLWY5ZRkRHiGhT5vV3NmztCWVl0tlrUjV9xpVQS69erXNi+Y4LKbwG87+NQyqv7XRzA1F8fLUtzgkIgNMAvszMMwBcBGA1Ec3IctwrzDwv8/r7aE3sOXFpwra0yMPOBQEB4lMoXUjhNcRFQFxJ4TXEcUoj5wSEmfcxc1Xm/YcA6gCMsWtV/sQlXr93r4ThXCqU27fL2BSfMQ9rmym8hsmTpcbsu4Ds3w98+KFbvvruu/EIYxucE5C2ENFEABcAeDPLx4uJqJqIfkdEMzs5xyoiqiSiyqamppAs7ZrBg4GSEv9rIK5ktRjKyuKxCJJJ4T3vPNuWSILE2LH+C4iLvhqXMLbBWQEhogEAfg3gr5n5aLuPqwBMYOa5AL4D4DcdnYeZH2bmCmauKCkpCc/gbhCHTCzXCqWxw/dC6Uqs3lBaGh8BceW+xsVX2+KkgBBRb4h4/IyZn2z/OTMfZeaPMu+fAdCbiIZHbGaPiYuA9O4tGVAuYEI+vhdKMzmlK0yZ4v89Tadlws0JE2xbIsTFV9vinIAQEQH4IYA6Zv5WB8eMyhwHIloI+R3OLy1UViZ9CD7PdJpOS4y8oMC2JcL550vIxefa8uHDwIEDbglIaan0IRxt3/b3iFQKmDhRKjwuUFwsoWyffbU9zgkIgIsB3AHg8jZpup8ionuJ6N7MMTcDqCGiagDfBnAbs/szIpmmtM8O5FJWCyCpvFOm+H1PXUrhNcQh3OKar5p5zny+p+2xPLzmXJj5VQDUxTEPAHggGouCo20q79y5dm3JBZMWuWyZbUvOxvdC6VIKr6FtKu8FF9i1JReMry5ebNuSsyktBTZutG1FcLjYAoktvufXv/++GzObtqe0VATE/TZodlxK4TUYW3z11YMHgSNH3PPVKVMkY7C52bYlwaACEiEDBwIjR/pbKF3LwDJMmQKcOAHs22fbktxIpYAxY2QmXFcYMAAYNcpfX3WxVQeIPadPAzt32rYkGFRAIsbnTCxXBcT3lp1rKbwGn1N51VejQQUkYnwvlC6lRRp8T490LYXX4HNyQjotCRaTJtm25Gx899X2qIBETGmpv6m86bRbaZGGCRNktlUfH3ZHjgBNTe62QPbulUW7fCOVEr8oKrJtydmMHi2zDfjoq9lQAYkY86DwsQbiWlqkobBQHhY+FkpXQy1Aq02NjXbtyAVXfbVXLxlH5WP5z4YKSMQYp/atH8S1mU3bYzKxfMMHAfFVmF28p4DfYez2qIBEjK+F0tW0SIOJ1/uWyutiCq/B11TeQ4dkgS4Xw4KA3NfGRlkawXdUQCJm0CBZncy3QulyTRkQu44c8W8VvVRKpmPp39+2JecydKhMv+Gbr7qawmsoLZUp3ffutW1J/qiAWMDHVF4fBATw72Hnagqvwcdwi+u+GqdMLBUQC/haKIncS4s0+FooXU3hNaivBo+vlZ1sqIBYoKwM2LPHr/TIdBoYPx7o08e2JdnxcRW9o0dlxluXWyBTpsio6ZMnbVvSfYyv9u1r25LsjB8vmYO+VXayoQJiAR9nOnU5qwWQh8WYMX4JiOuhFkBsYwbeece2Jd3H9VZdYaGMp/LJVztCBcQCPqbyui4ggH+pvK6tmJcNH8Mt6bSbWW1t8c1XO0IFxAK+FcoPPpA0XtcFxLepN1xO4TX45quHD4uvuizKgL9p5+1RAbHA4MFASYk/hdLUlFwXELOK3ocf2rake6RSMrWFiym8huHDJfXcN191WZQB8dWjR2UlSp9xVkCI6GoiaiCiNBHdl+XzPkT0y8znbxLRxOitzB2fUnl9iNUD/mViuZ7CC0higk/hFvXVaHFSQIioAMB3AVwDYAaATxPRjHaH3QPgA2YuBXA/gP8XrZX54VN6pLFz8mS7dnSFb8kJrnf2GnzyVfO/98VXfbmvHdGlgGR5cIOIloViTSsLAaSZuZGZTwH4BYCV7Y5ZCeDRzPsnAFxBRJ0uhesSZWXA7t1+pPKm0+4teJQNn6be+PBDWeHR9RYIIA+7HTv8WEUvnXY/LAjIGBUifyo7HdGdFsjjRPRVEs4jou8A+KeQ7RoDYFebv3dn9mU9hplPAzgCoLj9iYhoFRFVElFlU1NTSOb2HJ9mOvUhAwuQWL0vfUu+hFoAv1bR277dj3vaty8wdqwfvtoZ3RGQRQDGAVgPYAOAvQAuDtOoIGHmh5m5gpkrSkpKbJvzP/iUyuuLgADSCvGhVudDCq/Bp3CLDym8Bl98tTO6IyDNAE4AOA9AXwDvMHPY80jugYiWYWxmX9ZjiKgQwGAAB0O2KzB8KZQm1OKLgPgSr/chhdfgi68ePy4TFKqvRkd3BGQDREAuBLAU0qH9q1CtkmuWEdEkIioCcBuAte2OWQvgzsz7mwG8yOxPVvWQIZIi6boD+ZLCaygtlb6ljz+2bUnnpFLAqFHAgAG2LemaUaOk/8t1XzXhYB9EGRA7m5oknddXuiMg9zDz3zFzMzPvY+aVOPdhHiiZPo3PA1gHoA7A48y8lYj+nohWZA77IYBiIkoD+BKAc1J9XceHVF6fasqA2OnD1Bs+pPAaTCqv6wLiU78S4F/WYDYKuzqAmSuz7PtJOOacdY1nADzTbt/ftXn/MYBbwrYjTEpLgZdftm1F5/hcKKdPt2tLZ6RSwDXX2Lai+5SWArW1tq3oHF8GERra+uoFF9i1JVecHAeSFMrKgF27gBMnbFvSMek0MHIkMHCgbUu6hw+pvB99BLz3nj8tEEAedo2NwJkzti3pmHQaGDZMFsLyAR98tStUQCziQyqvTxlYgB9Tb/jWqgPE1lOnpH/JVXzz1YEDZXVSn0NYKiAW8SG7xbdCSeT+pIo+pfAafPDV7dv9CV8Zpkxxvx+0M1RALOL6WJATJ2ThK58EBHB/7ibz//bpYee6gJw6Bbz7rvpq1KiAWGToUKC42N1CaUJrvhXKsjLJwnJ16o10WlJjfelXAmQqmz593PXVd98FWlr881WTdu5yP2hnqIBYxuVUXh9j9YDYe+aMPFRcxLewIAD06uV2aNC3dHODCWP62gpRAbGMy/n1vhZK18MtvszC2x4ffNW3++q6r3aFCohlTCqviyOn02kJsfmSFmkwtToXW3bHjgH79vnVgW4w8fqWsCcyyoHt22UG3hEjbFvSM1zvB+0KFRDLlJbKyGkXU3l9DLUAMm5lwAA3a3W+1pQBsfnECRFA1zC+6s+CDoLr/aBdoQJiGZebsD7NbNoWM/WGi7U6H1N4DS77qo8pvIayMjfvaXdQAbGMq+GWkydl/Qcfa8qAu/F683/28b66KiBnzkgL3sd7Crhb2ekOKiCWGTpUpl9wrVDu2OFnWqTBpPKePm3bkrPxbWqYtowbB/Tu7Z6v7t4t40B89dXSUvenNOoIFRAHcDGV1+dYPdC6ip5rqby+9isBQGGhLMXqmoD4Nolie0wUwsV+0K5QAXEAF8MtcRAQwL37mkr52f9hUF8NHld9tTuogDhAWZn0N7iUyptOy6SEw4fbtiQ3XOxbOnbMrxXzsmEExKWl27ZvB4qKZLS8jxhfVQFRcsKk8rq0CJKvaZGGUaNkXIBLhdLnDCxDaalMR79/v21LWkmngcmTgYIC25bkhukHdamy012cEhAi+hciqieizUT0FBEN6eC4HUS0hYg2EdE5C175hotNWJ9j9YCbq+iZB4TvAgK4dV9991XAPV/tLk4JCIDnAMxi5jkAtgH4WifHXsbM85i5IhrTwsO1cEtzs2RhxaFQunJPAb9TeA2uCQiz32NADC4m0nQHpwSEmX+fWQ8dAN4AMNamPVFhVlFzpVDu3CkZTD4/6ACx36VU3lTKv1l42zNhgoSKXPHV99+XvqU4+KqrUxp1hlMC0o7PAPhdB58xgN8T0UYiWtXZSYhoFRFVElFlU1NT4EYGhUs1EN+zWgxlZdKa2rnTtiWC7xlYgHRWT5jgjoD4nsJrKCtzrx+0O0QuIET0PBHVZHmtbKoUA+cAABlsSURBVHPM3wI4DeBnHZzmE8w8H8A1AFYT0SUdXY+ZH2bmCmauKCkpCfS3BIlLMVBTKH0XENfCLakUUF5u24r8cclXt22Tre/C7OukipELCDMvZ+ZZWV5rAICI7gJwHYDbmbMnCzLznsx2P4CnACyMyPzQMKm8J0/atkQeDv36SbjFZ1xKjzx6VMItvj/ogNa+JRdSeVMpGeA4caJtS/LDJV/tCU6FsIjoagBfAbCCmY93cEx/Ihpo3gP4JICa6KwMh9JSmTrEhSas7ym8htGjRQhdqNXFIQPLUFoKHDkCHDpk2xJpgUyZIiLiM6Yf1AVf7QlOCQiABwAMBPBcJkX3QQAgovOJ6JnMMSMBvEpE1QDeAvA0Mz9rx9zgcCnc4ussvO1xKZU3bgICuHFft22LR1gQcMdXe4JTus3MWaPuzLwXwKcy7xsBzI3SrihwJZXXzGx63XV27QiK0lJg61bbVrT+X+MgzG0FZNEie3a0tMh9vfJKezYESVkZsH69bSt6hmstkMQybBgwZIj9Gojph4lTra6xUYTRJqmUzGbbr59dO4Jg0iRp3dn21T17JO01Tr7qSj9od1EBcQQiN1J5GxpkO3WqXTuCwqTy7tpl1444pPAa+vYVMbQtICYDKy4CUlbmTj9od1EBcQgXYqBxK5SupEfGSUAA9dUwcMVXe4IKiEOUlcn6FadO2bOhoQEYPBgYMcKeDUHiQnrkoUPAwYMqIEGTSklI8Pzz7doRFC4lJ3QXFRCHcCGV12S1+J7Caxg9GjjvPLu1ujhlYBlKS4EDB4DDh+3ZsG2b3NO4+GpxsRv9oD1BBcQhXKgtxyktEgB69ZLMJ5v3NI4CYnzEhJFsEDdfNWnnGsJScsJ2DPT4cckCiUsHuqG83O6DLpUSIZs82Z4NQTNtmmzr6+1cv7lZsuviJCCAVDJs+mpPUQFxiOJi6X+wVVs2141boZw6Veb3am62c/1USiYg7NPHzvXDYPJkGf1tsvaiZscOSc2OU6sOEGHeuVMqcz6gAuIQtlN545bCa5g2TaZ0b2y0c/24ZWABQO/eEhq01QKJWwaWYdo0mWPMlzCWCohj2MxuicvMpu2xGW5hbu3sjRtTp9prgcRZQAB7wtxTVEAco6xMmuc2UnkbGoCxY2Ut8ThhWlQ2CuV778lMvNOnR3/tsJk2TWrKNkb5b9smszcUF0d/7TAxWWUqIEpOmFTeHTuiv3bcsloMgwdLOq+NQmmuaWqWcWLqVKno2PDVOIYFAUk5nzDBXsuup6iAOIatVF5mcdo4CghgL9xSVyfbOAqIzXBLXCs7gNxXbYEoOWErldcMCotbB7rBFMqoF0Gqr5c10OMyWrotxleiFuaPPpK5zeLsqw0NEolwHRUQxxg+HBg0KPoWiHkIxLlW98EHQFNTtNetr5drx2W0dFuKi8Vfo64tG1+NY78SIP5y/Diwe7dtS7rGOQEhom8Q0Z7MglKbiOhTHRx3NRE1EFGaiO6L2s6wsJXKa7Ja4lyrA6J/2NXVxTN8ZbARGqytlW2cBQTwox/EOQHJcD8zz8u8nmn/IREVAPgugGsAzADwaSKaEbWRYWEjlbe+XnL7J0yI9rpRYSPc8uGHUouMs4DYiNfX1ckgxtKsy8/5j82swZ7iqoB0xUIAaWZuZOZTAH4BYKVlmwKjvFwyW6JcWKauThzX97WlO2L8eFnHIspCaVp1ca0pA+Iz+/dHO6liXZ2IR+/e0V0zSkaOlMxBFZDc+TwRbSaiR4hoaJbPxwBou0TQ7sy+cyCiVURUSUSVTVEHwHNkxgzJrY9yTpzaWrluXOnVSx52URbKOGdgGWyEW+rq4i3KRP5kYlkRECJ6nohqsrxWAvgegCkA5gHYB+Df8rkWMz/MzBXMXFFSUhKA9eFjHuQm1hs2x4/LFPIzZ0ZzPVtEXSjr64GCgnisg94RUYdbTp2S8G6cKzuAPwJiJWDBzMu7cxwRfR/Af2f5aA+AcW3+HpvZFwumTpUa89at0VyvoUHSW5NQKB9/XNbR7ts3/OvV14t4FBWFfy1bTJoU7aSK6bS0zuPcAgHkGfDoo9KPNnCgbWs6xrkQFhGNbvPnDQBqshy2AUAZEU0ioiIAtwFYG4V9UdCnj8R4o2qBmOvEXUCmThWhjCpBIe6hFkD6IUpLW8N1YRP3DCyDL5lYzgkIgH8moi1EtBnAZQD+BgCI6HwiegYAmPk0gM8DWAegDsDjzBxRfT0aZsyIrgVSWxvvrBaDKZRRPOxOn5ZU7Dj3fxhmzoyusmP+d3FNNzdE6av54JyAMPMdzDybmecw8wpm3pfZv5eZP9XmuGeYuZyZpzDzP9qzOBxmzpQHUBSTKtbWytiTOIdaACmUUYUGGxtl/ZEkCMisWdKqO3Ei/GvV1Umqedwm/GyPyTKLqhKZK84JiCJEmYkV9wwsw3nnScHcsiX8a9VkAq9xT0wAREBaWqLp9E1CWBAQ8Zg2rdWPXEUFxFHMgyfs0MDJk8nIajHMmhVNoaypkXTMJNzXWbNkG/Z9bWmRPoEkCAgAzJ4dTWUnH1RAHCWqTKxt26RgJuFBB0QXbtmyRZZ9jXuoBZBWXVFR+AKyY4f835IiILNmyfK2R4/atqRjVEAcpW9fSQENuwWSlAwsgwm3hN05WVMjNcgkUFgoD/WwBWTzZtnOmRPudVzBtOxc7gdRAXGYKDKxtm6Vlk5cZ+Ftj3moh/mw+/hjSYAwD4AkEEVo0IRzktCvBLT6qsthLBUQh4kiE2vLFglBRDGwzgWiCLfU10sCRFJaIEA04ZbNm6VVPmBAeNdwifHj5be63JGuAuIwM2bIeIIwM7Gqq4G5c8M7v2uYcEuYtTpz7qS1QIBwW8xbtiQnfAVIZGDmTG2BKDliHuwm9hs0R4/KHFhJEhBAWgZh1upqaqSVE8c1uzsi7EysEyekNZ6kVh3QmokV9Uqa3UUFxGGmTpVpTTZtCuf8pmaTNAGZNUvW6QhrCvItWySHP67TjWcj7HBLba0kPySpBQKIrx48KFPmu4gKiMP07i0OFJaAVFfLNokCAoT3sEtSBpYh7HBL0jKwDMZXXQ1jqYA4zrx5IiBhNGGrq4GhQ4GxY4M/t8uEmd1y+DCwa1ey+j8Ms2fLgz4MX928WWYSmDw5+HO7jAqIkhfz5gFNTcC+fcGfu7paanREwZ/bZcaNE+EMo2VnzjlvXvDndp358yXcsmtX18f2lC1b5GFaUBD8uV1m5Ehg1KjwohD5ogLiOOZBFLQDnTkjhTJp4StABHP+fGDjxuDPXVUl2/nzgz+365jfbO5BUDBLCyRpYUHDggXh+GoQqIA4jon5Bi0gjY2yEmESBQSQh92WLcGPsamqkpDgiBHBntcH5syRFkLQArJnj7TCL7gg2PP6wvz5MnPC8eO2LTkXFRDHGTRIBk8FLSBJ7UA3zJ8v4hH0VDFVVclsfQDSRzF9evACUlkp24qKYM/rC/PnSwZaWOn8+aAC4gGmIz1IKitbs7ySyIIFsg0yNHDsmIxCT6qAAPLbgxaQjRulZZPkyg4Q/H0NAqcEhIh+SUSbMq8dRJT1sZn5bEvmuMqo7YyaefNkBtkPPwzunBs2SMihT5/gzukTU6bIWtNBFsrqaonXJ11A9u0LNumjslJShM87L7hz+sS4ccDw4W72gzglIMz8p8w8j5nnAfg1gCc7OfyyzLGxb9jOny8PpqAedi0tUigvvDCY8/lIr14SUw9SQJLcgW4wv/3tt4M5H7M8OE2LMYmYpA9tgXQTIiIAtwL4uW1bXGDhQtm++WYw50ulZBqTJAsIIIWyulrmGwuCqirpPD///GDO5yMmazCoh92uXdKBntT+D8P8+TJA9eRJ25acjZMCAmApgPeZOdXB5wzg90S0kYhWdXYiIlpFRJVEVNnU1BS4oVEwfLiEXIISkA0bZJt0AVmwQOZYCmop1spKKehJG1fTloEDZWmAyoACyyZsk+QWCCC///Rp9wYURi4gRPQ8EdVkea1sc9in0Xnr4xPMPB/ANQBWE9ElHR3IzA8zcwUzV5SUlAT0K6Jn0SLgjTeCOdeGDUC/fslZ2a0jjIAGIcxHjkgNcfHi/M/lOxddBLz+ejAj0isrZQblpE1h0h7TAguqEhkUkQsIMy9n5llZXmsAgIgKAdwI4JednGNPZrsfwFMAFkZhu00WLQL27pVJAPNlwwaJ/xcW5n8unykvB4YNA9avz/9cb70lD0wVELkH+/fLTM/58sYbIh5J7UA3TJgAjB4djK8GiYshrOUA6pk566OSiPoT0UDzHsAnATi85EowXHSRbPOtgZw6JR2cSQ9fARJqWrIkmEL5+utyvkWL8j+X7yxZItt87+vp0+LvF1+cv02+QyT3QQWka25Du/AVEZ1PRM9k/hwJ4FUiqgbwFoCnmfnZiG2MnLlzZY2JfAVk40ZZcvUTnwjGLt9ZskT6QA4ezO88r78uqaaDBgVjl8/MnCl9Ifk+7KqrZWyNEaSks2QJsGOHRCJcwTkBYea7mPnBdvv2MvOnMu8bmXlu5jWTmf/RjqXR0qePhJ1efz2/87zyimxVQARTu83nvra0SKhFw1dCQYG0xPL11ddek622QARzH1xqhTgnIErHLF0qsfYTJ3I/xyuvSOx/5Mjg7PKZigrpC8qnUDY0yDTuKiCtLFkiU2/kM/j1tddkEN24ccHZ5TPz5gF9+6qAKDmybJn0YeRas2tpkUK5dGmgZnlNv37SsjO13Vz4wx9kq626VhYvFn/LNeTKDLz6qrY+2lJUJH2X+fhq0KiAeMQnPiEjqF9+Obfv19YCH3ygAtKepUvlQZdry+7FF6WWXFoarF0+c/HF0rJ78cXcvt/YKLF+FeWzWbpU+jGPHrVtiaAC4hGDB8uAopdeyu37RnhUQM5m+XIZ4Wv6h3pCS4v8Py6/PNkDCNszcKD0gzz/fG7fN99bvjw4m+LA8uWylo9p9dpGBcQzli2T2nIuawOsWydLgiZtWdCuuOQSCQ/k8rCrqQEOHBABUc5m+XIZCPjBBz3/7nPPyboq5eXB2+UzS5bImJjnnrNtiaAC4hnLlgHNzRIf7gmnTklN+aqrQjHLa/r3l4KZS6F84QXZXnZZsDbFgeXLpS+jpyHXM2ck9HXlldqqa0+fPlLhUQFRcmLZMsnEePrpnn1v/XrJqVcByc7y5bLmyv79Pfves89KLVkzhc5l4UIR554+7KqqpNVy5ZXh2OU7V14pY5eCmJUiX1RAPKNfP+CKK4Df/rZncw2tWyedmlpTzs4nPynbdeu6/52jR6VVt2JFODb5TlGRhPaefrpnvvr009LyuOKK8GzzGeOrv/udXTsAFRAvuf56mWeou8uxMgO/+Y10nutI6ewsWCDTsD/Z2Qo07Vi3TsKJK1d2fWxSufFGYOfOnk3v/tRTkn2VxHXlu8OsWcCkST3z1bBQAfGQ666T7dq13Tu+tlaavDffHJ5NvtOrlzzsnn0W+Oij7n1n7VqguFgHEHbG9dfLyPTuPuy2b5cBiDfcEK5dPkME3HST9L8dPmzXFhUQDxkzRlIkf9nhfMVn88QT4nQ33hiuXb5z000yT1h3QgPHjwNr1kj4qqAgfNt8pbgYuPRSEZDuhLGeekq2KiCdc9NN0vr97W/t2qEC4il33CGTzVVXd34cM/Dzn0v4atSoaGzzlaVLJWzy2GNdH7tmjUzTcccd4dvlO7feKi3grtb0ZgYefVQ63ydOjMQ0b1m4UBI3uuOrL74IfO5z+U8Ymg0VEE+57Tagd2/gxz/u/Lg//lHmarr77mjs8pmCAuCuu6RW19WMpz/+sRTgSy+NxDSvue02Gbvw/e93ftyGDTKu5p57orHLZ3r1kjK9bp3M0NsZ3/2utAAHDAjBjuBPqURBcbF03v7Xf3Ues3/oIRnBfuut0dnmM3/xFzIO4ZFHOj5m2zYpuHfdJQVZ6Rzjf4891vnkig89JFmGt90WnW0+89nPSmj6Bz/o+Jh335XW8t13yxiSoFH395gvf1ny5TtyoMZG4PHHxXn69YvWNl8pLZU0yQcekHEz2bj/fklRXb06Wtt8ZvVqqeh897vZP9+1C/jJT0SUNVOwe4wbJwk13/tex3Nj3X+/iMwXvhCSEcycmNeCBQs4blxyCfPo0cxHjpz72Z13Mvfty7x3b+Rmec1rrzEDzP/0T+d+tm0bc+/ezJ/7XPR2+c411zAXFzMfPnzuZ6tWMRcWMu/YEb1dPlNZKb76jW+c+9k770j5v/PO/K8DoJKzPFOtPMgB3AJgK4AWABXtPvsagDSABgBXdfD9SQDezBz3SwBF3bluHAXkjTeYiZhXrz57//PPy3/3q1+1Y5fvXHcdc//+zKlU677Tp5mXL2ceOJB53z57tvlKZSVzr17M99xz9v5XXxUf/uIX7djlOzffLEKxdWvrvjNnmK+9lrlfP+Zdu/K/hmsCMh3AVAAvtxUQADMAVAPokxGJ7QAKsnz/cQC3Zd4/COB/dee6cRQQZuYvfEH+k/ffz9zSwlxVJTW98nLmY8dsW+cnu3YxDxnCPH0687vvMp86JSINMD/0kG3r/OW+++QefvOb4qubNzOPGsU8aRLz0aO2rfOTffuYhw9nLitjbmxkbm5m/tKX5D7/x38Ecw2nBOR/Ln6ugHwNwNfa/L0OwOJ23yEABwAUZv5eDGBdd64XVwE5dYp5xQr5b44YIbW5MWPOrj0rPefll5kHDGAuKmIeNkzu75e/LA8+JTeam5lvuUXuZUmJtEhGjz679qz0nNdeYx40SMKrxcVyf1evDs5XOxKQwuB6UwJhDIA32vy9O7OvLcUADjPz6U6O+R+IaBWAVQAwfvz44Cx1iN69ZQDWT34iOd+TJgF/9VeSqaXkzqWXyqjo//xP4NAhGdxmZgFQcqOwEPjFLySD8LnngPHjpYNdl1jOjyVLJAX6gQeApia5vytWhD+bMYm4hHBioucBZBu69rfMvCZzzMsA/jczV2b+fgDAG8z808zfPwTwO2Z+os15h2eOKc38PS5zzKyubKqoqODKysr8fpiiKErCIKKNzFzRfn9oLRBmzmUtsT0A2k6MPTazry0HAQwhosJMKyTbMYqiKErIuDYOZC2A24ioDxFNAlAG4K22B2TicS8BMFMD3glgTaRWKoqiKHYEhIhuIKLdkA7wp4loHQAw81ZIhlUtgGcBrGbmM5nvPENE52dO8VUAXyKiNKRP5IdR/wZFUZSkE1ofiItoH4iiKErP6agPxLUQlqIoiuIJKiCKoihKTqiAKIqiKDmhAqIoiqLkRKI60YmoCcC7OX59OGQKlSShvzkZJO03J+33Avn/5gnMXNJ+Z6IEJB+IqDJbFkKc0d+cDJL2m5P2e4HwfrOGsBRFUZScUAFRFEVRckIFpPs8bNsAC+hvTgZJ+81J+71ASL9Z+0AURVGUnNAWiKIoipITKiCKoihKTqiAdAERXU1EDUSUJqL7bNsTNkQ0joheIqJaItpKRF+0bVNUEFEBEb1NRP9t25YoIKIhRPQEEdUTUR0RLbZtU9gQ0d9k/LqGiH5ORH1t2xQ0RPQIEe0nopo2+4YR0XNElMpshwZxLRWQTiCiAgDfBXANgBkAPk1EM+xaFTqnAXyZmWcAuAjA6gT8ZsMXAdTZNiJC/gPAs8w8DcBcxPy3E9EYAF8AUJFZwbQAwG12rQqFHwG4ut2++wC8wMxlAF7I/J03KiCdsxBAmpkbmfkUgF8AWGnZplBh5n3MXJV5/yHkodLhmvNxgYjGArgWwA9s2xIFRDQYwCXIrKXDzKeY+bBdqyKhEMB5RFQIoB+AvZbtCRxm/iOAQ+12rwTwaOb9owD+JIhrqYB0zhgAu9r8vRsJeJgaiGgigAsAvGnXkkj4dwBfAdBi25CImASgCcB/ZcJ2PyCi/raNChNm3gPgXwHsBLAPwBFm/r1dqyJjJDPvy7x/D8DIIE6qAqJkhYgGAPg1gL9m5qO27QkTIroOwH5m3mjblggpBDAfwPeY+QIAxxBQWMNVMnH/lRDxPB9AfyL6M7tWRU9mWfBAxm+ogHTOHgDj2vw9NrMv1hBRb4h4/IyZn7RtTwRcDGAFEe2AhCkvJ6Kf2jUpdHYD2M3MpnX5BERQ4sxyAO8wcxMzNwN4EsASyzZFxftENBoAMtv9QZxUBaRzNgAoI6JJRFQE6XBba9mmUCEigsTF65j5W7btiQJm/hozj2XmiZD/8YvMHOuaKTO/B2AXEU3N7LoCQK1Fk6JgJ4CLiKhfxs+vQMwTB9qwFsCdmfd3AlgTxEkLgzhJXGHm00T0eQDrIBkbjzDzVstmhc3FAO4AsIWINmX2/R9mfsaiTUo4/BWAn2UqR40A7rZsT6gw85tE9ASAKki24duI4bQmRPRzAMsADCei3QC+DuCbAB4nonsgS1rcGsi1dCoTRVEUJRc0hKUoiqLkhAqIoiiKkhMqIIqiKEpOqIAoiqIoOaECoiiKouSECoiiWCQzI+5f2rZDUXJBBURR7DIEgAqI4iUqIIpil28CmEJEm4joX2wboyg9QQcSKopFMjMe/3dmfQpF8QptgSiKoig5oQKiKIqi5IQKiKLY5UMAA20boSi5oAKiKBZh5oMAXiOiGu1EV3xDO9EVRVGUnNAWiKIoipITKiCKoihKTqiAKIqiKDmhAqIoiqLkhAqIoiiKkhMqIIqiKEpOqIAoiqIoOfH/Afic2K7B4BzZAAAAAElFTkSuQmCC\n" }, "metadata": { "needs_background": "light" } }, { "output_type": "display_data", "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEGCAYAAACO8lkDAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deXRX13XvvxvEDAYDAgNCTJIQILAZbGM5JDHItVsP2HHs2m1T10O8kjip0yTLdTq99/qaqUmzmqnvxW4c2+s5cVIHx6wExzbgIQzGFjJBAgnpJ0DMIGYxI9jvj/07ldCAftLv3nvOuXd/1tL6Sb/h3q372+d+z9l7n3OImaEoiqIorell2wBFURTFPVQcFEVRlHaoOCiKoijtUHFQFEVR2qHioCiKorQjx7YBQTBy5EieOHGibTMURVG8Yv369QeZObej12IhDhMnTkR5ebltMxRFUbyCiBo6e03DSoqiKEo7VBwURVGUdqg4KIqiKO1QcVAURVHaoeKgKIqitEPFQVEURWmHioOiKIrSjljMc1C6DzPw1lvAe+8BM2cCt90G9NKuQtYcPAj8138BZ84A99wD5OfbtigevPsusHo1MG0acMcdQO/eti2KPyoOCaS5GXjkEeCFF1qeu+UW4Fe/AgYNsmeX76xbB9x5J3DggPz9D/8AvPgicNdddu3ymYsXgc98BnjmmZbnFi4Efv1rYMgQe3YlAe0rJpCnnhJh+Md/BI4eBX70I+CNN4CHH5YRhdJ9duyQ0dfgwcD69UB9vYzI/vRPgQ8+sG2dv/zTP4kwPPkkcOQI8PTTwDvvAJ/6lPpq6DCz9z9z585lJTNWrWImYv7sZy99/utfZwaYX37Zjl0+c/Ei8+23Mw8cyFxb2/L8oUPM48Yxz5jBfPasPft85f33mXv1Yn7ooUuf/+53xVdffNGOXXECQDl3cl8ljoH8zps3j3Vtpa5hBhYsALZvB2pqpJdraG4G5s2TkURdHdCnjzUzvWP5cuDmm4Fvfxv4ylcufe3VVyWs9PTTwKc/bcc+XykrAyorxR+vuKLl+QsXgBtuAHbvBrZuBfr1s2ej7xDRemae19FrGlZKEMuXS1Lv7/7uUmEAgJwc4F/+BWhoAF56yY59vvKNbwBjxgBf+EL71+68E7j+euBrXwPOn4/eNl9ZtQpYsUJCoK2FAZBk9Ne/DuzZc2neTAkWFYcE8f3vA6NHSzK6I267TeLk3/mOxnMzpaICWLkS+NKXOu7BEokYNzQAv/lN9Pb5yve/D4wYIcnojli0SEa6//Zv6qthoeKQEHbvBpYtAx56qPNhOBHwuc8BGzdKUlXpmmeflev56KOdv+dP/gTIywN+/OPo7PKZAwekGukv/xIYMKDj9xABjz8ObNkCrFkTrX1JQcUhITz3nJQFXu4mBgD33w/07w/89KeRmOU1Z88CP/+55BSGDev8fTk5ct3feAPYuTM6+3zlhRckBNdVjuaTn5Tw6LPPRmNX0rAmDkQ0nojeIqLNRLSJiJ5IPz+ciN4korr045W2bIwTL70kyegpUy7/vmHDgLvvlvc3N0djm68sWwYcPgw8+GDX7/2zP5Pwx8svh2+X7/ziF8C118qEt8sxeDBw770y6fDs2WhsSxI2Rw7NAL7MzNMBzAfwOBFNB/AUgBXMXAhgRfpvJQvq6oCqKpmxmwn33is3vd//Ply7fGfJEomL33xz1+8tLASuvlrFoSt27ADKy7vnq01NkvdRgsWaODDzXmauSP/eBKAawDgAiwE8n37b8wB0fmmWvPKKPN59d2bv/6M/ktDSr38dnk2+09wM/Pa3ksTPyXCdgXvvlfj4rl3h2uYzxucy9dWFC2UEob4aPE7kHIhoIoDZANYBGM3Me9Mv7QMwupPPPEZE5URU3tjYGImdvrJkiVR2ZLrOz6BBIhC//rVWgnTGmjUyY/fOOzP/jLnhLVsWjk1xYMkSoKQEKCrK7P39+knC/9VXJaemBId1cSCiwQB+BeCLzHy89WvpGXwd3p6Y+WlmnsfM83JzcyOw1E8OHQLef18WK+sOd94pQ/yqqnDs8p2lS4G+fUVEM2XaNKlaev318OzymePHZX5DT3x1/36tsAsaq+JARH0gwvAiMy9JP72fiMakXx8D4IAt++LAW29J7z+TuHhrzPtXrAjepjjw2mvAxz7WvcXfiGSBwxUrNNnfEe++K7Ofu+urZWXyqL4aLDarlQjATwBUM/N3W720FICp/3gQwKtR2xYnli+XG9i113bvc/n5kkRdvjwcu3xm/35g8+aWm1J3uOUW4NgxWcFVuZTly2Veww03dO9zo0cDs2aprwaNzZHDjQA+BWAhEW1I//wJgG8CuJmI6gCUpf9Wesjy5cBNN2WeNG1NWRnw9tu67ENb3n5bHj/+8e5/tqxM9s14440gLYoHy5dLuXX//t3/bFmZhKROnw7erqRis1ppFTMTM89i5mvSP8uY+RAzL2LmQmYuY+bDtmz0nW3bZOnoRYt69vmyMuDkSdkQSGnhrbdkNDZnTvc/e+WVwDXXSAhFaWHvXmDTpp6NxgD53NmzsnaYEgzWE9JKeJgebk/FwfSM9UZ2KW+/DXz0oz0bjQHSO37vPeDcuUDN8hrjqwsX9uzzCxbIgnzvvBOYSYlHxSHGrFkjPdWuZpp2xvDhwPTp2htrzZ49sp5PT0JKhgULZBvRiorAzPKeNWukhPrqq3v2+cGDZUSmvhocKg4xZu1aSe5lszd0aakcR2vIhVWr5PFjH+v5MT7yEXnUGegtrF0rS5v3dDQGiK+uW6eVYEGh4hBTjh6VGG53Kz/acuONcqzq6mDs8p333pOE6TXX9PwYo0dLJZiKg3DyJLBhg9zcs+HGG4FTp4A//CEYu5KOikNMMUnkIBocoMsiG9atk0R0tjvlLVggoxCdgS57bJvd3bLB+Lr6ajCoOMSUtWslnHTdddkdp6AAyM3VWC4gJb0VFRL+yJb582X5jfr67I/lO+ZmPn9+dscZP15moKuvBoOKQ0xZs0aSe223A+0uRNKj096YbIJ05kz2NzGgZVLiBx9kfyzfWbtWiiaGD8/+WKWl6qtBoeIQQy5elPBHtsN0w7XXyrLfx44FczxfMaG6IEYOM2ZI7iLp4sAs1zVIX925E9C1OLNHxSGGpFKyxv28ecEcb+5cefzww2CO5yvr1kkyOdPVbS9Hnz7A7NkqDjt3AgcPBu+rughf9qg4xBBTP9+TGbwdoQ1OWLdORg1EwRzv2mvlu0py6WXQvmqOk3RfDQIVhxhSUSHLSU+fHszxRo2SZF+SG9zx40BtbfcXMLwc8+ZJ6WVNTXDH9I2KCpnZPGtWMMcbOlTKhJPsq0Gh4hBDKiqksWVbbtmauXOT3eA2bpTH2bODO6YmpcVXp02T1ViDIum+GhQqDjGDWRpcUMN0w9y50nM+frzr98aRDRvkMZvJb20pKpIF/MrLgzumb4Tlqzt2SC5D6TkqDjGjoUHq58NocEByk9IbNgAjRwJjxwZ3zF69pNw4qTN69+6Vn7B8VUcP2aHiEDOCTvAZkt7gNmyQUUNQyWjD1VdLyCqJa1eZjkbQvqpJ6WBQcYgZJsE3c2awxx01ChgzpiX2niTOn5e9tIMMKRlmzZKy44aG4I/tOqYjE/R1HToUmDQpmb4aJCoOMaOiomWCVdDMmpXMBrdli2wkE4Y4mCWqkxhaqqhoybsETVJ9NUhUHGLGxo3BlQW2ZeZM2Ts5aXX5YSSjDSUlEqpKojiE7au1tbLcidIzVBxixJEjwO7dwYeUDLNmSQ+6ri6c47vKhg1Av37A1KnBH3vQIFncMGm93JMnga1bw/XVCxd0qflsUHGIEZs2yWNJSTjHNw25sjKc47vKhg3yv2ezEc3lSGLFUnW1lF2rr7qLikOMqKqSx7Aa3LRpkuxOWoMLM/wBiDjU10tiOimE7asFBTLaS5qvBomKQ4yoqpLk3vjx4RzfhFaSFAJpbJSfsG5iQIvwmBtmEqiqkqKJKVPCOX5OjhRmJMlXg0bFIUZUVbUkOMNi5sxk9cY2b5bHoNap6ghTsZSkG1lVlVzT3r3DO0fSfDVoVBxiAnOLOITJrFnAtm3JCYGYPM6MGeGdIz9fNmVKUvI0Kl/du1eX0egpKg4xYf9+4NCh8Ko/DOb4SQmBbN4MXHEFMG5ceOcgknyOGaXEHVNVF7Y4aFI6O1QcYkLYCT5D0sRh0yYJf4QZqgOSJQ7qq36g4hATTO8o7AaXnw8MHJicEMimTeGGlAzTp0tvOglbsUYlDqNHA8OGJcdXg0bFISZUVcn6R7m54Z6nVy+pWEpCgzOVSmEmow3mHEm4rlVVEqrLywv3PCZcl4RrGgYqDjEhigSfISkNzoR5oho5tD5nnImiqs6QFF8NAxWHGMAsN5UoeriANLiGBlkCIc6YSqUoruvEiVL3nwRxiNJXi4ulWOPIkWjOFydUHGLA3r3AiRPSEKLAnGfLlmjOZwtTqRR2+AOQev/i4viLw6FDUloala9OmyaPOnroPioOMcBsUK8NLliiqlQyTJ8ef3EwHQr1VfdRcYgBpsGFsWpoRxQWSk837g3OiENUTJ8u4boTJ6I7Z9RE3ZGZOFGWfYm7r4aBikMMqKmRpZ/DnKjVmr59ZU2cODe4Q4eiq1QymHOZG2gc2bJF/GfixGjO17t3cqrrgkbFIQbU1EhPLKrwBxD/KpDaWnmMajQGJKNiqaamZeQZFXH31bCwKg5E9CwRHSCiqlbPDSeiN4moLv14pU0bfWDLlmhvYoA0uLo62V85jqRS8lhQEN05p0wB+vSJ941sy5boQkqGadOA7duB06ejPa/v2B45PAfg1jbPPQVgBTMXAliR/lvphFOnJE5to8E1N8s+BHGkrk4m/E2aFN05c3KAyZNbRi1x4/x58RcbHRnm+FfXBY1VcWDmdwEcbvP0YgDPp39/HsBdkRrlGeZGYkMcgPj2clMpWSqkX79ozzt1anzFob5eOhTqq35ge+TQEaOZeW/6930ARnf0JiJ6jIjKiai8sbExOuscI+pKJYNp4HFtcKmUxMajpqhIRi0XL0Z/7rCx5atFRTIKjKuvhoWL4vDfMDMD4E5ee5qZ5zHzvNywFxRymJoaSURHfSMbMkSqo+I4VGeWG3SU+QZDURFw9iywY0f05w4bU4UVtTj06yfVUXEdkYWFi+Kwn4jGAED68YBle5ympgaYMAEYMCD6cxcWyk00bhw+DBw9akcczI0zjjeymhpgzBhg6NDoz11UFM9rGiYuisNSAA+mf38QwKsWbXEeG9Ufhrg2OCN4tsJKQDxHZDaq6gzGV7nDOITSEbZLWX8OYC2AqUS0i4geAfBNADcTUR2AsvTfSgdcvGhXHAoLZbLY4bYlBZ5jo4zVMHq0hOziJrrMLfNxbFBYKAtF7ttn5/w+kmPz5Mz8QCcvLYrUEE/ZtUtKWW32xgDpaV9/vR0bwqCuTvI4kydHf26ieFYsNTbKyqg2R7mAXNcxY+zY4BsuhpWUDIl6EbO2tBaHOGGrjNVQVBS/sJKtSiWDCRHGzVfDRMXBY2xVfxgmT5YSwbj1cuvq7OQbDEVFUq0Upxm9tn01P1/WdIqbr/7N3wBLl4ZzbBUHj6mrAwYPBq66ys75zQJqcWtwqZSdfINh6lSJ0cdp9nldnfhLfr6d8/fuLd9pnEYOZ84A3/sesGFDOMdXcfAYM1ErygX32hK3ctbDhyU2bnvkAMQrtJRKydpRUS6415bCwnh1ZLZtk05EWB0ZFQePsd3DBeJXImiEzuZ1NcIUpxuZK76aSgEXLti1IyjCrqpTcfCU5mbpObjQ4E6ckH1644DNMlbDkCHA2LHxEQdmd8Th3Dlg5067dgSFioPSIQ0NIhC2G1zcerk2y1hbE6eKpT17JLmuvhosqRQwbBgwfHg4x1dx8BTTa7AZGwfiV86aSgHjxwP9+9u1I06zz9VXwyHs0ZiKg6e4EP4A4lciaLuM1VBUFJ/Z56746lVXSXVfXHxVxUHpkFQKGDjQXhmroXdvqULRBhcsxoY4lLOmUrLD3fjxdu0wqxfHYeRw7pzsbqfioLTD3MRslrEazB4EvnP4sPy4MHIwjd70un2mrk521MuxuliPEJdwXUODrK2m4qC0w5XwByB2pFL+b1DjSvgDaEmIx0EcbG2c1BGFhVLld+6cbUuyIwpfVXHwkAsXgK1b3biJAS0b1PheIuhK4hSQ/Tny8vwXB1fKWA1FRdKJ2bbNtiXZYcKNKg7KJezcKZu1u9Lg4lIi6EoZq6GgwH9x2L9flspWXw2WVEqS66NGhXcOFQcPcSn8AcQneZpKSW/ddhmroaAgHtcUcGM0BsTLV6dMCTfnqOLgITZ3KuuIsWPlhup7L9elPA4gN7L9+4GmJtuW9BwXliNpzYgRsk1pHMQh7Guq4uAhqZTEpF3ZtKRXLwnFaIMLlilT5NHn65pKSZXShAm2LRGI5Lr63JGJKueo4uAhZkjZy6Fvz/f4+JEjMunMtZED4Pd1TaVkWXcXylgNvvtqVDlHh24vSqa41sMFRKzq6/1dndW1PA7QMnLw+UbmUhmrYcoUmUDW3Gzbkp4Rla+qOHjGxYtyE3atwRUUyOJqe/fatqRnuJbHAWR11tGj/RUHZrmuLgkuIPY0N8tuez6i4qB0yK5dMqfAtQbne3zcNDhXylgNPodAGhslma6+Giz19VIAMnZsuOdRcfAMF8MfgP/x8bo6KWMdMMC2JZficzmr+mo4RJVzVHHwDFcbXH6+LMLn843MpZCSoaBARounT9u2pPu4NsfBMGaM9Lx99tUo2r+Kg2ekUkC/ftLLdYk+faQqxefemGuCC7SEQLZutWtHT6irkw6DK2Wshl69/C1nNTlH4xdhouLgGXV17pWxGnwNgRw9Chw86F4PF/A7BJJKiTD07WvbkvaY6jrf2Ls3ul31HLzFKJfD1R4uIA2urs6/clZXQ3WA/+Lg4jUFWjoy6qudo+LgEWZI6XKDO3bMv93LXCxjNVx5pewR7Js4mDJWF68pIB0ZH0uvVRyUDjFDSpcbHODfcN3VMlaDj+Wshw5JR8Hljgzg33WNcle9jMSBiL6VyXNKuLi2iFlbfG1wpox14EDblnSMj7kcl0N1gL++Wl8f3a56mY4cbu7guT8O0hCla1xvcJMmyaOPNzJXrykgtjU0+LV7meu+mp8vN1j11c65rDgQ0WeJqBLAVCLa2OpnG4CN0ZioGFIpqfywvVF7Z/i6e5mLSzy0pqBA8k3bt9u2JHNSKamoMx0G18jJ8a/0Oupd9boanPwMwGsAvgHgqVbPNzGzZ2lH/0mlJC7eu7dtSzrHtxJBl8tYDa0X4CsqsmtLptTVSe+8Xz/blnSOb75qliOJYo4D0HVYqTeA4wAeB9DU6gdENDxc05S2mGnzLuNb8jSKvXizxcf4uE++6ks5a9Shuq7EYT2A8vRjI4BaAHXp39eHa5rSGtc2au+MKVNk97ITJ2xbkhkul7EacnNlhVafxMHFlYPbMmWKX6XXTokDM09i5skAlgO4g5lHMvMIALcDeCMKAxXBtY3aO8O3PXpNg3O5l0vk14js6FEpZXX5mgL+jchMHmfixGjOl2m10nxmXmb+YObXAJSGY5JARLcS0RYiShHRU11/It6Ym60vDc4XcairA8aNc7eM1eBTOasPoTrAv3k5US9Hkqk47CGifyCiiemfvwewJyyjiKg3gB9BymWnA3iAiKaHdT4fcL000ODb7mU+hOoAsXHrVj92L/NhNAZIcQeR+mpnZCoODwDIBfAKgCXp3x8IyygA1wFIMfNWZj4H4CUAi4M+yapVwOLFwIEDQR85eOrr3Vzhsi1XXCExcl96Yy4v8dAan3Yvc33GuaF/fxk1+uKrUS+d09U8h68S0WxmPszMTzDzbGaew8xfDLmUdRyAna3+3pV+rrVtjxFRORGVNzY29ugkJ04AS5e2JCVdJpWS0kAXV7hsiy/LIR87JuWBvowcAD+ua3297JkwaJBtS7rGl1zO4cPy44w4ANgK4Aki+pCIniOiPyWiK6MwrCuY+WlmnsfM83Jzc3t0DJ8anC/hD8CfBudLqA5QXw0LX+Y62Mg5dlWt9Atm/itmng3gewAmA1hCRO8S0T8R0XUh2bUbQOt5wHnp5wJlwgQJ1fjQ4KLa4CMIpkwBdu6Uva5dxtWdyjpizBiZge6Lr/oiDgUFUgnY1GTbkstjoyOT8aqszPwhM3+DmW+ClLJuAvBoSHZ9AKCQiCYRUV8A9wNYGvRJfNm97MiR6IeU2VBQIPMytm2zbcnlMeFEH0TXl3LWkyeBPXv8uKaAP9V1NvI4ma7K+r+JqO1SG7cz82Mh2ARmbgbweQCvA6gG8Etm3hTGuXxocL6UBhp8KRFMpYCxY/2IjQN++KrZzlR9NVhSKVm3bMCA6M6Z6cghB8A6IppFRDdDevahzpBm5mXMXMTMU5j5a2Gdp6DA/d3LfCkNNPgSH/cpNg60zHW4cMG2JZ3jy3wcgy+l1zZ8NSNxYOavAngSwDoAzwO4jZl/GKZhUWF2Lzt0yLYlneNLaaBh5EhZ7sH13pgvZayGggJZtnt34Nm34PCtI+NL6bWz4kBEHwXwAwD/DOAdAD8gorFhGhYVPvRy6+sl/OH6LF6DD/Hx48dlfotvIwfA7etaXw+MGCHbm/qC677a1GTHVzMNK30HwD3phPQDAJ4BsDI8s6LDhwbnW/gDcL/B+VSpZPDFV30ZNRhcn5djK+fY1SS4LxHRlwD8AsAft/p7IoCnI7AvdCZNcn8KvU+lgYaCAtmcxtXlHnya42DIy5P9EVz2VV87Mrt2AWfO2LakY2z5alcjhyHpn7kAPguZpTwOwGcAzAnXtGjo109mHrva4E6eBPbu9a83VlAAnD/v7nIPPpWxGnr1kryTq7567px83z5dU0BGj8wtlVauYSvneNmd4Jj5fwEAEb0LYA4zm41+/ieA34ZuXUS4HALxrYzV0DoE4mIiPZWSiWWDB9u2pHu47Kvbt8t2pj776nQHl/dMpYDRo6XII0oyzTmMBtB6e/Nz6edigcsNzrfSQIOJ5bu6bpXr+0Z3hsu7l/nekXHVV22F6jIVhxcAvE9E/zM9algH4LmwjIqaggIpZT1yxLYl7fGtNNBw1VUyucxV0U2l/EpGGwoKgNOnJdToGr766vDh8uOyrzorDulJaA8BOJL+eYiZvxGmYVFibhIu1jrX18u8gWHDbFvSPUw5q4u9saYmWU/Htx4u4HbFUiolYbpRo2xb0n1cjR6cPi3zWpwVBwBg5gpm/l7658MwjYoa1xucbz0xg6sNzscyVoPLvmoWhySybUn3KSx0syNjczmSjMUhzpiEqYsNzsfSQIPZvcy15R58LGM15OcDOTnqq0FTUCCVVq6tJGzTV1UcIItZ5eW51+DOnpWlr30dORQWulnOanqIPt7IcnJkbo5rvnrhgqzC6+M1BdwtZ7WZx1FxSONifNzX0kCDqyGQVEoS5r6VsRpcDNft2iXzHHztyLjsq8OH21mORMUhjYsNztcyVoOr5ay+lrEaXCxn9TlUB7hbzmozVKfikKagQBa3On7ctiUt+N7gXN29zOfYOCC2NzXJ/teu4HtHxiwWqL7agopDGhd3hEqlZFZkD7fIto6L5axNTcC+fUBRkW1Leo6LIZBUSpaiycuzbUnPcS16cPas5OtUHCzjYoPzuTTQUFjo1jU1QuVjGavBVV+dPFnWf/IV18pZt26VnKMtX/X4qwwWF3eE8j38AbhXzhoHcZg4UW7CrvmqryElg2vlrLZ9VcUhzeDBUsHiSoMzpYG+N7jCQqli2bnTtiWCz2Wshr59gQkT3PFVZj+XlW9LQYH01Ldts22JoOLgEC7FHHfulDkCvouDayGQujrZVW/QINuWZIdLvrp/vywt77uvmpuwK9e1rk4S5cOH2zm/ikMrXGpwptfgc+IUcK9EsK7O/2sKuOmrPofqAPd8tbbW7jVVcWhFQQGwZ4/0gmxTWyuPvt/Ixo51q5y1rs7/mxggvnrkCHD4sG1L4uOrI0bIApfqq4KKQytMz8GFKfS1tS15EJ/p1UvCDS70xo4cAQ4ejI84AG7cyGprJQ+Sn2/bkuxwqfT61CmZda7i4AguNTjTa/C5jNXgSjlrXMIfgFu+WlsrHYDevW1bkj2u+KqZb2VzNKbi0ArXGpzvw3RDQYE4u+1y1jiJw+TJ0nFwwVfjkscBxFcbGqTCziYu+KqKQyuGDpXZyLYb3LlzUk4XlwZnyll37bJrR12d3FB9r6oBgP793VhJ+MIFsSEuvupKOauKg4O4EHO0PTMyaFwZkdXVSVy8f3+7dgSFCxVLO3fKpLG4iIMr5ay1tcDo0bJ8ji1UHNrgQoOLSxmrwZUSQdvVH0Hjgq/GpVLJoL7agopDGwoKpDd0+rQ9G0yDs+0cQTFunPTWbd7ImN1ocEFSUCArsx47Zs8GF8IfQTJypISXbYuuC3kcFYc2uFDOWlsrTmprZmTQuFDOeugQcPRofG5igBvhuriUXBtcKGc1Kwfb9lUVhzZMnSqPW7bYsyFOlUoG2yWCcRuNAS0+4oKvxqHk2mA7XGfObdtXVRza4EKDi1v4A2gRB1vlrHELfwByE+vVyw1xiBNFRbJFr63VWV3pyKg4tGHIEFnyoabGzvlPnAB2745fgysulnLW7dvtnL+uTiZpTZpk5/xh0L+/LN9tSxzM92n7JhY0U6dKtaCtjb9cWTlYxaEDiovtNTgzpIybONgO19XVyY20b1875w+L4mJ7HRlTch03Xy0ulkdb17WuTuawDBxo5/wGFYcOmDpVbmI2NnCPW2mgwXaDs73CZVhMnSr/28WL0Z87rr5qOjI2xcEFX7UiDkR0LxFtIqKLRDSvzWtfJaIUEW0holts2FdcLJUtBw5Ef25XhpRBM2KEVGDZaHAXL4rYG4GKE8XFUnZtYzOlOOZxAKm+ysuz46vMcl4jUDaxNXKoAvAJAO+2fpKIpgO4H8AMALcC+A8iinw5L5shkNpaN4aUYWArBGLmrUybFv25w8a2r44cCVx5ZfTnDhtbvtrYKKsHu+CrVsSBmauZuSN3XgzgJWY+y8zbAMkyOnwAABTASURBVKQAXBetdXZDIHGs/jCYcF3UmO8xriMHQH01aGyFlqur5dEFX3Ut5zAOQOsB8q70c+0goseIqJyIyhsbGwM1Yvx42aAm6huZS0PKMCgullBd1BvUuNTggmbUKJnRa0t04+yrx4/LZLQocakjE5o4ENFyIqrq4GdxEMdn5qeZeR4zz8vNzQ3ikP9Nr17SI4q6N7Z/v+Q6XBhShoFx+KhvZDU1Mts8YDdxAiI7IZCjR+XGGXdfjfq61tTI/uZ5edGetyNCEwdmLmPmkg5+Xr3Mx3YDGN/q77z0c5FjIwRierja4IKlpkbOHadZvK1RXw0eW75aXS3fZy8HYjoOmHAJSwHcT0T9iGgSgEIA79swZOpUWdM9ylmScW9wZp6BLXGIK8XFMnGyqSm6c27eLI9x9dVx46QHb2OU68o1tVXKejcR7QJwA4DfEtHrAMDMmwD8EsBmAL8D8DgzW1lwobhYSiCjXGNl8+aWGdpxJCdHSnSjbHBHjki4Ls7iYOL+Zt5BFFRXA/36ieDHESK5rlF2ZE6dkl3oXPFVW9VKrzBzHjP3Y+bRzHxLq9e+xsxTmHkqM79mwz7ATolgdbX0GuIa/gCij4+bc7nSGwsDGyEQE/6Iw77RnRG1r5p7TaLFwQdszJI04hBniotlzZrz56M5n0vVH2ExZUr0C/AlxVcbGqRHHwWudWRUHDph8GCJO0bV4I4dA/budccxwqK4GGhujm5Rs5oayXPENfwBSHhn8uToOjKnT8uCe0nwVSC6cF1NjYi8K6sjqDhchuLiliRx2MQ9GW2Iupy1ulrKknNyojmfLaKMj5vJYdOnR3M+W0QdWq6uFpHv1y+a83WFisNlmDFDksRRLGqWFHEwDS4q0Y17pZJhxgy5iTU3h3+uuFcqGQoLJf8Xpa+6dE1VHC5DSQlw8qTEHcNm82bpMcRpv4GOuOIKmeCzaVP45zp7VpaVToI4lJTI/gpRVNdVV0v4I24L7rVlwADJ50ThqxcuSPjKJV9VcbgMJSXyWFUV/rmSEv4AgJkzgcrK8M+zZYs0uhkzwj+XbaL21SlT3Al/hMnMmdFc0/p66czoyMETzE0lqgbnkmOESUmJ/L9hh0CMAM2cGe55XKC4WHrz6qvBUlIiS5OfORPuecz35pKvqjhchiuuAPLzw29wp07JbOwkNbgoQiBVVUCfPvFdObQ1AwZIlUvYvnr+vNwsk+SrFy6En+yvrJT8hktJfhWHLigpCT/muHmzVH+41GsIE/N/hn0jq6yUHnWfPuGexxVKSsK/plu2iEDMmhXueVwhqnBdZaWIu0v7uKg4dEEUIZCNG+UxKQ3OhEDCzjtUVrY07iQQRQjE+GpSOjKFhTJPJgpxcM1XVRy6IIoQyMaN0mOYPDm8c7hEFCGQ48eBHTuScxMDxFcvXgw3BLJxo4zE4rqPQ1v69JHOTJgdmdOn5f7imq+qOHRBFEnpyko5T5zXqWlL2FUg5tiu9cbCJIoQSGWl5Bv69g3vHK4RdriuulpEXcXBM8xCeGE5BzPwhz8kJ6RkKCmR3tLp0+EcP4niUFAgN+0we7kbN7p3EwubmTNlFHr8eDjHN9+Xa76q4tAFYYdA9u0DDh1KpjhcvBje7NPKSlkfa8KEcI7vIiYEEpavHj4M7NqVTF8FwruulZUyZ8SVNZUMKg4ZUFISXm8sSbX4rQm7wVVVSajOhR21oiQKX1VxCJaqKilhdW0CbMKaTs+4+mqpAjl5MvhjJ636w1BQIL2lMG5kzG5Wf0TBzJnAzp2yyVHQJLUjM2GCbMIVpui66KsqDhkwe3ZLbiBoNm6Und9Gjgz+2C6TkyM3mQ8/DP7Yu3ZJqG727OCP7Tpz5sjjhg3BH3vjRmD48PjuVNgZRDJaCsNXDxwA9uxxczSm4pABpsFVVAR/7CQm+Axz5sg1ZQ72uOZ7Mt9bkjCCGKavxnmnws6YO1cE90LAmxYbwZk7N9jjBoGKQwaMGwfk5gbf4M6ckdnX11wT7HF9Yc4cCX9s3x7scSsqJNfgYm8sbHJzgfHjg/fV5mYRhySOxgDx1ZMnJbwcJOZ7cvG6qjhkAJE4R9DDyspKaXTz5gV7XF8Ia0RWUSFVO4MGBXtcX5g9O/hrWl0tZcdJ99X164M9bkWFTH4dNizY4waBikOGzJkjVQVnzwZ3zPJyeUxqg5s5U3IPYYhDEkNKhjlzZA2kEyeCO6bxVRfDH1EwbRrQv384vurqNVVxyJA5c6SXH2Q52/r1kuBLUi1+a/r3l3LTIHtj+/ZJgs/VBhcFc+YEX0Cxfr3MG0nCCrcdkZMjVYtBisORI7IZlasdGRWHDAkj0bd+vYwakpjgMwSdlDahP1cbXBSY/z3IMOj69XLcpM0baY3x1aC2DXbdVxP8VXePyZOBoUODE4czZ2QUktSQkmHuXKCxEdi9O5jjme8nqUl+QEpNR40Kzlebm6VSR31VltDYujWY47mcjAZUHDLGJKU/+CCY423cKI0uyeEPIPhE3wcfyDLLV1wRzPF8xPiqyRNky+bN0plRX5XHoHx1/XqpLMvNDeZ4QaPi0A3mz5c47qlT2R8r6clow9VXSzx33brsj8UMrF0L3HBD9sfynfnzZWQaxGJx6qvCjBmSJwvCVwHx1fnzgzlWGKg4dIPSUuntB9Eje+89GfqPH5/9sXxm4EAZVq9Zk/2xtm2TGacqDuKrzMHcyN57T0otXVsYLmr69gWuvTYYX92zB2hocNtXVRy6gVH5tWuzP9bq1cCNNyY7GW0oLQXef1+2n8wG87243OCi4vrrxbeCuJGtXi3XNMnJaENpqeQKsl1q3gdf1a+7G4wcKfHsbBvcvn2S1LrxxmDs8p3SUmls2ZZerl0r5ZYuLmIWNVdcIdch247M4cOSc1BfFUpLpROTbd5h7VpZeNLVZDSg4tBtSkvli82m9HL1annUBieUlsqjuS49Ze1a4LrrkrWj3uUwvppN6aURF/VVwfhqth3ENWskwd+vX/Y2hYWKQze54QYpvcymnG31aklsuVrfHDV5eUB+fnYN7uRJGXm4PEyPmtJSSUhv3tzzY6xeLQUD110XnF0+M3KkTATMpiNz9qyMPFz3VRWHbmJ6DqtW9fwYq1dLYitJ+/B2RWlpdg1uzRpZMfMjHwnOJt8xN59sruvq1RL6GDgwGJviQGmp+FtPowfvvw+cO+f+aEzFoZvMmCG9hxUrevb5piZJaOlN7FIWLJCJcPX1Pfv8ypXSw12wIFi7fKagQFYUXrmyZ58/c0ZuZOqrl7JgAXDwYM9HZCtXSrHAxz4WrF1Bo+LQTXr1AhYuFHHoSc/hnXekHHbRouBt85myMnl8882efX7lSqkmS+pKrB1BJH62YkXP8g6rVolAqK9eSra+umKFhJSHDw/OpjBQcegBZWVSp1xT0/3Pvvmm5BtcH1JGTWGh5B160uCOHZO5JwsXBm+X75SVya54PakEe/NNoE8f93u4UZOfL3mH5cu7/9mTJ2XeiA+Ca0UciOjbRFRDRBuJ6BUiGtbqta8SUYqIthDRLTbs6wrzxfYktLR8OfDRj4pAKC0QATffLCOA7u629e670jNWcWhPtr56ww1SHqxcSlkZ8PbbkjvoDqtWSSmsD75qa+TwJoASZp4FoBbAVwGAiKYDuB/ADAC3AvgPInKuMHHyZGDSpO73HHbvljilGZYql1JWBhw92v0a8tdek4Tp9deHY5fPjB0LTJ/e/RHZwYOyaqj6asfcfHPLKKA7vP66FKL4kMexIg7M/AYzN6f/fA9AXvr3xQBeYuazzLwNQAqAk0V0t94qDa47MyWXLZPHW5wcD9ln0SIZQbz2WuafYQaWLpVrqqOxjrnlFunlNjVl/pnXXpNrq77aMTfdJPNpeuKrixb5kRtzIefwMABziccB2NnqtV3p59pBRI8RUTkRlTc2NoZsYnvuvlsW4Hvjjcw/88orMuKYOTM8u3wmN1dyMUuWZP6ZDz+UEdmdd4Znl+/cfbeEP7pzI3vlFRl1JH2xvc4YOlQEYsmSzAtTNm+WarzFi8O1LShCEwciWk5EVR38LG71nr8H0Azgxe4en5mfZuZ5zDwv18Katx//uCxG9sormb3/2DEJQ33iE7qe0uW45x5ZzjzTjdyXLpUKsttuC9cunyktFeHNVHRPngR+9zsRFV1PqXPuuQeorQU2bcrs/UuXyuMdd4RnU5CE9tUzcxkzl3Tw8yoAENFfAbgdwJ8z/7f27gbQep3SvPRzztGnj3zJS5dmtmDcb38r77v77vBt85lPfEIef/Wrrt/LDLz0ksRvXV0T3wV69wbuukt88MyZrt//+usSLlVfvTx33SUdvUx8FQB++UuZaT52bLh2BYWtaqVbATwJ4E5mbr07wlIA9xNRPyKaBKAQwPs2bMyE++6TfWBNLuFyPP+8LM/t8vrtLpCfLw3oZz/rerj+wQfAli3Apz4VjW0+c999wIkTwKuvdv3e558HRo+Wqjqlc666SibEZeKrGzfKbno++aqtQeMPAQwB8CYRbSCi/wsAzLwJwC8BbAbwOwCPM3M3Cxuj49ZbpRfwzDOXf19DgySvH3pIF4XLhEceASoru96L4IUXJAl9773R2OUzCxcCEyZ07at798oI48EHZXSsXJ6HH5bQ0jvvXP59zz8v1/OBB6KxKwhsVSsVMPN4Zr4m/fOZVq99jZmnMPNUZu5GCi16cnLkhv/aa8DOnZ2/7yc/kceHHorGLt954AGprf/xjzt/z7FjIg733CPJQeXy9OolortiBZBKdf6+n/5U5pk8/HB0tvnMffdJ7vFyvnryJPDcc1I0MWJEZKZljaabsuTTn5aG96//2vHrTU3AD38ojjFxYqSmecuQITL8/tnPOhfdp5+Wa/ulL0Vrm8888ojU2H/rWx2/fuoU8O//LiPiqVOjtc1XBgyQTt/LL3e+LthPfiL7Ynz5y9HaljXM7P3P3Llz2SaPPsrcty9zQ0P71/75n5kB5vffj94un2lokGv66KPtXzt8mHnECOaysujt8p3Pf545J4c5lWr/2re+Jb76+99Hb5fP7NnD3L8/81/8RfvXjh1jvuoq5o98JHq7MgFAOXdyX7V+Yw/ix7Y4bN/OPHAg8623Ml+82PJ8TY04zSc/ac82n3niCWYi5nfeufT5xx5j7tWL+cMP7djlM7t3Mw8ZwnzTTcwXLrQ8X18vPnz77fZs85m//Vu5m77xxqXPf+EL8vy6dXbs6goVhwj4wQ/kan7xi8znzzNv28ZcVCQ93F27bFvnJ01NzAUFzKNHM2/YIML77W/LdX7ySdvW+cszz8g1/NznmM+dY96xg3n6dOZhw6Sjo3SfkyeZp01jHjmSubxcfPX735fr/Nd/bdu6zlFxiICLF1t6CcOHM/fpwzx0qA7Rs2XzZuYxY2SkMGqUXN977hEBVnrGxYvMX/mKXMsrr5Tw3ZAhzCtX2rbMb2prmfPyZLRrfPWOO5jPnrVtWedcThxIXvebefPmcXl5uW0zwAz85jcyazo3F/jsZzUJHQSNjZLUb2iQdWn+/M915m4QLFsmE7iGDwc+8xlgyhTbFvnP4cPAD34g2wjfdJMUVrhcvk5E65m5w0VSVBwURVESyuXEQftfiqIoSjtUHBRFUZR2qDgoiqIo7VBxUBRFUdqh4qAoiqK0Q8VBURRFaYeKg6IoitIOFQdFURSlHbGYBEdEjQAaevjxkQAOBmiOD+j/nAz0f04G2fzPE5i5w012YyEO2UBE5Z3NEIwr+j8nA/2fk0FY/7OGlRRFUZR2qDgoiqIo7VBxAJ62bYAF9H9OBvo/J4NQ/ufE5xwURVGU9ujIQVEURWmHioOiKIrSjkSLAxHdSkRbiChFRE/ZtidsiGg8Eb1FRJuJaBMRPWHbpiggot5E9CER/ca2LVFBRMOI6GUiqiGiaiK6wbZNYUJEf5P26Soi+jkR9bdtUxgQ0bNEdICIqlo9N5yI3iSiuvTjlUGcK7HiQES9AfwIwB8DmA7gASKabteq0GkG8GVmng5gPoDHE/A/A8ATAKptGxEx3wPwO2YuBnA1Yvz/E9E4AH8NYB4zlwDoDeB+u1aFxnMAbm3z3FMAVjBzIYAV6b+zJrHiAOA6AClm3srM5wC8BGCxZZtChZn3MnNF+vcmyA1jnF2rwoWI8gDcBuA/bdsSFUQ0FMBHAfwEAJj5HDMftWtV6OQAGEBEOQAGAthj2Z5QYOZ3ARxu8/RiAM+nf38ewF1BnCvJ4jAOwM5Wf+9CzG+UrSGiiQBmA1hn15LQ+XcATwK4aNuQCJkEoBHAT9PhtP8kokG2jQoLZt4N4DsAdgDYC+AYM79h16pIGc3Me9O/7wMwOoiDJlkcEgsRDQbwKwBfZObjtu0JCyK6HcABZl5v25aIyQEwB8D/YebZAE4ioFCDi6Rj7IshojgWwCAi+gu7VtmBZW5CIPMTkiwOuwGMb/V3Xvq5WENEfSDC8CIzL7FtT8jcCOBOItoOCRsuJKL/Z9ekSNgFYBczm1HhyxCxiCtlALYxcyMznwewBECpZZuiZD8RjQGA9OOBIA6aZHH4AEAhEU0ior6QBNZSyzaFChERJA5dzczftW1P2DDzV5k5j5knQr7flcwc+x4lM+8DsJOIpqafWgRgs0WTwmYHgPlENDDt44sQ4wR8BywF8GD69wcBvBrEQXOCOIiPMHMzEX0ewOuQ6oZnmXmTZbPC5kYAnwJQSUQb0s/9HTMvs2iTEg5fAPBiuuOzFcBDlu0JDWZeR0QvA6iAVOR9iJguo0FEPwfwcQAjiWgXgP8B4JsAfklEj0C2LrgvkHPp8hmKoihKW5IcVlIURVE6QcVBURRFaYeKg6IoitIOFQdFURSlHSoOiqIoSjtUHBQlJNIro37Oth2K0hNUHBQlPIYBUHFQvETFQVHC45sAphDRBiL6tm1jFKU76CQ4RQmJ9Mq3v0nvMaAoXqEjB0VRFKUdKg6KoihKO1QcFCU8mgAMsW2EovQEFQdFCQlmPgRgdXrTe01IK16hCWlFURSlHTpyUBRFUdqh4qAoiqK0Q8VBURRFaYeKg6IoitIOFQdFURSlHSoOiqIoSjtUHBRFUZR2/H8INgSzDlVovgAAAABJRU5ErkJggg==\n" }, "metadata": { "needs_background": "light" } } ], "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 541 }, "id": "vtf1Dwi2VSR1", "outputId": "2b7c39ac-8532-4f72-dfe1-d86871330fde" }, "execution_count": 3 }, { "cell_type": "markdown", "source": [ "It is important to check that the time step used was small enough to get accurate results. If $dt$ is small enough, changing it slightly won't change the results." ], "metadata": { "id": "iQV_e8alVSR2" } }, { "cell_type": "markdown", "source": [ "Note that this method can be used for ordinary, second-order differential equations, even if the function depends of something other than time. For example, it could be used to solve \n", "\n", "$$\n", "\\frac{d^2z}{dx^2} = 3z - \\frac{dz}{dx}\n", "$$ \n", "\n", "for $z(x)$." ], "metadata": { "id": "a0kGcPSlVSR2" } }, { "cell_type": "markdown", "source": [ "##Reference" ], "metadata": { "id": "8vdkesv4VSR2" } }, { "cell_type": "markdown", "source": [ "[Alan Cromer, \"Stable solutions using the Euler approximation,\" *Am. J. Phys.* **49**, 455-459 (1981).]( http://scitation.aip.org/content/aapt/journal/ajp/49/5/10.1119/1.12478) \n", "(Note: As mentioned in the article, this method was discovered accidentally by a high school student.)" ], "metadata": { "id": "A0B7kGCLVSR2" } } ] }