{ "metadata": { "name": "", "signature": "sha256:0fade634d3cf2129a907c7c3366975fb52c3afd1b28c9df07d6dcb2e2bd1c30d" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# The Explicit Euler Method and Order of Accuracy\n", "\n", "By Thomas P Ogden (<t.p.ogden@durham.ac.uk>)\n", "\n", "_Here we introduce numerical integration of ODEs using finite difference (FD) methods with the Explicit Euler method, use it to solve a first-order ODE and compare with the analytic known solution to check the order of accuracy_." ] }, { "cell_type": "code", "collapsed": false, "input": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 6 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "Say we have an ordinary differential equation (ODE) $y' = f(t,y(t))$ with an initial condition $y(t_0) = y_0$ and we want to solve it numerically. If we know $y(t)$ at a time $t_n$ and want to know what $y$ is at a later time $t_{n+1}$, the fundametal theorem of calculus tells us that we find it by integrating $y'$ over the time interval,\n", "\n", "$$ y(t_{n+1}) = y(t_n) + \\int_{t_n}^{t_{n+1}} \\! y'(t) \\, \\mathrm{d}t = y(t_n) + \\int_{t_n}^{t_{n+1}} \\! f(t,y(t)) \\, \\mathrm{d}t.$$\n", "\n", "The idea behind any ODE solver is to compute the right-hand-side integral for some numerical approximation of $f$. The problem is then computed over a series of steps $n = 1, 2, \\dots N$ to give a sequence of points $y_n$ which approximate $y(t)$ to some order of accuracy as a function of the stepsize. The method is **consistent** if the local error (i.e. the error from step $n$ to step $n+1$) goes to zero faster than the stepsize ($t_{n+1} - t_n$) goes to zero." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The Explicit Euler Method\n", "\n", "The first numerical approximations to $f$ we're going to look at are based on Taylor series expansions. We know we can expand any infinitely differentiable function $f$ in a Taylor series around the point $t_n$ as \n", "\n", "$$ y(t_{n+1}) = y(t_{n}) + \\frac{y'(t_n)}{1!} h + \\frac{y''(t_n)}{2!} h^2 + \\frac{y'''(t_n)}{3!} h^3 + \\dots $$\n", "\n", "where $h := t_{n+1} - t_{n}$ is the distance between the points. We'll focus on methods where $h$ is the same for each step (i.e. $\\forall n$).\n", "\n", "The most basic approximation then, called the **Explicit Euler** method, is to truncate the series at the linear term, discarding quadratic and higher terms. Putting in our ODE $y' = f(t,y(t))$ gives us the finite difference step\n", "\n", "$$\n", "y(t_{n+1}) \\approx y(t_n) + hf(t_n, y_n) + \\mathcal{O}(h^2)\n", "$$\n", "\n", "and so our sequence of approximation points $y_n$ is calculated as\n", "\n", "$$\n", "y_{n+1} = y_n + hf(t_n, y_n).\n", "$$\n", "\n", "As we decided to truncate the Taylor expansion after the linear term, the error at each step, called the **local error**, will be on the order of $\\mathcal{O}(h^2)$ if $h$ is small and the cumulative effect of these errors over many steps, called the **global error**, will be on the order of $\\mathcal{O}(h^1)$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Implementing the Method in Python\n", "\n", "We'll now define a function to implement the Explicit Euler method for a first-order ODE given an initial condition. Below we'll try it out on a test problem." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def ode_int_ee(func, y_0, t, args={}):\n", " \"\"\" Explicit Euler approximation to an first-order ODE system with initial\n", " conditions.\n", "\n", " Args:\n", " func: (callable) The first-order ODE system to be approximated.\n", " y_0: (array) The initial condition.\n", " t: (array) A sequence of time points for which to solve for y.\n", " args: (dict) Extra arguments to pass to function.\n", "\n", " Out:\n", " y: (array) the approximated solution of the system at each time in t,\n", " with the initial value y_0 in the first row.\n", " \"\"\"\n", "\n", " y = np.zeros([len(t), len(y_0)]) # Initialise the approximation array\n", " y[0] = y_0\n", "\n", " for i, t_i in enumerate(t[:-1]): # Loop through the time steps\n", "\n", " h = t[i+1] - t_i # size of the step\n", " y[i+1] = y[i] + h*func(t_i, y[i], args) # Euler step\n", "\n", " return y" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example 1: An Exponential\n", "\n", "To see that our Explicit Euler solver actually works, let's take a simple ODE: \n", "\n", "$$y'(t) = ay(t)$$\n", "\n", "with intial value $y(0) = 1$ and $a \\in \\mathbb{C}$. We know the analytic solution of this equation is $y = \\mathrm{e}^{at}$ so we can check the accuracy of the method against this.\n", "\n", "Here's that ODE written as a Python function." ] }, { "cell_type": "code", "collapsed": false, "input": [ "def exp(t, y, args):\n", " \"\"\" An exponential function described as a first-order ODE. \"\"\"\n", " \n", " dydt = args['a']*y\n", " return dydt" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "I've gone for a one-line docstring in this case as the ODE function is short. We would document the arguments and return variable for a more involved system. It might seem odd to have $a$ passed as an item of the dict `args`, but for more involved problems we'll want to pass more arguments, and this is a tidy way. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Solving the Problem\n", "\n", "Let's set our initial condition $y_0$ and make our list of timesteps $t_1, t_2 \\dots t_N$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "y_0 = np.array([1.]) # Initial condition\n", "\n", "N = 10 # Num of steps to take\n", "t_max = 5.\n", "t = np.linspace(0., t_max, N+1) # Array of time steps" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "And we need our factor $a$, which we'll just set to $1$ for now." ] }, { "cell_type": "code", "collapsed": false, "input": [ "solve_args = {'a':1}" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we're ready to solve using the Explicit Eulder integrator we wrote, and we'll plot the results alongside the known analytic result $y = \\mathrm{e}^{t}$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "y = ode_int_ee(exp, y_0, t, solve_args) # Solve the ODE with Explicit Euler\n", "\n", "plt.plot(t, y[:,0], 'b-o', label='Explicit Euler')\n", "plt.plot(t, np.exp(solve_args['a']*t), 'k--', label='Known')\n", "plt.xlabel(r'$t$')\n", "plt.ylabel(r'$y(t)$')\n", "plt.legend(loc=2)\n", "\n", "plt.show()" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAAEiCAYAAABgAHyVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlcVOX+B/DPsCibC7iALIWBEmLEIoqoMCqKS2K0XVwu\nGkreuil667a4FFou91qRYpmZYGoi15YfKkiKipiGmGnuCSg5jJIKuLAzw/z+mBhFBmVg4MwMn/fr\nNS+ZM+ec+Z5x5OPznOc8R6RQKBQgIiIyYEZCF0BERNTaGHZERGTwGHZERGTwGHZERGTwGHZERGTw\nGHZERGTwGHZERGTwGHZERGTwBAu7yMhI2Nra4qmnnqq3vKysDNOmTYO3tzf69euHo0ePAgAyMzPh\n4+MDT09PxMXFCVEyERHpKZFQM6gcOnQIVlZWiIiIwOnTp1XLp02bhqCgIERGRkImk6GsrAxWVlZw\nc3NDeno6HBwc4Ofnh8TERLi7uwtROhER6RnBWnbDhg2DtbV1vWW3b9/GoUOHEBkZCQAwMTFBly5d\nkJ2dDVdXVzg7O8PU1BTh4eFITk4WomwiItJDOnXO7vLly+jRowemT5+O/v37IyoqChUVFZBKpXBy\nclKt5+joCKlUKmClRESkT0yELuB+MpkMx44dw8KFC7F27VrMmjUL27dvh6WlpUb7cXV1RV5eXitV\nSUREusLFxQW5ubmPXE+nWnaOjo7o1q0bJkyYAHNzc0yaNAm7d++Gg4MDJBKJaj2JRAJHR8dG95OX\nlweFQsFHEx/vv/++4DXo04OfFz8vfl6682hqw0anws7Ozg6urq44evQoamtrkZKSguDgYPj5+SEn\nJwf5+fmorq5GUlISQkNDhS6XiIj0hGBhN2nSJAQEBODixYtwcnJCQkICAODrr79GdHQ0+vbtC6lU\nivDwcBgbGyM+Ph5hYWHw9fVFZGQkR2ISEVGTCXbpQWsSiUQwwMNqNRkZGRCLxUKXoTf4eWmGn5dm\n+Hlppqm/7xl2RESkt5r6+16nztkRERG1Bp269KC12djYoKSkROgySEdYW1ujuLhY6DKIqA20q25M\ndm/S/fh9INJ/7MYkIiL6C8OOiIgMHsOOiIgMHsOunZo+fTqioqJUz/v374/t27c3adtOnTqp7jPY\nVmJiYjBq1Kg2fU8iMhztajSmPhGLxcjKyoKpqWm95VlZWfDw8Gjx/kUiEUQiker5mTNnmrzt3bt3\nVT9nZGRg1KhRqKmpeeg2zs7O+PPPP2Ficu8rJxKJIJVK0alTJw0qJyLSHMPuLykpmVi9eg+qqkzQ\nsaMMc+aMxvjxgYLtRyQS4b333sP8+fM13rap2nIkokgkwoYNGzB58uQ2e8/71dTUNPiPAxG1H+zG\nhDKgoqN/xJ49H+LgwRjs2fMhoqN/REpKpiD7eZTS0lL069cPS5cuVS374IMP0K9fP1RUVAAAjIyM\nEB8fjyFDhsDGxgYjRox46Ozgzs7O+Oabb1TPT506hTFjxqBnz57o1q1bvS5EIyMjHDlyBFevXsXY\nsWMhl8vRqVMndOrUCZs3b27WMT34/vn5+TAyMsLVq1fVrl9eXo4333wTTzzxBKytrTF27Nh6xycW\nizF//nxMmTIFvXr1QmxsbLPqIiLDwJYdgNWr9yAvb2m9ZXl5S/H3vy+Cp2fTW2WnTu1BSUnD/cTF\nLWpW666xlpeVlRW2b9+OgIAADB06FLW1tfjoo4+QlZUFc3Nz1XofffQRdu3aBXt7e8ybNw+hoaE4\nc+ZMve7LOvd3a167dg1BQUF444038MMPP8DExASHDh1qsI29vT3S0tIQHBxcr2tT0+N58P2bIioq\nCteuXcP27dvx+OOP4/3338czzzyDM2fOwNjYGACwatUqrF27FvHx8aitrW3yvonI8LBlB6CqSn3m\n19Yaa7Sf2lr1+6ms1Gw/gDIYli5dCmtra9XDxsZG9bqHhwdWr16N8PBwTJkyBWvWrGlwJ4iXXnoJ\nTzzxBMzMzPDvf/8b58+fVw0seVjwbN68Gaamppg/fz7Mzc1hamqKESNGNFpnU49n1qxZ9Y7Hy8ur\nSds+6ObNm0hMTMSKFSvg6+uL7t2744MPPsClS5fqDZwZPHgwIiIi0LFjx3r/CSAi/Xf+/Hk8++yz\nTV6fLTsAHTvK1C7395cjLa3p+wkJkWHPnobLzczkGtckEomwcOHCh56ze+mll/D222/D0tISU6dO\nbfC6n5+f6ue67r6CgoJHvnd+fj4CAgJgZKS9/wuJRCJ8+eWXWjlnd/nyZQBASEhIveWmpqaqm/yK\nRCIMGzasxe9FRLpp9erVSNPgFzRbdgDmzBkNF5cF9Za5uMzH7NmaDXXX1n6aavbs2XB3d4eVlRVi\nYmIavJ6dna36OS8vDyUlJao7vD+sy7B37944cuQI5PJHh7S2ArFTp071zs8dO3as0XUff/xxAMC+\nfftQUlKiepSWluJvf/ubaj0OSCEyTHK5HPv371f7n/zGMOwAjB8fiFWrQhASsghBQTEICVmEVavG\naHyeTVv7qfOwLsJNmzYhJSUFSUlJSEpKwqpVq7Bv375663z77be4dOkSKisr8fHHH8Pd3R2DBg16\n5L6nTp2K6upqvPnmm8jLy0N1dTXS09PVrmtnZwe5XI7jx4+36Hh8fX3xv//9D3/88Qd+/fVXrFmz\nptF1e/bsicmTJ+Odd97Bzz//DLlcjlu3buGHH35AWVmZ6r047yWRYTI2NsaZM2ewcuXKpm+kMECN\nHZY+Ha5YLFZ07NhRYWVlVe+RkpKiOHv2rKJz586K/fv3q9bfsmWLwtbWVlFYWKhQKBQKkUikiI+P\nVwwZMkRhbW2tEIvFipycHNX606dPV0RFRameOzs7K7755hvV8xMnTiiCg4MV1tbWChsbG8WYMWNU\nr4lEIsXhw4dVz1977TWFjY2NomvXrootW7aoPR5nZ2eFmZlZg+M5c+aMQqFQKAoKChQjR45UdO7c\nWTFo0CDF999/rzAyMlJIpVKFQqFQxMTEKEaNGqXaX3l5uWLhwoWKPn36KKysrBROTk6KKVOmKMrL\ny1Wf39KlSx/6GevT94GI1Gvqv2Pe9cBAGRkZ4aeffkJAQIDQpeis9vR9IDJUvOsBERHRXxh2RERk\n8HjpgYHiRdREZGjeeustuLu74+WXX9Z4W8FadpGRkbC1tcVTTz3V4DW5XA5vb29MmDBBtSwzMxM+\nPj7w9PREXFxcW5ZKREQCu3LlCj755BOcO3euWdsLFnYvv/xyoxcErlq1Cv369VNdCyaXyxEZGYnv\nv/8ex48fx4YNG3D+/Pm2LJeIiAT02WefAVBeX9wcgoXdsGHDYG1t3WB5QUEBUlNTMXPmTNUIm+zs\nbLi6usLZ2RmmpqYIDw9HcnJyW5dMREQCKC0txZdffonnnnsOjz32WLP2oXMDVObNm4eVK1fWm5lD\nKpXCyclJ9dzR0RFSqVSI8oiIqI1t2rQJt27dwty5c5u9D50aoLJr1y707NkT3t7eyMjIaNG+7p8+\nSywWQywWt2h/REQkjP3798PPzw+DBw9GRkZGs/JBp8LuyJEj2LFjB1JTU1FZWYk7d+4gIiICr732\nmmqCXwCQSCSqOR4bo26uSCIi0j/bt29HUVERRCJRg8bL4sWLm7QPnerGXLZsGSQSCS5fvoxt27Zh\nxIgR2LRpEwYMGICcnBzk5+ejuroaSUlJCA0NFbrcViMWi+vdmLWkpASBgYEICgrC7du3BayMiKjt\niUQidO/evUX7ECzsJk2ahICAAFy8eBFOTk5ISEhosE7daEwTExPEx8cjLCwMvr6+iIyMbHDvNkNy\n/41MJRIJhg0bBltbW+zduxddunQRuDoiIv0jWNglJibi6tWrqKqqgkQiaXCRYFBQEHbs2FHv+YkT\nJ3D69GnMmTOnrcsVxJkzZzB48GCMGDEC27dvR4cOHbBx40b06dMHGzduhLu7O2xsbPCPf/yj3kXk\np06dwogRI2BjYwMXFxcsXbpU9fqcOXMwa9Ys1bqBgYFwdnZWPf/vf/+L8ePHA1B2BY8aNQorVqxA\n7969YWtry+5hItJLOtWNSfccOXIEgYGBeP3117F69ep6r/3xxx84evQoDhw4gF27dmHTpk3Ytm0b\nAOD27dsYNWoUgoKC8OeffyIlJQVffvklPvnkEwDAqFGjVLfrKS0txYkTJwAAOTk5AIC9e/di1Kh7\n99/LzMxERUUFjh8/jnXr1mHJkiU4cuRIqx8/EZE26dQAFV3Q2KjNxkb/aLp+UygUCvz0008wNzfH\npEmT1L6+dOlS2NjYwM7ODkOHDsUvv/yCyZMnIyUlBbdv38Zbb70FU1NTPPnkk5gyZQq++uorvPnm\nmwgKClKdFz137hwGDhyIPn36YO/evXjsscdw5MgRxMbGqt6rS5cueP/992FkZIRnn30Wrq6uOH78\nOO+mQEStas2aNbh+/TpiYmK0cpNotux0kEgkwhtvvIGQkBAMGzYMFy9erPe6s7MzbGxsVM8dHBxQ\nWloKQHmOz8PDA+bm5qrXBw0apBrN2rlzZ/j5+SE9PR379u3D6NGjERwcjL179+Knn35Cp06d0L9/\nf9W2Tz/9dL0vmoODA+7evdsqx01EBADV1dVYtmwZsrOztRJ0AFt2DWjaImvp9YCNMTY2xsaNG/HP\nf/4TgYGB2Lt3r9p5RB/k5OSEc+fOoby8HBYWFgCArKyserMO1IXbhQsXsHHjRjg7O+OVV15B3759\nERwc3CrHQ0TUVP/73/9w7do1xMfHa22fbNnpqLqp0j777DNERERALBbj2LFjja5bt/748ePRuXNn\nfPzxx6ipqcHvv/+Obdu2YcaMGar1g4ODkZaWhsLCQvj4+MDGxga9e/fGunXrmhR2vOEpEbUWhUKB\n2NhYPPnkkxg9erTW9suw01F1lx4AyhGS0dHRCA4OhkQiqfda3bp1y7p06YI9e/bgwIEDsLW1xZgx\nYxAZGYl58+ap1vf394dCocCIESNUy4KDg3H37t16YXf/fhurjYhImw4fPoxff/0Vc+fO1VoXJgCI\nFAb43/TGbtPe1Nu3U/vA7wOR7vnwww8RGxsLiUSiOhXzME39d8ywo3aL3wci3VRUVIRu3bo1aV2G\nHcOOHoHfByL919R/xzxnR0REBo9hR0REBo9hR0REBo9hR0REgkpNTcUrr7yCW7dutdp7cIAKtVv8\nPhDphhEjRiA3NxeXLl2CiYlmE3txgAoREem8kydP4sCBA5g9e7bGQaeJdjU3prW1NWf/IBVra2uh\nSyBq91atWgULCwvMnDmzVd+nXYVdcXGx0CUQEdFf/vzzT2zduhUzZ85s9f98shuTiIgEcfjwYQBA\ndHR0q79XuxqgQkREuqWkpKRFrTpOF2Z4h0VERA/gaEwiIqK/CBZ2kZGRsLW1rXf3bYlEguHDh8PD\nwwNisRgbN25UvZaZmQkfHx94enoiLi5OgIqJiEhfCdaNeejQIVhZWSEiIgKnT58GABQWFqKwsBBe\nXl64efMm+vfvj4yMDPTp0wdubm5IT0+Hg4MD/Pz8kJiYCHd3d7X7ZjcmEVH7oPPdmMOGDWtwUtLO\nzg5eXl4AgO7du8PPzw9SqRTZ2dlwdXWFs7MzTE1NER4ejuTkZCHKJiKiFsjLy8Ozzz6L3NzcNn1f\nnT1nl5ubi7Nnz8Lf3x9SqRROTk6q1xwdHSGVSgWsjoiImiMuLg6pqamwtLRs0/fVyYvKS0tLER4e\njtjYWFhaWjZr1pOYmBjVz2KxGGKxWHsFEhGRxm7fvo0NGzYgPDwcvXr1atY+MjIykJGRofF2Ohd2\nNTU1eP755zF16lRMnDgRAODg4ACJRKJaRyKRwNHR8aH7uT/siIhIePHx8SgtLW3RReQPNl4WL17c\npO10qhtToVBgxowZ8PDwwNy5c1XLBwwYgJycHOTn56O6uhpJSUkIDQ0VsFIiItKEXC7H6tWrMWzY\nMPj6+rb5+wsWdpMmTUJAQAB+//13ODk5ISEhAUeOHMGWLVuwf/9+eHt7w9vbG2lpaTAxMUF8fDzC\nwsLg6+uLyMjIRkdiEhGR7jl//jyKi4vrNWTaEmdQISKiNnH37l2Ym5tr9VY+nC7M8A6LiIgeoPPX\n2REREbUVhh0RERk8hh0RERk8hh0REbWKiooKjBs3DpmZmUKXwrAjIqLWsXXrVuzevRu1tbVCl8LR\nmEREpH0KhQKenp4wNjbGiRMnmjXtY1M09fe9zk0XRkRE+m///v04c+YMEhISWi3oNMGWHRERad0z\nzzyDY8eO4cqVK+jYsWOrvQ+vsyMiIkGUlpbizJkzeO2111o16DTBlh0REWmdTCZDdXU1LCwsWvV9\nOF2Y4R0WERE9gN2YREREf2HYERGRwWPYERGRwWPYERGRVkydOhXffPON0GWoxbAjIqIWy87Oxjff\nfIOioiKhS1GLozGJiKjFJk+ejJSUFBQUFKBTp05t9r4cjUlERG2ioKAA27dvx4wZM9o06DTBsCMi\nohb5/PPPUVtbi9mzZwtdSqMEC7vIyEjY2triqaeeqrc8MzMTPj4+8PT0RFxc3COXExGRsLKzs/Hs\ns8+id+/eQpfSKMHO2R06dAhWVlaIiIjA6dOnAQByuRxubm5IT0+Hg4MD/Pz8kJiYiL59+6pd7u7u\nrnbfPGdHRNR2FAoFSktLBenC1PlzdsOGDYO1tXW9ZdnZ2XB1dYWzszNMTU0RHh6O5ORkHDt2TO1y\nIiISnkgk0tlzdXV06pydVCqFk5OT6rmjoyOkUmmjy4mIiJpCp27eqs0b/MXExKh+FovFEIvFWts3\nEREJIyMjAxkZGRpvp1Nh5+DgAIlEonoukUjg6OjY6PKHuT/siIjIMDzYeFm8eHGTttOpbswBAwYg\nJycH+fn5qK6uRlJSEkJDQxtdTkREwnj33XexZMkSoctoMsHCbtKkSQgICMDFixfh5OSEhIQEmJiY\nID4+HmFhYfD19UVkZCTc3d0bXU5ERG2vqKgIq1atqtfjpus4XRgREWlkxYoVePfdd3H69Gn0799f\n0Fp4p3LDOywiIsHV1NSgd+/ecHd3x969e4Uup8m/73VqgAoREem2b7/9FlKpFOvWrRO6FI3o1AAV\nIiLSbcePH0ffvn0xduxYoUvRCLsxiYhII6WlpbCyshK6DAA8Z8ewIyJqB3R+bkwiIqK2wrAjIiKD\nx7AjIiKDx7AjIqKH+uqrrzBz5kxUVlYKXUqzcYAKERE1SiaTwcPDA127dkVWVpZW706jDbyonIiI\nWmzJkiW4ePEifvjhB50LOk2wZUdERGplZmZi+PDh+Pvf/46NGzcKXY5avM7O8A6LiKjNlJSUwNPT\nE2ZmZvj111/RqVMnoUtSi9fZERFRs5mZmWHixIlITEzU2aDTBFt2RESkt9iyIyIi+gvDjoiIDB7D\njoiIAADV1dVCl9BqGHZERISUlBT0798fly5dErqUVtHkASrV1dXYuXMndu7cCZlMBmNjY9y9exc2\nNjYYPXo0XnjhBRgZ6UZ2coAKEVHTXbt2DZ6ennBwcEBWVhbMzMyELqnJtHqd3eHDh5GSkoKpU6ei\nT58+MDU1Vb1WUVGBM2fOYMuWLZg2bRp8fHxaVrkWMOyIiJqmtrYWISEhOHz4MI4fPw53d3ehS9KI\n1kZjVlVVoUuXLli2bBn69etXL+gAwNzcHH5+fli1ahXkcnnzK77P+vXrERAQAF9fX8ydO1e1PDMz\nEz4+PvD09ERcXJxW3ouIqD37+OOPkZ6ejk8//VTvgk4TGl9nV1hYCDs7OwBAeXk5LCwstFpQcXEx\nfH19cebMGZibm+OZZ57B3LlzMXLkSLi5uSE9PR0ODg7w8/NDYmKi2r8ctuyIiB4tPz8fffv2xYQJ\nE/Dtt9/q5dyXWr/ObtmyZdi9ezd27typWnb27FkcOnSoeRU2wtzcHAqFArdv30ZFRQXKy8vRtWtX\nZGdnw9XVFc7OzjA1NUV4eDiSk5O1+t5ERO2Js7Mztm3bhvXr1+tl0GmiyWEXFhaGy5cv44svvsCE\nCRMQFRWFkydPYu/evVotyNzcHGvXroWzszPs7OwwZMgQDBw4EFKpFE5OTqr1HB0dIZVKtfreRETt\nzXPPPQcbGxuhy2h1Tb7Fj7u7O9zd3eHi4oKQkBAUFhbi2LFjGDJkiFYLunHjBl599VWcO3cO1tbW\nePHFF5GSkqLx/zpiYmJUP4vFYojFYq3WSUREbS8jIwMZGRkab/fIc3ZVVVW4e/cuunfv/sid5eXl\nwcXFReMi7peSkoLNmzdj27ZtAIC1a9ciPz8fYWFhiImJQVpaGgBg+fLlMDIywttvv91gHzxnR0TU\nPmjtnF3Hjh2RlZWFrVu3oqKiQu06N27cwIIFC7RyMeKwYcPwyy+/oLi4GFVVVdi9ezdGjx4NPz8/\n5OTkID8/H9XV1UhKSkJoaGiL34+IqL2oqanBuXPnhC5DEE0ejXnt2jUkJCTg+vXrqKysRGVlJe7c\nuQMzMzN4eXlh1qxZ6NKli1aK2rhxIxISElBeXo4xY8Zg8eLFMDIywsGDBzF37lzIZDJERUVhzpw5\n6g+KLTsiogYWLFiAjz76CGfPnoWrq6vQ5WhFq928ddq0aejZsycCAgIwePBg1WUIuoRhR0RU34ED\nBzBy5Ei8/PLL2LBhg9DlaE2r3qn8woULyMrKws8//4wTJ07ghRdewJtvvsnpwoiIdFBRURGefvpp\nWFpa4tdff4WlpaXQJWlNq4XdiRMnUFlZicGDBwMAdu/eDRcXFxw6dAgzZsxoXrVaxrAjIlJSKBR4\n7rnnkJKSgqysLJ2Y0lGbmvr7vsmXHtRJS0uDkZERPv30U1hYWMDBwQFmZmbo0KFDswolIqLWk5ub\ni3379mHFihUGF3Sa0Lhld/78edy9excDBw5ULdu4cSN69OiB8ePHa73A5mDLjojonoKCAtjb2+vM\nqSZtatVzdrqOYUdE1D5ofW5MIiIifcWwIyIig8ewIyIyIFevXsWXX37JUzkP4Dk7IiIDIZfLMXr0\naGRlZeH8+fN47LHHhC6p1bXapQdERKSbVq5cif379+Orr75qF0GnCbbsiIgMQHZ2NoYMGYKwsDAk\nJSUZ/M1Y6/DSA8M7LCIite7cuQNvb2/IZDL89ttv6Nq1q9AltRleekBE1E7cunUL3bp1w9atW9tV\n0GmCLTsiIgOgUCjaTdfl/diyIyJqR9pj0GmCYUdERAaPYUdERAaPYUdEpGf279+P8PBw3L59W+hS\n9AYvKici0iM3b97E1KlT0aVLF5iY8Fd4U/GTIiLSEwqFApGRkSgqKkJqaiosLS2FLklvMOyIiPTE\n559/jp07dyI2NhZeXl5Cl6NXdPKcXVlZGaZNmwZvb2/069cPR48eBQBkZmbCx8cHnp6eiIuLE7hK\nIqK2c+7cObzxxhsYO3YsoqOjhS5H7+jkReXTpk1DUFAQIiMjIZPJUFZWBisrK7i5uSE9PR0ODg7w\n8/NDYmIi3N3dG2zPi8qJyNBUVlZi8eLFmDdvHnr27Cl0OTpDb+fGvH37Nry9vXHp0qV6y3/++Wcs\nXrwYaWlpAIAVK1YAAN55550G+2DYERG1D3o7g8rly5fRo0cPTJ8+Hf3790dUVBQqKioglUrh5OSk\nWs/R0RFSqVTASomISF/o3AAVmUyGY8eOYeHChVi7di1mzZqF7du3azzqKCYmRvWzWCyGWCzWbqFE\nRNTmMjIykJGRgYsX/8DRo3lN3k7nws7R0RHdunXDhAkTAACTJk3Cpk2bEB0dDYlEolpPIpHA0dGx\n0f3cH3ZERPpGLpejvLwcnTp1EroUnSIWi1FWZoQtW37EpUsJAJo2J6jOdWPa2dnB1dUVR48eRW1t\nLVJSUhAcHAw/Pz/k5OQgPz8f1dXVSEpKQmhoqNDlEhG1iv/85z/w8vLCjRs3hC5F56xevQd5eUs1\n2kbnwg4Avv76a0RHR6Nv376QSqUIDw+HsbEx4uPjERYWBl9fX0RGRqodiUlEpO+ysrLw3nvvwc/P\nD927dxe6HJ1TXKx5p6TOjcbUBo7GJCJ9VTcivba2FidPnuTNWB+weTMwffpC1NZ++NcSPR2NSUTU\nXikUCrz22mu4cuUK7zr+gPJyYMYMICICcHcfjccfX6DR9gw7IiIdkZWVha1btyImJgYBAQFCl6Mz\nLlwABg0CEhKAhQuBkycD8dlnIQgJWdTkfbAbk4hIh+zbtw9isRjGxsZCl6ITtmwB/vEPwMJC+fPo\n0fVf19sZVLSBYUdEpN/Ky4E5c4ANG4DAQCAxEbC3b7ie3s6gQkRE7Vtdt2V8PLBgAbBvn/qg0wTD\njohIAMXFxTh27JjQZeicb74BBgwACguB3buBDz8EtHGPWoYdEVEbu3jxIvz9/REaGoqKigqhy9EJ\nFRVAVBQwdSrg6wucPAmEhGhv/ww7IqI2tH//fvj7+6OkpATbt2+Hubm50CUJrq7b8quvgPnzld2W\nDg7afQ+GHRFRG1m3bh1CQkLQq1cvHD16FEOHDhW6JMHVdVteu6bstly6VDvdlg/iaEwiojZw8eJF\n9OvXD6NGjcK2bdvQpUsXoUsSVEUFEB0NrF8PDB2qHG35kLn9G8VLDwzvsIhIz/3000/w9/eHSWs0\nXfTI778DL74InD4NvPsusGRJ81tzDDvDOywiIr23dSvwyiuAmZnyIvExY1q2P15nR0REOqOiQhly\nU6YA3t7K0ZYtDTpNMOyIiLRs8+bNSElJEboMnfH774C/v/L83DvvAAcONO/8XEsw7IiItKS2thbz\n589HREQE1q9fL3Q5OiExUTnaUioFUlOB5ctbZ7TlozDsiIi0oKysDC+88AKWL1+OV155Bdu3bxe6\nJEFVVACzZgGTJwNPP63sthw7Vrh62veQICIiLSgoKMCECRNw6tQpxMbGIjo6GiKRSOiyBHPxonK0\n5alTym7LJUsAU1Nha2LYERG10JUrV1BQUICdO3di3LhxQpcjqMRE5UCUjh2BlBRAVz4OXnpARKQF\nZWVlsLS0FLoMwVRUAPPmAevWAQEBwLZtgJNT678vLz0gImpD7TnoLl4EBg9WBt3bbwMZGW0TdJrQ\n2bCTy+V/+fGXAAAXr0lEQVTw9vbGhAkTVMsyMzPh4+MDT09PxMXFCVgdEbVXtbW1QpegU7ZtU96l\nQCIBdu0CVqwQ/vycOjobdqtWrUK/fv1UJ3nlcjkiIyPx/fff4/jx49iwYQPOnz8vcJVE1J4UFhZi\n6NCh7X6kJQBUVgKvvgpMmgQ89ZRytOX48UJX1TidDLuCggKkpqZi5syZqr7Y7OxsuLq6wtnZGaam\npggPD0dycrLAlRJRe3Hq1CkMHDgQv/32W7uf2zInR3mR+BdfAP/+N3DwoO51Wz5IJ8Nu3rx5WLly\nJYyM7pUnlUrhdN+n6ejoCKlUKkR5RNTO7NixAwEBAaitrcWhQ4cQFhYmdEmCSUoCfHyU3ZY7dwL/\n/a9udls+SOf+e7Jr1y707NkT3t7eyMjIaPZ+YmJiVD+LxWKIxeIW10ZE7c+6devw6quvwtfXF8nJ\nybC3txe6pDaVkpKJ1av3oKLCBJcvy1BQMBqDBwdi2zbgscfavp6MjIxmZYPOhd2RI0ewY8cOpKam\norKyEnfu3EFERARee+01SCQS1XoSiQSOD5lc7f6wIyJqLj8/P/z973/H2rVrYWFhIXQ5bSolJRPR\n0T8iL2+palnXrgvw9tvAY48FClLTg42XxYsXN2k7nb7O7uDBg/joo4+wc+dOyGQyuLm5Yd++fbC3\nt8fAgQORmJgId3f3BtvxOjsiopYbPnwhMjI+bLA8JGQR0tI+EKCihpr6+17nWnYPqhuNaWJigvj4\neISFhUEmkyEqKkpt0BERUcvcuAF8/DGQmak+Iiorjdu4opbT6bALCgpCUFBQvecnTpwQsCIiMmQ5\nOTno06eP0GUIpi7k1qwByssBW1sZCgsbrmdmJm/74lpIJ0djEhG1tTVr1sDd3R3fffed0KW0uZs3\nlRM29+6tHF0ZGgqcPQt89dVouLgsqLeui8t8zJ49SqBKm0+nW3ZERK1NJpMhOjoan3/+OUJDQxES\nEiJ0SW3m5k1lSy4uTtmSCw8HFi0C6s4QubsrB6HExS1CZaUxzMzkmD17DMaPF2ZwSkvo9ACV5uIA\nFSJqilu3buGll17C3r178e9//xvLly+HsbH+nY/S1KNCTp809fc9w46I2q3x48dj7969+OKLLxAZ\nGSl0Oa3uwZD729+UIdevn9CVNR/DzvAOi4i07Ny5c7hx40a9gXCG6OZN4JNPlCFXVmYYIVeHYWd4\nh0VEpBFDDrk6BnOdHRFRS12+fBkdOnSAg4OD0KW0iaKie92VhhpymuKlB0RksCQSCWbNmoW+ffu2\niykEi4qABQsAZ2flfeXGjwdOnwYSE9t30AFs2RGRAbp69SqWLVuG9evXQ6FQ4JVXXsH8+fOFLqvV\nFBUpuytXr1a25F56SdmS8/AQujLdwXN2RGRQ7ty5AycnJ5SXlyMyMhILFizAY0JMz98GGHIcoMKw\nI2rHNm/ejCFDhuCJJ54QupRWUVQExMYqQ660tH2GXB2GneEdFhG1c8XF91py7T3k6nA0JhEZrFu3\nbiE2NhZ5eXnYsmWL0OW0ugdD7sUXlSHXv7/QlekPjsYkIr1x584dfPjhh+jduzeWLFmCyspKVFdX\nC11WqykuVoaaszOwbBkwdixw6hSQlMSg0xRbdkSkF9asWYP3338fxcXFCA0NxeLFi+Hl5SV0WVqT\nkpKJ1av3oKrKBEZGMvTsORqpqYFsyWkJw46I9ML169fh7++PxYsXY8CAAUKXo1UpKZmIjv4ReXlL\n71u6AEOHAmvXBjLktIADVIhIL9TW1sLIyPDOvJSXAwEBC/Hbbx82eC0kZBHS0j4QoCr90dTf94b3\nzSEivVVVVYWkpCS1v7wMKejKy4HvvlNO49WjB/Dbb+o72SorDf92Q23FcL49RKS3ampq8OWXX6JP\nnz4IDw/HTz/9JHRJWlcXcOHhQM+ewAsvAAcOABERwIABMrXbmJnJ27hKw8WwIyLByGQyJCQkwM3N\nDbNmzYK9vT327NmDoUOHCl2aVqgLuP37galTgX37gKtXgbVrgZiY0XBxWVBvWxeX+Zg9e5RAlRse\nDlAhIsFs3LgRUVFR8PX1xZo1azB27FiIRCKhy2qRigpg927gf/8Ddu1STuPVo4cy4F58EQgKAkwe\n+M07fnwgACAubhEqK41hZibH7NljVMup5XRugIpEIkFERASuX7+OHj16YPr06Zg+fToAIDMzE3Pn\nzoVMJkNUVBRmz56tdh8coEKkHyorK5Geno7x48frdcipC7ju3YHnn2884Eg79Ha6sMLCQhQWFsLL\nyws3b95E//79kZGRgT59+sDNzQ3p6elwcHCAn58fEhMT4e7u3mAfDDsi3VJbWwvAsAaZ1AXc9u3A\nzp33Au6555TTeDHg2obejsa0s7NTXSjavXt3+Pn5QSqVIjs7G66urnB2doapqSnCw8ORnJwscLVE\n9DAKhQLJycnw8fHBN998I3Q5LVZRAXz/PTBpkrJr8vnngfR0YMoUYO9e4No1YN06YORIBp2u0em/\njtzcXJw9exb+/v7YvXs3nJycVK85Ojri6NGjAlZHRI25cOECkpOTsW3bNpw8eRKurq7o0qWL0GU1\nS0UFkJZ2r4uytBTo1k0ZcC++CIjFDDZ9oLN/RaWlpQgPD0dsbCwsLS017s+//67EYrEYYrFYuwUS\nkVpHjx6Fv78/AMDb2xsJCQmYOnUqTPQoEeoCrq6Lsi7gJk1SdlEy4ISTkZGBjIwMjbfTuXN2gPKa\nm2eeeQZjx47F3LlzAQBZWVmIiYlBWloaAGD58uUwMjLC22+/3WB7nrMjan1VVVXo2LFjg+VyuRzr\n16/HuHHjdPKmqffPQdmxowxz5ozG+PGBjQbcc88pW3DDhzPgdJHeDlBRKBSYNm0aunfvjk8++US1\nXCaTwc3NDfv27YO9vT0GDhzIASpEbezmzZvYtWsXkpOTsWfPHvz+++9wdHQUuqwmUzcHpa3tAri5\nheDXXwMbBJxYDJiaClcvPZre3s/u8OHD2LJlCzw9PeHt7Q1A2YobM2YM4uPjERYWprr0QF3QEZH2\nff3119iwYQMOHz6M2tpaODo6Yvr06apRlvpi9eo9D0y2DPz551IUFS3Cyy8HMuAMmM6F3dChQxv9\nBxQUFIQTJ060cUVEdPLkSdy+fRsLFizAxIkT4ePjoxfXxdXUKO//9vPPQFYWcPCg+l95gwcb48sv\n27g4alM6F3ZE1PYqKyuxf/9+GBkZYcyYMQ1eX7lypV4MMLl+XRlsdY9fflFO2QUAvXoBXbrIcP16\nw+0sLDgHpaHT/W8vEbWK4uJipKSkIDk5GWlpaSgrK4NYLFYbdroYdDLZvVZb3ePSJeVrJiaAtzcw\ncyYweLDy8dhjQGrqaERHL6jXlamcg7LhMZNh0bkBKtrAASpED5eTkwN3d3fI5XL06tULoaGhmDhx\nIkaMGKF2hKUuuHGjfrAdO3av1WZndy/UBg8GfH0Bc3P1+0lJyURc3N775qAcxTko9ZjejsbUBoYd\nkZJCoVB7bk2hUGD58uUIDg7GgAEDdG4aL5kMOH26frjl5Slfq2u11QWbvz/w+OOAHpxCpFbAsDO8\nwyJqkurqamRkZCA5ORk7duzAvn370LdvX6HLeqgbN5QDSO5vtZWVKV/TpNVG7Y/eXnpARM3z448/\nIiEhAbt378adO3dgYWGBkJAQVFdXC1JPYxdvy2TAmTP1W225ucptTEwALy8gMvJeuLHVRtrAsCPS\nIwqFAjU1NejQoUOD1zIzM3HgwAG8+OKLePbZZzFy5EiYC9QEUnfxdnb2Ajg5AZcuBTZotUVF3Wu1\nWVgIUjIZOHZjEukouVyOixcv4uTJk/Ue4eHhWLVqVYP1y8rKYGZmBmNj4zavtbYWkEqBnBzg4kVg\n+fKFuHLlwwbrde68CNOmfcBWG2kNuzGJ9FxSUhKmTJkCAOjQoQM8PDwwfvx4BAUFqV3f0tKy1Wsq\nKlKG2YOPnBzl5Ml1RCL1v1q8vY2xenWrl0nUAMOOqA0pFApcu3atXkvNysoK8fHxDdYVi8X4+uuv\n4eXlhSeffFJt12VrKCtTnkNTF2rFxffWMzEBnngC6NsXCA5W/ln3iIyUYc+ehvs2M+PF2yQMhh1R\nG7l06RL8/f1x48YN1bLevXtj+PDhate3t7dHREREq9RSUwPk56sPtIKC+us6OioD7KWX6geas3Pj\nc0jOmTMaeXm8eJt0B8/ZEWlBaWkpTp06hZMnTyI3N7feHTvqVFdX49VXX8XTTz8NLy8veHp6omvX\nri1+78ZGPSoUyvNo6gLt8mXltWx1rK0BN7f6Yda3L+DqCjS3d5QXb1Nb4HV2hndYpGMUCgUmT56M\n48ePIzc3V/Wds7GxwaVLl9rkztw//JCJefN+xB9/3GtBWVouQI8eIbh+PVA1wwigvDatT5+Ggda3\nr/K2NkT6iGFneIdFbaCmpgZXr16FRCLBlStXIJFIIJFIsGjRItja2jZYf9y4cbCwsICXl5fq4eDg\n0OI7AlRUANeuAVevKv9s7Ofi4oUAGo567N59EaZO/aBeoDk4ADo2UQpRi3E0JtEDamtrUVhYCIlE\nAnd3d3Tu3LnBOoMGDWpwG6muXbtixowZ+OWX3xt0F6ampmpUQ2npw8Or7ufbtxtua2KinLm/Vy9l\n92JgIJCaaoI//mi4roeHMWJjNSqNyKAx7MigLV68GPv27cOVK1cglUoh++tE1d69exEcHNxg/Tff\nfBMVFRVwcnJSPaysrNReJJ2XtwAAMG5cIO7ceXSIXbsG3L3bsMYOHQB7e2WIubsDI0feC7Veve69\n1q1bw5ZZSIhMbdhx1CNRfQw70ivZ2dk4efKkqnux7rF69Wq1t6YpKSmBSCTC0KFD6wWYl5eX2v1P\nnjwZcrmyZVVSAly4ANy6Bcyf3/AO13l5S/Hii4sgEtU/N1bH3PxeWD39NDBmzL3guj/ErK2bf2E1\nRz0SNQ3DjgSVn5+PS5cuobi4GCUlJSguLsaxY79BIukAc3PneqMLAeDzzz/H119/DSMjI9jb28PJ\nyQne3t4NuiRrapQh9dprn6KkRPlzSYnycfo0kJl57/mDr6vrQmzsn4qVlTGmTlUfYp07t/7sIHWf\nS1zcovtGPY7hqEeiB3CASjvW2JD1lvj555+RnZ2tCq66EJs5cybCwsIarP/666/js88+e2CpMYDP\nAMwCADg5LcDs2SHw8AhEbq4Et24pUFtrjzt3TNSGVUmJ8tzYw5iZKVtU9z+6dm24rO7xzjsLceRI\nw4EgISGLkJb2QfM+LCJqMY7GNLzD0qp756A+BHAXgAguLiuwalVIvcDbsWMHdu7cWS+4iouL8fbb\nb+Of//wnFArlyMGyMuVjyZK3kJCwEgBgadkVFhbWsLCwQWDgv/Dkk5NRWnpv3bIy4Nq187h9+zpk\nMhtUV1vjjz/iUFW1AsCDTaJFAOqHiqVl4+H0sODq2lUZds37vOp3F65axVYUkZDafdiNHr1AKy0V\nbWhJC6pulnuFQqH2DtJnz57Fr7/+ivLycpSVlaGsrAzl5eUIChJj+PAQVFUBlZVQ/ZmQsAabNsXi\n+vU/UVNjBKAMQC2AxQDeQ9++i/DSSx+owujo0WW4cCEOJiY2MDKygUhkg9paa5iYvASZbBzKyoD6\n36Bbf+2vC5QttPpMTJQhZWWl/LPuUff80KEYXL8e02C7p56Kwfr1MfUCq7HZO1oLL5Im0j0GGXaZ\nmZmYO3cuZDIZoqKiMHv2bLXrKa9xUsDFZUGDlkpLVVVVoaysDFVVVaiurkZVVRWqqqrRqZMNevSw\nR00N6j0SEjYiLu47FBc/A6AcQBksLXdj4sSJeO65txqEUWbmFzh48CPU1JSjuroMNTVlUCjk8PBY\nADe3DxusL5EsxZ9/LnygSlMACwC8r+YIvgPwPYALAAIBWAIoAjAdwCAAMRCJYh4aSC15/qjpHUNC\nFmLPHt3uLszIyIBYLBa6DL3Bz0sz/Lw0Y3DX2cnlckRGRiI9PR0ODg7w8/NDcHAw3N3dG9niA+Tl\nGWHKlHfg6joQ9vZiPP74sw3C6NKlROTmfga5vApyeTXk8irU1laje/co9OjxboP1b936FKWl76h5\nv7cA/EfN8s0A9gPYpVpSVmaKrVtF2Lr1LTXr9wIwCEZGljA2toS5uQVMTCxx69ZQXLwIdOyo7ILr\n2FE5AKJnz1dgZPQSLCwsYWlpCUtLC1hamtZbr/6fz8PM7HnExCzE8eN1oRIDZdABwcFy7Nkj3G1X\n9GF0IX8ZaYafl2b4ebUOvQm77OxsuLq6wtnZGQAQHh6O5OTkh4TdewBEuH3bGMePn8WpU51gZfUs\nTE1R71FeboTq6o4wMuoMU9MOMDfvCGPjDujevTecnNBg/Vu3RuLGjU/RoUNHdOjQQfVnr17ucHRs\nuP7y5QOQm/sVgI4ALKBsSZnC2zsGGzeqC6SJ6NBhogYzXfT466EZkWg0oqMbhsrcuWMEvb8YRxcS\nUWvQm7CTSqVwcnJSPXd0dMTRo0cfskUNAGOEhLz3iO6vv/31aKoBfz2aJinJFLm5vRss79lTDk9P\nDd5Wy+4PlQsXDuHJJxfpTKiMHx+oE3UQkeHQm3N23333HdLS0rB+/XoAwJYtW3D06FHExcU1WFck\ncgWQ18YVEhFRW3NxcUFubu4j19Oblp2DgwMkEonquUQigaOjo9p1FYpHHzgREbUfejMH+oABA5CT\nk4P8/HxUV1cjKSkJoaGhQpdFRER6QG9adiYmJoiPj0dYWJjq0oPGB6cQERHdozfn7IiIiJpLb7ox\niYiImsugwi4zMxM+Pj7w9PRUO0qT6ouMjIStrS2eeuopoUvReRKJBMOHD4eHhwfEYjE2btwodEk6\nrbKyEoMGDYKXlxf8/f0RyzvJNolcLoe3tzcmTJggdCk6z9nZGZ6envD29sbAgQMfub7BdGPK5XK4\nubnVm2ElMTGR5/Ue4tChQ7CyskJERAROnz4tdDk6rbCwEIWFhfDy8sLNmzfRv39/HDhwgN+vhygv\nL4eFhQWqqqrg6+uL//u//4Orq6vQZem0Tz75BMePH8fdu3exY8cOocvRab1798bx48dhY2PTpPUN\npmV3/wwrpqamqhlWqHHDhg2DtbW10GXoBTs7O9UNX7t37w4/Pz9cvXpV4Kp0m4WFBQCgtLQUMplM\n7UTmdE9BQQFSU1Mxc+ZM3rWliTT5nAwm7NTNsCKVSgWsiAxVbm4uzp49C39/f6FL0Wm1tbV4+umn\nYWtri9dff73ev09qaN68eVi5ciWMmj5XYLsmEokwYsQIeHt7qyYbeRiD+VRFQk7oSO1GaWkpwsPD\nERsbC0tLS6HL0WlGRkb47bffkJubi88//xwnTpwQuiSdtWvXLvTs2RPe3t5s1TXR4cOH8dtvv2Hr\n1q1YtmwZDh069ND1DSbsNJlhhag5ampq8Pzzz2Pq1KmYOHGi0OXoDWdnZ4wbNw4HDx4UuhSddeTI\nEezYsQO9e/fGpEmTsH//fkRERAhdlk7r1asXAMDd3R1hYWHIzs5+6PoGE3acYYVak0KhwIwZM+Dh\n4YG5c+cKXY7Ou3nzJm7dugUAKCoqwu7duznq9yGWLVsGiUSCy5cvY9u2bRgxYgQ2bdokdFk6q7y8\nHHfv3gUA3LhxA6mpqY/8funNDCqPwhlWNDdp0iQcPHgQRUVFcHJywpIlS/Dyyy8LXZZOOnz4MLZs\n2aIa6gwAy5cvx5gxunOfPV1y7do1TJs2DXK5HHZ2dvjXv/6FkSNHCl2W3uBpmYf7888/ERYWBgDo\n1q0b5s2bh9GjRz90G4O59ICIiKgxBtONSURE1BiGHRERGTyGHRERGTyGHRERGTyGHRERGTyGHRER\nGTyGHRERGTyGHRERGTyGHZEBOH/+PJYtWyZ0GUQ6i2FHZAAOHDigmsaMiBpi2BHpuR9//BEbNmxA\nQUEBCgsLhS6HSCdxbkwiAzBhwgTs3LlT6DKIdBZbdkR6rrCwEHZ2dkKXQaTTGHZEeu7YsWMYOHAg\njh07hvLycqHLIdJJDDsiPefi4gKJRILKykpYWFgIXQ6RTuI5OyIiMnhs2RERkcFj2BERkcFj2BER\nkcFj2BERkcFj2BERkcFj2BERkcFj2BERkcFj2BERkcH7fy2i34hOL4n5AAAAAElFTkSuQmCC\n", "text": [ "<matplotlib.figure.Figure at 0x10828e850>" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, so with $N = 10$ steps the numerical solution (blue dots) follows the known exponential behaviour, but is underestimating the known result (black dashed line). We know from the Taylor series truncation we made that the error in the approximation scales with the stepsize, so if we want a more accurate solution the first thing to do is to make the stepsize $h$ smaller. Let's try that." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Checking Accuracy\n", "\n", "To see that the accuracy improves with smaller $h$, we'll solve the system for different numbers of steps: $N = 5, 10, 20, 40$ (corresponding to $h = 1, \\frac{1}{2}, \\tfrac{1}{4}, \\tfrac{1}{8}$)." ] }, { "cell_type": "code", "collapsed": false, "input": [ "for N in [5, 10, 20, 40]:\n", "\n", " t = np.linspace(0., t_max, N+1) # Time steps\n", " y = ode_int_ee(exp, y_0, t, solve_args)\n", " plt.plot(t, y, '-o', label='%d steps'%N) \n", " \n", "plt.plot(t, np.exp(solve_args['a']*t), 'k--', label='Known')\n", "\n", "plt.xlabel(r'$t$'), plt.ylabel(r'$y$')\n", "plt.legend(loc=2)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 12, "text": [ "<matplotlib.legend.Legend at 0x108406210>" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAbYAAAEiCAYAAACV/vclAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVNX7wPEPq+CCgqgokCS4oqgoroHIF7fMtTRXQNTS\nCsX2XNJKM1MzpdIyFbfU/OUK7ilSWWLuC4YiKJKJCAmKLMPc3x/kJDIkM8Mqz/v14hXc5blnZnrN\n47n3nOcYKYqiIIQQQjwhjMu6AUIIIURxksQmhBDiiSKJTQghxBNFEpsQQogniiQ2IYQQTxRJbEII\nIZ4oktiEEEI8USSxCSGEeKKUSmILDAykXr16tGrVKt/2e/fu4e/vT9u2bWnRogVHjx4FIDIyEnd3\nd9zc3AgJCSmNJgohhHhCGJVG5ZGffvqJ6tWr4+fnx9mzZzXb/f396datG4GBgahUKu7du0f16tVp\n2rQpBw4cwN7eHg8PDzZs2EDz5s1LuplCCCGeAKXSY/P09MTa2jrftjt37vDTTz8RGBgIgKmpKTVr\n1iQqKgoXFxecnJwwMzNj2LBhbN++vTSaKYQQ4glQZs/Y4uLiqFOnDgEBAbRs2ZLx48dz//59EhMT\ncXR01Bzn4OBAYmJiWTVTCCFEBWNaVhdWqVQcO3aM6dOns3TpUl5++WU2b95MtWrVdIrj4uJCbGxs\nCbVSCCFEeeHs7Mzly5cfe1yZ9dgcHByoXbs2/fr1w9LSkuHDh7N7927s7e1JSEjQHJeQkICDg0Oh\ncWJjY1EURX6K+DNz5swyb0NF+pH3S94veb/Kz09ROzFlltjs7OxwcXHh6NGjqNVqwsPD8fX1xcPD\ng0uXLhEfH092djabNm2if//+ZdVMIYQQFUypJLbhw4fTpUsXYmJicHR0ZNWqVQCsXr2ayZMn06RJ\nExITExk2bBgmJiasXLmSQYMG0a5dOwIDA2VEpBBCiCIrleH+JcnIyIgK/hJKVUREBN7e3mXdjApD\n3i/dyPulG3m/dFPU73tJbEIIISqEon7fS0ktIYQQT5QyG+5f0mxsbEhNTS3rZohHWFtbk5KSUtbN\nEEI8wZ7YW5Fyi7J8ks9FCKEvuRUphBCiUpLEJoQQ4okiiU0IIcQTRRKbEEKIJ4oktnIgICAAc3Nz\natSooflZtmyZ3rHGjx9fzC0UQoiK44kd7v9fwsMjWbJkH1lZplSpomLSpJ707etVZnGMjIwICAjg\nm2++0flcIYQQj1AquMJeQmHbw8IOK87OUxVQND/OzlOVsLDDOl23uOIoiqL4+/sr48aNK/LxJ0+e\nVIYNG6Y0aNBAsbGxUbp06aKkpqYq8+bNU8zMzBQzMzOlevXqSo0aNRS1Wq0oiqJs3bpVcXd3V6ys\nrJTmzZsr69ev18RbtWqV4uLioixdulRp2rSpUrduXeWNN95QcnJyFEVRlMzMTOWtt95S2rRpo9So\nUUNp3LixsnnzZp1fp6IU/rkIIcTjFPX7o9L12JYs2Uds7Jx822Jj5zB69Azc3Ire2zpzZh+pqQXj\nhITM0LnXZmRkxA8//MCWLVuwtbVlwIABzJw5s9C16V599VU6d+5MTEwMFhYWHD9+HHNzc95++22i\no6MxMzPL1/vbv38/o0aNIjQ0lN69e3PgwAFGjhyJo6Mjnp6eAFy9epUzZ85w+vRpEhMT6dmzJ7Vr\n1+a9995jzZo1bN26lYMHD+Lo6EhiYiJpaWk6vUYhhCgtle4ZW1aW9lyuVpvoFEet1h4nM1O3OABB\nQUH88ccf3Lhxg0WLFrFr167/fE7m4ODAjRs3+PPPPzExMaFDhw5UrVoVQLNu0cMWL16Mv78/L7zw\nAtWrV2fgwIH06dOHNWvWaI5RqVS8/fbbVKlShUaNGjFixAhCQ0M118vOziYmJgaVSoW9vb2suCCE\nKDVXrlzhhRdeKPLxla7HVqWKSuv2Tp1y2bOn6HF69VKxb1/B7RYWuTq3yd3dXfP7s88+S1paGn5+\nfqxevRozM7MCx3/22WcsXrwYHx8fzMzMGDVqFB9++CGQ1/szMjLKd3xcXBwRERF89913mm25ubl4\nef3bs7SxscHJyUnzd/v27Vm4cCEAvXv35v3332fGjBlcuHCB//3vfyxYsICnn35a59cqhBC6Wrdu\nHVu2bCny8ZWuxzZpUk+cnafl2+bsPJWgoB5lEkebB4np0Z7XA/b29nz66ackJCSwfft2QkJC2PNP\nVjYyMkKtVuc73snJiTFjxpCamqr5SUtLIywsTHNMSkoK8fHxmr+PHTuGo6OjJubYsWM5cuQIV69e\nBeC9994z+HUKIcTjKIrCmjVr6N69e5HPqXQ9tgfPv0JCZpCZaYKFRS5BQb11fi5WXHEANm7cSJ8+\nfahatSoRERHMmTOHAQMGYG5urvX41atX06NHD+rXr4+xsTEmJibcunULgPr167N//34yMzOxsLAA\nIDg4mCFDhtCtWzd8fX2pWrUqZ8+eBaBdu3YAmJqaMn/+fBYuXMiff/7Jxo0bCQwMBODQoUNYWVnh\n5uZGbm4u5ubmJCUl6fw6hRBCVzk5OUycOJGmTZty8ODBop1UggNYSkVhL6EivTRvb2/FxsZGqVq1\nqvL0008rb7zxhpKenl7o8f7+/kr9+vWVGjVqKF27dlVmz56t2XflyhWlY8eOipWVlWJtba0ZFRke\nHq506tRJqVWrllK7dm2lW7duyuHDeSM4H4yKXLZsmdKkSROlTp06ypQpUxSVSqUoiqJs2LBBadGi\nhVK9enXFxcVFGTdunHL58mW9XmtF+lyEEOVLUb8/pLq/IDQ0lDlz5nDp0qUSv5Z8LkIIfUl1fyGE\nEJWSJDahdSSlEEJUVHIrUpQq+VyEEEWRlZWFqakpJib/zg0uV7ciAwMDqVevHq1atSqwLzc3l7Zt\n29KvXz/NtsjISNzd3XFzcyMkJKQ0miiEEKIc+frrr3FycuL27ds6n1sqiW3MmDGaeVaPWrx4MS1a\ntNDcCsvNzSUwMJAtW7Zw/PhxVqxYQXR0dGk0UwghRDmxZs0a6tSpQ+3atXU+t1QSm6enJ9bW1gW2\nX79+nV27djFu3DhN9zIqKgoXFxecnJwwMzNj2LBhbN++vTSaKYQQohw4f/48x48fx8/PT6/zy3Tw\nyJQpU5g/fz7Gxv82IzExUVPxAvLqFCYmJpZF84QQQpSBtWvXYmJiwvDhw/U6v8wqj4SFhVG3bl3a\ntm1LRESEQbFmzZql+d3b2xtvb2+D4gkhhCgbubm5rFu3jt69exMdHc3SpUt1jlFmie3IkSPs2LGD\nXbt2kZmZqSn8+8orr5CQkKA5LiEhAQcHh/+M9XBiE0IIUXElJSXh5OSEn59fgY7KBx98UKQYZXYr\n8uOPPyYhIYG4uDg2btyIj48Pa9asoX379ly6dIn4+Hiys7PZtGkT/fv3L6tmloqNGzfi6elJzZo1\ntVbzB9i2bRutW7fGxsYGHx8fvauEhIaG0rhxY0OaK4QQJaZ+/fr8/PPPDBkyRO8YpZLYhg8fTpcu\nXYiJicHR0ZFVq1YVOObBqEhTU1NWrlzJoEGDaNeuHYGBgU/82l82Nja89tprfP7551r3nzhxgtGj\nR7Nw4UJu3rxJnz596N69O9nZ2aXcUiGEKB0GFY0o/jKVpauwl/BfLy1sX5jSM6Cn0s2/m9IzoKcS\nti9Mr2sXV5wHDh06pJiamhbYPm7cOMXPzy/ftoYNGyobNmzQGiclJUWZOHGi0qxZM6VmzZqKq6ur\n8tNPPylHjhxRLCwsFGNjY6V69epK9erVNYWQz549q/Ts2VOpXbu24ujoqLz33ntKTk6OoiiKEhcX\npxgZGSmbNm1S3N3dFWtra2XAgAFKUlKS5prLli1TvL29lZo1ayr29vbK1KlTtbbtCfhfTghRRor6\n/VHplq0J3x/O5C8nE9s2VrMt9su83/v26FvqcYoiOjqaoUOH5tvWpk0bzp8/r/X4BQsWcPbsWQ4f\nPkzdunW5fPkypqamODk5sWzZMmbPnp3vVmZSUhJeXl5MmTKFdevWERcXR2BgIJaWlsyYMUNz3Bdf\nfMHu3buxtLRk5MiRjBo1ir179xITE0NQUBCHDh2iS5cupKeny9xDIUSZqXSJbcl3S/IlI4DYtrGM\nXjAat0S3Isc5s+EMqV1SC8QJ2RBS7IktOTmZmjVr5ttWq1YtzRpsj3JwcODu3bvExsZia2uLi4uL\nZp+ipRzNmjVrqF+/viaJ1alThwkTJrB48eJ8iS0wMJC6desCeWu8+fr68tdff2Fra4u5uTmXL1+m\nVatWWFlZ0bFjR4NftxBC6KPSFUHOUrK0blcbqbVuL0xhx2eqM3Vu0+PY2try999/59uWmpqqSTKP\nGjNmDP7+/kyYMIG6desSEBDwn2Vp4uLiiI2NxdraWvMzffp0bt68me84Dw8Pze/t27cH8ibZ29jY\nsGPHDrZu3YqjoyOenp5FXxBQCCGApUuXEhQURG5ursGxKl2PrYpRFa3bOzXoxJ4A7WW/tOl1uBf7\n2Fdgu4Wxhd5tK0yLFi04ceKE5m9FUTh16hQjRozQeryFhQXBwcEEBwdz8+ZNBg8ezPz58/nkk08w\nNjYu0GtzcnKiSZMmnDlz5j/bERUVhaurKwDHjh0D0EzF8PHxwcfHB5VKxcKFC/H39+fq1av5Jt8L\nIYQ2iqLwxRdfYGNjk6/osb4q3bfOpBGTcD7pnG+b8wlngoYHlUkcALVaTWZmpmaUY1ZWFpmZ//b8\nJk6cyNatWzl48CBZWVnMnz8ftVrN4MGDtcYLCwsjOjqa3NxcFEXBzMxMc9vSzs6Oa9eu5euN+fn5\ncePGDT755BOuXr2KWq3mypUr7N27N1/c0NBQkpKSSEtLIyQkhB49emBnZ0dMTAx79uwhIyMDlUqF\nubk5qampqFQqnd8LIUTlc+LECS5cuMDo0aOLJ2DJjV8pHYW9hP96aWH7wpReY3op3fy7Kb3G9DJo\nVGRxxFm1apViZGSkGBkZKcbGxpr/Xr16VXPMtm3bFDc3N8Xa2lrx8fFRYmJiCo33+eefKy4uLkr1\n6tUVV1dX5bXXXlNu3bqlKIqi5OTkKM8//7xSq1YtpVatWkpkZKSiKIpy4cIFpX///kq9evWUmjVr\nKq1bt1aWLl2qKMq/oyK///57zajI/v37Kzdv3lQUJW9EZZcuXZSaNWsq9evXV1544QVl//79Wtv2\nBPwvJ4QoZpMmTVKqVKmipKSk/OdxRf3+kPXYxGPFx8fTqFEjrl+/ToMGDQyKJZ+LEOJhOTk5NGjQ\ngO7du/P999//57Hlaj02IYQQQpvff/+dlJQUvSv5a1PpBo8I/RhUBUAIIQrRuXNnrl+/jq2tbbHF\nlFuRolTJ5yKE0JfcihRCCFEpSWITQgjxRJHEJoQQ4okiiU0IIUSpmzFjBkePHi2R2JLYhBBClKpT\np04xe/Zsjhw5UiLxJbEJIYQoVV999RWWlpYEBASUSHxJbEIIIUrN33//zfr16xkxYgTW1tYlcg1J\nbOXAO++8Q8uWLbGyssLe3p6XXnqJ1NT8a71t27aN1q1bY2Njg4+PT76FQnURGhpK48aNi6PZQgih\ns9WrV5ORkcErr7xSYteolIktMjyc6b16Mcvbm+m9ehEZHl6mcUxNTVm/fj2pqans27ePc+fO5eui\nnzhxgtGjR7Nw4UJu3rxJnz596N69u2Y1ACGEqCjWrl1Lp06dcHd3L7mL6FuNubwo7CUUtv1wWJgy\n1dlZUUDzM9XZWTkcpltl/uKKo82mTZsUKysrzd/jxo1T/Pz88h3TsGFDZcOGDVrPT0lJUSZOnKg0\na9ZMqVmzpuLq6qr89NNPypEjRxQLCwvF2NhYqV69ulK9enXl8OHDiqLkVejv2bOnUrt2bcXR0VF5\n7733lJycHEVR/q3uv2nTJk11/wEDBihJSUmaay5btkzx9vZWatasqdjb2ytTp07V2rYn4H85IYQB\n7ty5o1y8eFGvc4v6/VHpakXuW7KEObGx+bbNiY1lxujReLm5FT3OmTPMeeR24ZzYWGaEhODVt69B\nbfzpp5/y/WsmOjqaoUOH5jumTZs2nD9/Xuv5CxYs4OzZsxw+fJi6dety+fJlTE1NcXJyYtmyZcye\nPTvfrcykpCS8vLyYMmUK69atIy4ujsDAQCwtLZkxY4bmuC+++ILdu3djaWnJyJEjGTVqFHv37iUm\nJoagoCAOHTpEly5dSE9PJzo62qD3QAjxZLKyssLKyqpEr1EqtyIDAwOpV68erVq10mxLSEige/fu\nuLq64u3tTWhoqGZfZGQk7u7uuLm5ERISUqxtMc3K0rrdRK3WLU4hx5s8tECoPjZv3kxoaChffvml\nZltycjI1a9bMd1ytWrU0i4c+ysHBgbt37xIbG4tarcbFxQUnJycArXXW1qxZQ/369ZkxYwZ16tSh\nQ4cOTJgwgTVr1uQ7LjAwkLp161KjRg2Cg4PZv38/f/31F7a2tpibm3P58mXS09OxsrKiY8eOBr0P\nQgihr1LpsY0ZM4agoKB8yxKYmZmxaNEi2rRpQ3JyMi1btqRTp040btyYwMBADhw4gL29PR4eHvj6\n+tK8efNiaYuqShWt23M7dYI9e4oep1cv2LevYBwLC73btnnzZiZMmEBYWBgtWrTQbLe1teXvv//O\nd2xqaiqtW7fWGmfMmDFkZWUxYcIEEhMTee6551i4cCG1a9fWenxcXByxsbH5RigpioL6keTt4eGh\n+b19+/YAXL9+nfbt27Njxw6WLFnCpEmTcHNz44MPPsDHx0e3N0AIIYpBqfTYPD09CwzrtLOzo02b\nNkDeF7eHhweJiYlERUVpehhmZmYMGzaM7du3F1tbek6axDRn53zbpjo70yMoqEziPBAaGqpJat26\ndcu3r0WLFpw4cULzt6IonDp1CldXV62xLCwsCA4O5vTp05w/f55Lly4xf/58AIyNjQv02pycnGjS\npAmpqaman7///pu0tLR8x0VFRWl+P3bsGJDXOwTw8fFh27Zt3L59m+eeew5/f/8CiVEIIUpDuRgV\nefnyZc6fP0+nTp1ITEzE0dFRs8/BwYHExMRiu5ZX3770WryYGb16MatbN2b06kXvxYt1fi5WXHEA\nQkJCePPNN9m3bx+dO3cusH/ixIls3bqVgwcPkpWVxfz581Gr1QwePFhrvLCwMKKjo8nNzUVRFMzM\nzDS3Le3s7Lh27Ro3b97UHO/n58eNGzf45JNPuHr1Kmq1mitXrrB37958cUNDQ0lKSiItLY2QkBB6\n9OiBnZ0dMTEx7Nmzh4yMDFQqFebm5qSmpqJSqXR+L4QQT56DBw+yefNmcnNzS+V6ZT545O7duwwb\nNoxFixZRrVo1vRa0nDVrluZ3b29vvL29//N4r759DR7gUZxxJk+ejJmZWb52GxkZaXpMbdu2Ze3a\ntUyZMoWEhATatm3LwYMHMTMz0xovNjaWKVOm8Ndff9GwYUO6d+/OzJkzgbyeVf/+/WnWrBkAO3bs\nwNPTk8jISN59910+//xzMjMzcXJyYsKECfnivvbaa/Tp04e4uDg8PT1Zu3YtANnZ2Xz00UecP3+e\nqlWr0rVrV7Zt24a5ubnB740QouKbNWsWf/75J88//7xO50VERBAREaHz9UptodH4+Hj69evH2bNn\nNdtycnJ47rnn6NOnD8HBwQD89ttvzJo1iz3/PO+aO3cuxsbGvPPOO1rjykKjJS8+Pp5GjRpx/fp1\nGjRoYFAs+VyEqFzOnj2Lm5sbCxYs4I033jAoVrlfaFRRFMaOHYurq6smqUHeoIRLly4RHx9PdnY2\nmzZton///mXVTCGEEAb46quvsLCwYMyYMaV2zVK5FTl8+HAOHz5McnIyjo6OfPjhhzRp0oR169bh\n5uZG27ZtgbzeWe/evVm5ciWDBg1CpVIxfvz4YhsRKfSnzy1iIUTldufOHdauXcvw4cOxsbEpteuW\n2q3IkiK3IisW+VyEqDyWLVvGxIkTOXbsmGaKkCGK+v0hiU2UKvlchKg8srOz+fHHH+nTp0+xxJPE\nJl+g5ZJ8LkIIfZX7wSNCCCFESZDEJoQQ4okiiU0IIcQTRRKbEEKIYnPu3DlCQkK4f/9+scUMP3iQ\nXpMmFfl4SWzliFqtpkuXLhgbG/Pnn3/m27dt2zZat26NjY0NPj4++dZT00VoaCiNGzcujuYKIUQB\nH3/8MVOnTiXTwCW8Hgg/eJDJGzawr5DauNpIYitHCquXeeLECUaPHs3ChQu5efMmffr0oXv37mRn\nZ5dRS4UQoqDY2Fg2bdrExIkTC6zooq8l27YRO3KkTudUysT2oFvrPXkyvSZNIvzgwTKNAxATE8PS\npUtZsGBBgeGsS5cuZfDgwfj6+mJmZsZbb72FqakpW7Zs0RorNTWVV155hebNm1OrVi1atmzJzz//\nzK+//srEiRO5cuUKNWrUoEaNGkRGRgJ5tw969eqFra0tTz31FFOnTtVU54+Pj8fY2Jjvv/+edu3a\nYWNjw8CBA/MtdPr111/TvXt3atWqhYODA9OmTdP7vRBCVEzz58/H1NQ0X5lEQ2XpUfWozKv7l7YH\n3dqH/wUQu349AH11WBizuOJA3i3IwMBAFi5cWGClbIDo6GiGDh2ab1ubNm04f/681ngLFizg7Nmz\nHD58mLp163L58mVMTU1xcnJi2bJlzJ49O9+tzKSkJLy8vJgyZQrr1q0jLi6OwMBALC0tmTFjhua4\nL774gt27d2NpacnIkSMZNWoUe/fuJSYmhqCgIA4dOkSXLl1IT08nOjpap/dACFGx3bhxg1WrVhEQ\nEGBwsfSHqfVY6qbSJTZt3drYkSMZvXo1bjp0nc+sXUuqv3+BOCFbt+qc2BYvXkyDBg0YMGAA8fHx\nBfYnJycXSHi1atXK12N6mIODA3fv3iU2NhZbW1tcXFw0+7RNblyzZg3169fXJLE6deowYcIEFi9e\nnC+xBQYGUrduXQCCg4Px9fXlr7/+wtbWFnNzcy5fvkyrVq2wsrKiY8eOOr0HQoiKbePGjahUKt56\n661ijWvZsSN8+y2MG1fkcyrdrcjCurVqY93eisKO1/Vx6eXLl/nss88ICQnJt/3hBGRra8vff/+d\nb39qaqomyTxqzJgx+Pv7M2HCBOrWrUtAQAC3b98utA1xcXHExsZibW2t+Zk+fXq+xUgBPDw8NL8/\nqPt2/fp1bGxs2LFjB1u3bsXR0RFPT08OGnBbVghR8QQHB3Py5Ml8/5A21KWMDH50dOS53r3ptXVr\nkc+rdD22KoWUY+lUvTp7/llloCh6VavGPi3bLXRsz88//8ytW7do2bIlkHdbEsDNzY05c+YwYcIE\nWrRowYkTJzTnKIrCqVOnGDFihNaYFhYWBAcHExwczM2bNxk8eDDz58/nk08+wdjYuECvzcnJiSZN\nmnDmzJn/bGtUVBSurq4AHDt2DMjrHULeAqY+Pj6oVCoWLlyIv78/V69exVjHfzAIISomIyMj3Nzc\nijXmB/HxmBsbs3zIEOxGjcJoyZIinVfpvnUmDRyI8z/Pwh5wXreOoAEDyiTOiy++yJUrVzh9+jSn\nT59m165dAOzfv5/Ro0cDMHHiRLZu3crBgwfJyspi/vz5qNVqBhcy/DUsLIzo6Ghyc3NRFAUzMzPN\nbUs7OzuuXbuWrzfm5+fHjRs3+OSTT7h69SpqtZorV66wd+/efHFDQ0NJSkoiLS2NkJAQevTogZ2d\nHTExMezZs4eMjAxUKhXm5uakpqZqBp8IIYSuzt29y3dJSUyyt8euShWdzq10PbYHz79Ctm4lk7we\nVtCIETo/FyuuOJaWllhaWmr+zs7OxsjICDs7O6pVqwZA27ZtWbt2LVOmTCEhIYG2bdty8OBBzMzM\ntMaMjY1lypQp/PXXXzRs2JDu3bszc+ZMIK9n1b9/f5o1awbAjh078PT0JDIyknfffZfPP/+czMxM\nnJycmDBhQr64r732Gn369CEuLg5PT0/Wrl2rafNHH33E+fPnqVq1Kl27dmXbtm2Ym5vr9F4IIcQD\nM+PjqWFiwltPPaXzuVLdXzxWfHw8jRo14vr16waPdpLPRQjxOMfT02l//DiznJyY6eSk2S7V/YUQ\nQpSoH374genTpxdblZEHpsfFYWNqypR/nuHrShKbKJJHq6EIISo3tVrNzJkz2b59e7E+dvj577/Z\nk5LCO089hZWpfk/LKt0zNqE7JycncvWYJCmEeHKFh4dz/vx51q5dW2yjnxVFYXpcHPXMzHjN3l7v\nOPKMTZQq+VyEqPgURaFr167cuHGDS5cuYapnz+pRB1JS6HHmDEtcXAjSchuyXD1jCwwMpF69erRq\n1Srf9sjISNzd3XFzc8s3Qbmw7UIIIcpeZGQkv/76q6ZubXFQFIVpcXE4VqnCSwYOUiuVxDZmzBj2\n7NmTb1tubi6BgYFs2bKF48ePs2LFCs3cK23bhRBClA8RERHUq1ePMWPGFFvMnbdvE5WezvsNG1LF\nwFubpZLYPD09CyxhEBUVhYuLC05OTpiZmTFs2DC2b9/OsWPHtG4XQghRPsycOZPz58/nm4NrCLWi\nMCMuDhdLS/zt7AyOV2aDRxITE3F0dNT87eDgwNGjRwvdLoQQovyoXbu2wTHCDx5kybZtJKhURKen\n88agQZgVQwH1MktsxTl8fNasWZrfvb298fb2LrbYQgghip+2pb+2rV9P91q1NBWcIiIiiIiI0Dl2\nmSU2e3t7EhISNH8nJCTg4OBQ6Pb/8nBiE0IIUf4VtoTYw0t/PdpR+eCDD4oUu8wmaLdv355Lly4R\nHx9PdnY2mzZton///oVuf1J5e3szZ84czd+pqal4eXnRrVs37ty5U4YtE0KIf92/f79Y4xW2hFhx\n1DAplcQ2fPhwunTpQkxMDI6OjqxatQpTU1NWrlzJoEGDaNeuHYGBgTRv3rzQ7U8qIyMjzW3ZhIQE\nPD09qVevHvv379e6mrYQQpS2u3fv0qRJE7744otii1nYEmK6Lv2llVLBFfYSKspL8/b2VubMmaOc\nPXtWsbe3V4KCgjT7Vq1apbi4uCirVq1SmjVrplhbWysvv/yykpubqznm9OnTSvfu3RVra2ulUaNG\nyuzZszXQt5pZAAAgAElEQVT7g4KClJdeeklzrKenp9KwYUPN3/PmzVOeffZZRVEUZebMmYqvr68y\nd+5cxcnJSalbt64yc+bMYn+9FeVzEUL868MPP1QA5ddffy22mN/s2qUwcqTCoUOaH+exY5WwH38s\n9Jyifn9U6pJahQ0yKexhpa7HF9WRI0dYsGABb7/9Nu+++26+fVevXuXo0aMcOnSIK1eu4Ovri5eX\nFyNGjODOnTv06NGDV155hb179xIbG0uvXr2oUqUKb775Jj169CA4OBjI+xfXyZMnqV27NpcuXaJx\n48bs37+fvn37aq4VGRlJly5dOH78OJGRkQwePJiePXvSpUsXg16fEKLiSk5OZv78+QwcOJBOnToV\nW9ydDg5YdOxIpx9+QDE21nvpL20qdWIrDxRF4eeff8bS0pLhw4dr3T9nzhxsbGyws7PjmWee4fff\nf2fEiBGEh4dz584d3n77bczMzGjWrBkjR47k22+/5c0336Rbt24kJCQQFxfHhQsX6NChgyahPfXU\nUxw5coRFixZprlWzZk1mzpyJsbExAwcOxMXFhePHj0tiE6IS+/jjj7l37x4ff/xxscXcc/s2O2/f\nZt6AAbytx3prj1OpE5uuPS1De2baGBkZ8cYbbxAbG4unpycHDhygSZMmmv1OTk7Y2Nho/ra3t+fu\n3btA3jM5V1fXfJMkO3bsyOLFiwGwsrLCw8ODAwcOEB0dTc+ePXF2dmb9+vU0bdqUGjVq0LJlS825\nrVu3zlfM1N7envT09GJ/zUKIiuHatWt8+eWXBAQEFNtYh2y1muDLl3GxtGSynsvSPI4sW1MOmJiY\nEBoaSr9+/fDy8uLs2bNFOs/R0ZELFy6QkZGh2fbbb7/x1EP/AvL19WX//v0cOHCAHj164OPjw+HD\nh9m3bx++vr7F/lqEEE8OMzMz/Pz8inVK1ReJifxx/z6LnJ0NLp1VGEls5YDyz+igL7/8Ej8/P7y9\nvTl27Fihxz44vm/fvlhZWbFw4UJycnL4448/2LhxI2PHjtUc7+vry549e/jrr79wd3fHxsaGp59+\nmq+//rpIiU2RSvxCVFr169dn+fLl+apBGeJmdjYfxMfT28aGvsVQuaQwktjKgYersHz66adMnjwZ\nX19fEhISClRoeXh6QM2aNdm3bx+HDh2iXr169O7dm8DAQKZMmaI5vlOnTiiKgs9DD2R9fX1JT0/P\nl9gejltY24QQwhDTrlwhQ63mcxeXEv1ukfXYRKmSz0WIyul4ejoex4/zuoMDC1xc9IpRrtZjE0II\nUXkpikLQpUvUMTNjhpNTiV9PEpsQQgiN3377jZMnTxZrzPU3b/JrWhpzGzWiZjEtTPpf5FakKFXy\nuQhRfuXm5uLu7s79+/e5ePFivuk/+rqrUtE0KooGVapw1N0dYwOerRX1+6NSz2MTQgjxr2+++YYz\nZ86wadMmg5Pag7XWorOy+DMjgylDh2Lcrl0xtfS/SY9NlCr5XIQon5KSkmjatCnu7u4cOHDAoFGL\n2tZac16/nsXDhxtUMksGjwghhCiyd955h3v37vHll18aPBS/0LXWtm/XK15keDjTe/Uq8vGS2IQQ\nopJLTk5m+/btvPHGGzRr1szgeMW51lpkeDh7J09m9r59RT7niX3GZm1tLZOLyyFra+uyboIQ4hG2\ntrZcvHiRatWqFUs8Y7Va63Z91lrbt2QJc2JjdTrniU1sKSkpZd0EIYSoMOrWrVtssYw9PODbb2Hc\nOM0253XrCBoxQudYpllZup+j8xlCCCFEIfbcvs2PTz3F0Gef5c7WrWSCQWutqVJTdT7niR0VKYQQ\nonSlqVS0PHaM6iYmnGzf3vDq/VeuENmyJXvVauZkZWFE0Qqzy+ARIYSohBITE1EX8ixMX+9cucL1\nrCxWNm1qeFLLzoYXX8SrShV6ffUVM3QYFSk9NiGEqGSys7Np3bo1bdq0YcOGDcUSMyI1le6nT/O6\ngwML9SxynM8bb8Bnn8EPP8DgwYDMYxNCCFGIRYsWcfHiRUaNGlUs8e7l5jL2jz9wtrDgo6efNjxg\neHheUnvlFU1S00WZJ7bly5fTpUsX2rVrR3BwsGZ7ZGQk7u7uuLm5ERISUoYtFEKIJ8e1a9f48MMP\nGTBgAH379i2WmDPi4riSmcmKZs2oamJiWLDERPD3h9atYeFCvUKU6a3IlJQU2rVrx7lz57C0tOS5\n554jODiY//3vfzRt2pQDBw5gb2+Ph4cHGzZsoHnz5gViyK1IIYQouueff57du3cTHR1Nw4YNDY73\n6507dD15kokNGvBlkyaGBcvNhf/9D37/HY4fh6ZN8+0utiLIAQEB1KlTh65du9K5c2fq1aunf6Mf\nYWlpiaIo3LlzB4CMjAxq1apFVFQULi4uOP2zbs+wYcPYvn271sQmhBCiaA4fPsyWLVuYM2dOsSS1\nzNxcAv/4A8cqVfikUSPDG/jRR3D4MKxeXSCp6eKxtyJDQ0MJDAwkNTWV999/nw4dOvDpp58Wy2ga\nS0tLli5dipOTE3Z2dnTt2pUOHTqQmJiIo6Oj5jgHBwcSExMNvp4QQlRmzzzzDKtWreKNN94olngf\nXr3KxYwMljdtSg1D11mLiMhLbH5+eT8GeGxLTp48SWZmJmPGjGHMmDHs3r0bZ2dnVq1axdixYw26\n+K1bt5g4cSIXLlzA2tqaIUOGEB4ernMprFmzZml+9/b2xtvb26B2CSHEk8jExISAgACDYjxYjiZZ\nrebEnTv49uhBT0O/c2/dghEjoHFj+PJLzeaIiAgiIiJ0DvfYxLZnzx6MjY35/PPPqVq1Kvb29lhY\nWGBubq7zxR4VFRVFp06dcPlnaOiQIUOIjIxk0KBBJCQkaI5LSEjAwcGh0DgPJzYhhBAlQ9tyNLHr\n1hHu4KD/cjRqdd5gkZQU2L0bqlfX7Hq0o/LBBx8UKeRjb0UOHDiQ7t27s2nTJlatWsXs2bO5evUq\nNjY2ur+AR3h6evL777+TkpJCVlYWu3fvpmfPnnh4eHDp0iXi4+PJzs5m06ZN9O/f3+DrCSGE0J+2\n5WjiRo3SezkaIG9Y/+7def9t3drAFuZ5bI9N24ANQ7uyD1hZWTF9+nQGDRpERkYGvXv3pnv37hgb\nG7Ny5UoGDRqESqVi/PjxMnBECCH0cPfuXao/1AsyRHEuRwPA0aPw3nswaBBMnKh3ux4llUeEEOIJ\ndfjwYQYPHszu3bvp0KGDwfG6BwUR8fzzBbb32rqVPYsX6xbs77+hbVtQFDh5EoqwpJVUHhFCiEos\nLS2NgIAAbGxscHV1NTieWlFIb9MmbzmahzivW0fQgAG6BVMUGD8eEhJgw4YiJTVdyLI1QgjxBHr9\n9de5du0aP/30U7EsILro+nWOOzszwcqKOEOXo/n6a/i//4NPPoHOnQ1u26PkVqQQQjxhdu7cSf/+\n/Xn33XeZO3euwfGi0tLoevIk/WrX5gdXV52nZOVz5gx06ADe3rBrF+iwCkBRv+8lsQkhxBMkOzsb\nZ2dnbGxsiIqKokqVKgbFu6NS0fb338lVFE61b4+1mZn+we7dg/bt856vnT4NOq7aXWwltYQQQlQc\n5ubmbN26FQsLC4OTmqIojP/jD65lZhLZtq1hSQ0gKAj++AP279c5qelCEpsQQjxh2rdvXyxxvrlx\ng823bjH36afpUrOmYcHWr4dVq2D69LxCxyVIbkUKIYQo4Ozdu3Q4cQKvmjXZ7eaGsSHP1S5dAnf3\nvOH9Bw+CnnUl5RmbEEIIvdzLzcXj+HFSVSpOtW9PPUNKKGZl5Y18vHoVTp2Chwrc60qesQkhRCWR\nmJiIvb29wXEeFDg+ff8+NzMz+eiFF/ROapHh4exbsgTTc+dQ/fknPadPx8uApKYLSWxCCFGBnT9/\nno4dOzJv3jxeffVVveNoK3Acun49bWvU0HmeWmR4OHsnT2ZObKxm27QNG6BTJ7yKadXu/yKVR4QQ\nooJKS0vj+eefp1q1agwaNMigWNoKHMeOHKlXgeN9S5bkS2oAc2Jj2R8SYlAbi0p6bEIIUQEpikJg\nYCCXL1/mwIEDNGjQwKB4dwvZrk+BY9M7d7RuN8nUu1yybtcvlasIIYQoVosWLeKHH37g008/NXhx\n5Wy1muj0dK37LHQNlpSE6uxZrbtyLXSOphe5FSmEEBVMVlYWX3/9NYMHD+bNN980KJaiKLx26RKp\nbdpQd/XqfPt0LnCclgZ9+tBTpWLaIz3Iqc7O9AgKMqitRSXD/YUQogJKTU3FxMQEKysrg+IsuX6d\nyZcv895TT9E1Pp6Q7dv/LXA8YEDRB45kZkKfPvDzz7B9O5GKwv6QEEwyM8m1sKBHUJDBA0dkHpsQ\nQoj/tDclhWfPnKFf7dpsadlS/0nYKhUMGQLbtsG6dfDIIJTiIuuxCSGEKNTFe/d48fx5Wlarxrrm\nzfVPaooCL72Ul9QWLy6xpKYLSWxCCFHJpOTk0O/cOcyNjdnRqhXV9SxxBcA77+TVgHz/fZg0qfga\naQBJbEIIUc5FR0fz+uuvk52dbXCsHLWaoefPczUzk60tW9LQkJGKn34K8+fDq6/CrFkGt624yDM2\nIYQox9LT0+nQoQMpKSmcPHlSr/lqD0plZRkZceXePRJatWLV0KEE1K+vf8O+/RbGj4fhw/Oeq+mw\nYKi+pFakEEJUcIqiMHbsWGJiYvSehK2tVFbN0FDqtGoF+ia2H36Al1+G3r0hNLRUkpouyrw19+7d\nw9/fn7Zt29KiRQuOHj0KQGRkJO7u7ri5uRFSSmVYhBCiPFmwYAGbN29m7ty5dO/eXa8Y2kpl3QkI\n0KtUFgA//ggjRkDHjvB//weGVP4vIWXeY3vllVfo1q0bq1evRqVSce/ePXJzcwkMDOTAgQPY29vj\n4eGBr68vzZs3L+vmCiFEqdi1axdvv/02Q4YM4a233tI7TlYhox31Km517BgMHAhNmkBYGFSrpne7\nSlKZ9tju3LnDTz/9RGBgIACmpqbUrFmTqKgoXFxccHJywszMjGHDhrFd339dCCFEBeTp6cnbb7/N\nmjVrMDJgkc8clUrrdp2HjFy8mDcB29YW9u4FGxu921TSyjSxxcXFUadOHQICAmjZsiXjx4/n/v37\nJCYm4vjQuj0ODg4kJiaWYUuFEKJ01ahRg3nz5mFhwKjFuPv3udiiBSYrVuTbrnOprGvXoEcPMDGB\n/fvBwILLJa1Mb0WqVCqOHTvG9OnTWbp0KS+//DKbN2+mmo7d21kPDTP19vY2uCCoEEJUdDeysvA9\nfRqldWsWOzmxc+vWf0tljRhR9FJZt25Bz555dSAPHwYXl5Jsdj4RERFEREQQExvD0bNHi3xemSY2\nBwcHateuTb9+/QAYPnw4a9asYfLkySQkJGiOS0hIwMHBodA4s8rR/AkhhChrKTk59DxzhpvZ2fzY\npg0drax4VZ86jenp8OyzcPVq3u3HNm2Kv7H/wdvbm3s591h3ah1XBl2B00U7r0xvRdrZ2eHi4sLR\no0dRq9WEh4fj6+uLh4cHly5dIj4+nuzsbDZt2kT//v3LsqlCCFFicnJy+PTTT8kshvXK7qpUPHvm\nDDEZGWxv1YqO+hZJzsrKGyhy8iR8/z14eRncNn0s+W4JsW1jH3/gQ8p8VOTq1avx8/MjOTmZVq1a\nMW/ePExMTFi5ciWDBg1CpVIxfvx4GREphHgiKYrCxIkTWbFiBY0bNzZoJewstZqB585xLD2d/3N1\n5X/W1jqdHxkezr4lSzDNzER18SI9k5LwWrMG/rmrVhZSslJ0PqfME1uTJk347bffCmzv1q0bJ0+e\nLIMWCSFE6Zk9ezYrVqxg+vTpeiW1B1VFMo2MuJieTlLr1qwaOpRBderoFCcyPJy9kyczJ/bf3tE0\nW1uwsaFs+mqw9vRaTiSegKa6nVfmE7SFEKKyCg0N5f3338fPz48PP/xQ5/MfVBXZN3gwkYMGkeTn\nR+1Tp6gTHa1zrH1LluRLagBzkpPZXwYFMjJyMhi7fSx+2/xo3qk5DY831Ol8SWxCCFEGfvnlF8aP\nH4+vry/Lly/Xa66atqoit/399aoqYlrI8z2TYnjup4uLyRfp+G1HVp1axXTP6Zyae4ovg76k19Ve\nRY4hiU0IIcqAu7s7r7/+Oj/88APmepalyiyuqiLp6ajOn9e6K9eQ6v86WndmHe2/ac/NuzfZM2oP\nH/l8hKmxKX179GXPyj1FjiOJTQghyoClpSXz5s3DSs9Ri7mKQuzdu1r36ZSKEhPB05OeKSlMe+S5\n3FRnZ3oEBenVPl1k5GQwbsc4Rm8dTbsG7Tg14RQ9nXvqHa/MB48IIYTQTY5ajd/FiyS6uWEdGkpq\nQIBmn/O6dQSNGFG0QGfO5M1TS0vDa/duUKmYERKCSWYmuRYW9A4Kwkuf+W86uJh8kSGbh3A+6TzT\nPKcxy3sWpsaGpSZZj00IISqQzNxcXrxwgR23bzOvUSNcL18mZPv2f6uKDBhQtKoie/fCkCFgZQXh\n4dC6dUk3vYD1Z9bzctjLWJpZsm7QOnq5/PdztKJ+30tiE0KIEhYVFcXcuXNZt26dziUDH3YvN5eB\n585xIDWVLxs35hV7e/0CffstTJgALVvmJTV94+jpfs59Ju2exLcnv8WroRffDf4Oe6vHt0EWGhVC\niHIgKiqKHj16UKdOHe7cuaN3YrujUtH3zBl+TUsjtFkz/O3sdA+iVsP06TB3bt4iod9/DzVq6NUe\nfV1MvsjQzUM5m3SWqc9M5YPuHxh86/FRktiEEKKEPEhqtra2HDp0SKcVsB9MvM4yMsJYreZay5Zc\nbdaMjS1aMKRuXd0bk5kJY8bAxo3w0kvw5ZdgWrop4OFbj7tH7qa3S+8SuY4kNiGEKAFRUVH07NkT\nW1tbIiIi8i3F9TgPJl7nm6P27bfMtLPTL6ndvp1X9/Hnn2HePHjrLTBgjTdd3c+5z+Q9k1l+YjnP\nPPUMG57fgINV4YXtDSXD/YUQogR89dVX1K5dm0OHDumU1ED7xGvGjeO3gwd1b0hsLHTunLf69aZN\n8PbbpZrU/kj+g47fdmT5ieW898x7HPI/VKJJDaTHJoQQJeKbb77h9u3b1K9fX+dzs4pr4vWvv0L/\n/qAo8OOP0LWrzm0xxHdnv+OlnS9hYWpRorceHyU9NiGEKAHm5uZ6JTWAv7OytG7XaeL1//0f+PhA\nrVp5Ca4Uk9r9nPu8tPMlRm4ZSdv6bTk14VSpJTWQxCaEEOWGoijMvXqV082bU2Xlynz7nNetI2jA\ngKIEgQUL8uaoubvnJbXGjUuoxQX9kfwHnVZ0YvmJ5bzb9d1SufX4KJnHJoQQBlAUha+++orRo0fr\nXR4L8tZSe/mPP1h98yYj6tblhT//5OudO3WbeK1SwaRJsHQpDB0Kq1dDKdZ63HB2Ay+FvUQVkyqs\nHbSWPo37FGt8mccmhBAlLCcnh5deeonQ0FBycnIIDg7WK05ydjaDz5/npzt3+MDJiRkNG2LUogWD\nfH0fe65mcdCMjLzFQZOT8XrnHfj4YzAunZty93PuE7wnmG9OfENXx65sfGFjqffSHiaJTQgh9JCe\nns6QIUPYu3cvM2fOZPLkyUU+9+E5aiqVilhXV1JdXdnQvDnD6tUrchyti4PWqQOenniVUlKLuR3D\nkM1DOHPzDO92fZcPu3+ImYlZqVy7MJLYhBBCRzdu3KBv376cOXOGb7/9lrFjxxb5XG1z1IxXrOBT\nR0edkhoUsjjorVvMCAkp8eLFkP/WY/iIcJ5t/GyJX7MoZPCIEELoKCQkhJiYGHbu3KlTUgPtc9TU\nY8ey/8AB3RqRmYnphQtad5X04qD3c+4zIWwCI7aMwK2eGydfPllukhpIYhNCCJ19+OGHHDt2jD59\ndB8ckVEcc9ROnQIPD1TXr2vdXZKLg8bcjqHzis58ffxr3un6DhH+ETjW1G0CekkrF4ktNzeXtm3b\n0q9fP822yMhI3N3dcXNzIyQkpAxbJ4QQ+ZmamtK8eXOdz7t47x6n7tzRuq9IqSg3N68kVocOkJxM\nz5kzmebsnO+QklwcdOO5jbT7ph0JaQmEDQ/jE99Pyvx5mjbl4hnb4sWLadGiBenp6UBeogsMDOTA\ngQPY29vj4eGBr6+vXv8jCSFEebD+5k1e/uMPTNq3p/6aNdzw89PsK9LioFeugL9/Xr3H55+HZcvw\nsrUFD48SXxw0U5XJlD1TWHZ8GZ0dOrPphU3lrpf2sDJPbNevX2fXrl1MmzaNzz77DMgrHuri4oKT\nkxMAw4YNY/v27ZLYhBCl6urVq7z++ut888031K5dW68Y93NzmXT5Mt/euIFXzZpsGD+ek7/8QsjW\nrf/OURsxovA5aooCK1dCcHDe8P01a2DUKE29R6++fUt0oMil25cYsnkIp2+e5q0ubzHHZ0657KU9\nrMwT25QpU5g/fz5paWmabYmJifmKhjo4OHD06NGyaJ4QopLau3cvI0aMQKVSER0dzTPPPKNzjD8y\nMhhy/jxn791j6lNP8YGTE6bGxjTw8SnaKtdJSXlLzGzfDt7eeROun3pK9xejp03nNjFu5zjMTczZ\nOXwnzzV5rtSubYgyTWxhYWHUrVuXtm3bEhERoXecWbNmaX739vbG29vb4LYJISontVrNxx9/zPvv\nv4+rqytbtmyhcRFLUj08Py0lK4tLrq5Ua9uW3a1a0VvXHt+OHTB+PNy5A599BpMnl/iE6/D94Sz5\nbgn31feJux3H9TrX6fxMZza+sJGnapZeQn0gIiJCr9xQpontyJEj7Nixg127dpGZmUlaWhp+fn68\n8sorJCQkaI5LSEjAwaHwWewPJzYhhNBXbm4ugwcPZseOHYwcOZKvv/66yCtea5ufZrFyJQsaN9Yt\nqaWnw+uvw7ffQuvWeVX5W7bU9aXoLHx/OJO/nExs23/mxTWCWr/U4h37d8okqUHBjsoHH3xQpPPK\nTa3Iw4cPs2DBAnbu3IlKpaJp06b8+OOPNGjQgA4dOrBhwwatz9ikVqQQojhNnTqVBg0a8Oqrr2Kk\nw7plvSZNYt/gwQW3b93KnsWLixbkl19g9GiIj4d33oEPPgBz8yK3wRDd/boT4RxRYHuvq73Ys3JP\nqbThcSpkrcgH/xOZmpqycuVKBg0ahEqlYvz48TJwRAhRKj7++GOdz8lSq4kuZKmZIs1Py86GmTPh\n00+hYUOIjAQ9nunp49a9Wyz8dSGRCZHgXHB/prpkJ3uXhHKT2Lp160a3bt3y/X3y5MkybJEQQjze\nsbQ0Ai5eJCEjQ+t+bfPTNIWLs7JQ5eTQ88YNvOLiYOxYWLQIatQo2Ubzb0L7IuoLMnIyqGdZj7/4\nq2D7jUtvdYDiUm4SmxBClKatW7dSp04dvUY7Ql4v7cP4eOZdu4aduTkzX3iBdevX53vGpm1+mtbC\nxcbGMH06Xh99pN+L0UFyRjILjizQJLRhLYcxw2sGV7pcyf+MDXA+4UzQayUz2bskSWITQlQqd+7c\nYfLkyaxevZpBgwYVObE9POIxS6XihpsbV5s2ZYydHZ85O1PLzAyPGjUeOz9Na+FitZoZx47hVUyv\nUZvkjGQWHllISFRIvoTWvE7eY57mPfL+G7IhhEx1JhbGFgS9FkTfHiVfTLm4SWITQlQahw4dIiAg\ngOvXrzN9+nRmzJhRpPO0jXg0WbGCmXZ2zHpo1F7fx81P+/NPTM+c0bqrpAoXPy6hPaxvj74VMpE9\nqlzUihRCiJL20Ucf4ePjQ5UqVfjll1/46KOPMC/iiENtFflzx47lt4MHi3bxu3dh1ixo3BjVXwWf\nY0HxFy5OzkjmvQPv4fS5E/N+mUe/pv0498o5vnv+O61J7UkiiU0IUSm0adOGV199lZMnT9KpU6ci\nn3cjK4tT9+9r3ffYPlZuLqxYAU2a5A3d79uXnsuXl2jh4uSMZKb+OJWnFz+dL6FteH4DLeq0KJZr\nlLbw8Eh69Zpe5OPlVqQQolLo169fvhVEHidHrWZJYiKz4uO5V8htwv/sY+3dC2++CefOQefO8MMP\n0Llz3nO0+vWLvXBxckYyn/36GSFRIdzLvseLLV9khteMCpvMHggPj2Ty5L3Exs4B5hTpnHIzQVtf\nMkFbCPEwRVHIzc3F1FT/f7f/mJpK0KVLRGdk0NfGhgGJiczbsqXAiMfF2ooXnzkDb70F+/ZBo0Z5\ny8w8/7ymaHFxu51xm4W/LnziEtoDvXpNZ9++2f/8VQEnaAshhCFiYmJ49dVX8fT05P333y/SOQ+P\ndlTn5qJu145fnn6aRhYW7GzZkudsbcHNjQZVqvz3iMc//4QZM2DVKqhVK28+2sSJUKVKibzW2xm3\n+ezXz1gStYR72fcY6jqUGV4zcK3rWiLXK02KAidPwrZt8PPPuqcpSWxCiAovIyODuXPn8umnn2Jp\nacmQIUOKdJ620Y58+y2j+vZl+dChWJiYaDYXOuLx7l2YPx8WLICcHJgyBaZPB2trQ1+WVk9qQsvJ\ngZ9+yktm27fDtWt5NZ+trFQUMve9cEoF9wS8BCGEAXbu3Kk4OTkpgDJq1Cjlxo0bRT63R1CQwqFD\nBX56TZr0+JNVKkX55htFsbNTFFCUoUMVJTbWgFfy35LvJStTD0xVqn9cXTGaZaS8uPlF5dzNcyV2\nvdJw966ibNmiKH5+imJtnfc2WlgoSv/+irJypaIkJSlKWNhhxdl5qpLXjyva97302IQQFdrSpUup\nWrUqhw4dKvKSVYqisC05mZ/v3tW6X9tQEU0ZrMxMVHfv0vP2bbyuXoUuXWDrVtBhpKUubmfcZtFv\ni1hydAl3s+9W+B7arVuwc2dez2z/fsjMzOvc9usHAwdCz57w8IIKffvmTVsPCZnB3r1Fu4YkNiFE\nhbZ69WqsrKyKPCftYGoq7125QlR6OlVzc7Ue8+hoR61lsExN4b338Jozp0QGhqTcT8m75fgEJLQr\nV/IS2bZteQsYqNV566W+9FJeMnvmGTD7j0W5+/b1om9fL4yMZhd+0EMksQkhKgSVSqV1pKOtrW2h\n54HqxbwAACAASURBVDw8MCQzJ4dMd3dOu7jgWKUKK5o2xXb0aF5/XH1HRWHfzJkFy2CpVMw4cQKv\nYk5qjya0Ia5DmOE1g5Z1S35NtuLy8OCPbdvg7Nm87W5ueY8fBw6ENm1KbKCoJDYhRPmWnZ3N8uXL\n+eSTTzh8+DCNGjUq0nnaBoYYr1jBuBo1CHnhhbyBIfXrY2JsrH20Y1YWfPcdLFqE6YNv5kcUZxms\nlPspLPp1EYuPLq6QCa2wwR/PPJO3APiAAXmzH0qDJDYhRLmkVqv5/vvvmTZtGleuXKFbt25k6pBI\nZv/wQ4EyWOqxY0nYuhWLF1/UbCsw2vHWLfjoI/jyS7h5E9zcULm6wvnzBa6hbxms8P3hLPluCVlK\nFsaKMXVb1WVXzq4Kl9Du3cubrrdtW95zs9RUsLDIe042axY89xzUqVP67ZLEJoQod06fPk1gYCAn\nTpzAzc2NXbt20bt378euaK0oChF//83ca9f47d49rccUmhovXIDPP4e1a/NGNDz7LLz+Ovj40HPX\nLqY98oxtqrMzvfUogxW+P7zA8jDsgGd8nmHpK0vLfUK7dQvCwvKS2b59jx/8URYksQkhyp1atWpx\n79491q5dy4gRIzA2zl/W9uFnZ1UUhdcGDEDl5sYn164RlZ5OPTMzGpubc0lL7Hx9LEWBAwfyJlPv\n3p3X3fD3h+BgaNZMc9iDcleGlsHKyMlg2vJp+ZMawP+g2tVq5TapXbmSd3sxb8K07oM/SpuU1BJC\nlEuKomjtoWl7dma2YgU57dvTqGNH3n7qKfzr1ePHw4cLHKcpg9WlC2zYkJfQzp4FOzt49VWYMAH+\nYzCKPjJyMth9aTffX/iesJgwMvZlQPeCx3WL60ZEaESxXltfigKnTv07+OPBSjutWuUlsoEDoW3b\nkhv8UZiift9Lj00IUSbu37/PypUr6dSpE+3atSuwv7DbjtqWkMkZOxa3TZs4HhSE6T+9u74+Ppw7\ndoxvpk0j18wMk5wcxv/vf/T9+WcYPhySkvKG6YWGwrBhxVr66kEy23xhM2ExYdzLuUedqnXwc/Pj\n9zO/8zu/FzjHwrh4l63RlUr17+CPbdvKdvCHoSSxCSFK1Z07d/jqq6/4/PPPSUpK4t1339Wa2LT5\nPS2N04UsIWNtbq5JapA39yxt+XJiH557duQIkfxza/H116F792LrdhSWzEa5jWKo61C8GnphamxK\neJWCz9icTzgT9FrxLFuji4cHf4SFQUpK+Rj8YShJbEKIUnH79m3mzp3L8uXLSUtLo1evXkydOhVP\nT898xz36/GxC//6ktWzJV4mJRKWnY1zEJWT2LVpUcO4ZMOOZZ/AKCyuW13Q/5z67L+/m+/PfF0hm\nQ1oMoZtTN0yN83/NPlihOmRDCJnqTCyMLQh6LajUVq6uCIM/DFWmiS0hIQE/Pz+SkpKoU6cOAQEB\nBAQEABAZGUlwcDAqlYrx48cTVEyL8AkhyoapqSkrVqygd+/evPPOO7i7uxc4RtvzswMrVqBu355m\nnTsT4uJCXX9/phY2qVqlyqvTtH49pocOaW2HyUOFjfWhLZnZVrX9z2T2qL49+pZaIoOKN/jDUGWa\n2MzMzFi0aBFt2rQhOTmZli1b0qlTJxo3bkxgYCAHDhzA3t4eDw8PfH19ad78yV7OXPx/e3ceF3W5\n/338xQCyb0LKekJUFDEFXJA8eiMupBy3PHnUygq3un+B2v0r2+wolWaekyll52iiFeaeaeKWiUuI\nuFSmFiGgNqAIDKAM62z3H18YGGfAJZWBrufjMY+Tw3e+czV05u11XZ/ruoS2zMXFBblcjqOjY5PX\nfGhi/kw7bRr9tmwhIz5emnfz9cXh5kXVYWHEbN8uzZUVFYGbG2pvb8jLM3qPu1l7Vh9mW37Zwje/\nfaMPsycfeZKJwRNvK8wepOaKP954o+WKPx6UFv1NeHp64unpCUjb4vTr14/8/HxKS0vp0qUL/v7+\nAEyaNIkdO3aIYBMEM3fq1ClWrlzJ2LFjGTt2rNHPTYWaTqfjVHk5nxUUkHrjhsn72ltZGRSTxERF\nEePjA+vXS7uDrFghTQ6NHg1PPgmPPcaIAwf+0Nqz5sLsieAniPSPNKswa6r4Y+BA+Pe/peKPzp1b\nupUPhtn8VrKzszl//jwDBgxgz549+Pn56X/m6+tLRkZGC7ZOEISmlJSUsHnzZtauXcuJEydwcHAg\nJCTE6Lqb586mxMRQEBTEZwUF/FpZia1MxkOWlhSYeA99H+vKFdi4UQqz06elb+6oKGkDwscfB2dn\n/WvuZu1ZlaqKvdl79aX5ylol7nbuZhtmzRV//POfUvFHhw4t3coHzyx+Q0qlkkmTJrFs2TIcHBxu\nubvAzRYsWKD/58jIyNs+ukIQhD/mxIkTDBo0iNraWoKDg1m+fDnPPPMMLi4uBteZmjvb/+mn0K8f\nAwcNYlVgIBM7dOB7tZqZq1ZxZeZM/XXe//kPcQEBMGwYpKZKE0R9+kg16JMmgZdXk+0bHBNzy0XU\n9WG25ZctfJP1jT7MJveczMTgiWYXZn+G4o96hw4d4tChQ3f8uhb/balUKiZMmMBTTz2lH7rw8fFB\nLpfrr5HL5fj6+jZ5j8bBJgjCgxMSEsLs2bOZPHkyISEhJv9SqlSreX3zZqO5M6ZPZ9C2bRyZPVv/\nlFNVFQMzMrhx+jTVdnbYVlXhLJfjVFUlLaJ64w1pqLFbt9tqX+M9GW0sbIifEk/M8Jhmw+yJHk8w\npNMQswqzP1vxR72bOyoLFy68rde16G9Op9Mxbdo0goODmTNnjv75vn37cuHCBS5duoS3tzebNm1i\nw4YNLdhSQfjzysvLY8OGDcyYMQNXV1eDn7Vr147333+flIMHeXX2bP0w4/QxY1A/8ghbi4rYU1JC\nVRMl+gZbZV25wv5589icm2t03fzwcAanp99RtYOpPRnPLDtDt7Ru/GDzg1GYRfpHYm1pHunwZy/+\n+KNaNNjS0tJITk6mV69ehIaGArB48WIee+wxkpKSGD9+vL7cXxSOCMKDk5uby7Zt29i2bZt+fjsg\nIIAJEyYYXdvcMKNXv35M8/Iiw96ekybex7a0FF59FfbuhTNnmvxCsrS1veNv8RVfrjDak/Fa+DUU\nhxQ899JzLRpmKSlHWLFiPzU1VtjYqImPH0F09GBR/HGPiL0iBUEwkJCQwD//+U8AwsLCmDBhAk88\n8QRdu3Y1ujarspIxc+bwW+ODOev037KF9I8+QmZhYTL8Or/9NsvT04lRqaRv75EjeXPHDt5JTze6\n1/zoaN7eu7fZdqs0Kn6+9jPpeekczzvO1pVbqRlUY3TdoNxBHPnsyC0/h/slJeUIs2fvIyfnXf1z\njo5vANEolYP1xR9jx/55iz+aIvaKFAThrowYMQIHBwcef/xxfrl4kRVff83+jz7CRqfjhTFjcAwL\nI6WkhF0KBdlVVVBba/I+dlZWyFQqSEsjZu9e+P57Eo8ckebOdDriAgKIWb8ehg7VVzOO6Nnztkv0\nCysKSZenk54nPU5dOUWlqhIAL0cvXKxdKKTQ6HX2lvb34mO6Y6WlkJkJ8+btNwg1AKXyXby95/P5\n54PbVPFHSxHBJgh/IuXl5Xz33Xfs3r2biooK1q9fb3TNgAEDGDBggMle1refforu/HlsQkIY4ubG\nbB8fVtfW8rOJ96o6fhzc3UGpBGtrYgYNImbkSBg5Enr0MDm02FSJ/qMjo/nh6g8GQZZbKs3FWcms\nCPUMZXrodCL8IojwjeAvLn9hd6/dD3xPRo1GGkLMzDR+FOoz1vTXbteulowff9+a9qcihiIFoY2r\nrq7mo48+Ys+ePRw9ehSVSoWTkxMhffpg27MntTIZNjod8ePGERMVhUKl4lBZGXNfew35U08Z3S90\n82aOJibiULc11T/69+e0gwM5dcOXAJ0XLqTvjz+ycfJkKciGDAEnp9tuc1FFkRRgdUF28spJfW/M\n09GTCF8pwCL8Iujj1Qc7azuT90n5NsVwT8bJ92ZPxooKyMoyDq+sLKn8vp6Hh3SsW+PHe++9yfff\nv2N0z+jo+ezd+/Yfbltbdrvf9yLYBKGN02q1eHl50aFDB0aNGsXIkSO5XlvL/9uyxaA35rJuHe4R\nEVzs1g0dIFu3Dm3d3q2N/Z/kZA716QPHjsGxYyzIzaWfnR2Jfn76Ev04uZyT/fqx4PDhW7ZPrVVz\n9tpZfU8sXZ5OTqnUy6rvjdWH2ADfATzs8vAdr3W9GzodFBSY7n39/nvDdTKZtBLh5gDr1s300W6m\n5tg6d36d5csfIyZm8H3/92rNxBybIPxJlJWVceTIEQ4dOkRcXJx+Xqy+9D5+3DiysrL0i6ZLVCqG\nxccbrSu7/uyzyD77jIXR0Qx1c2PBV1/xrYn3sz16FNaskQ7nHDgQtZUVMVlZxGRlGVx33M50L6qo\noojjecf1QXYy/yQVqgqgoTc2q8+sW/bG7pXaWsjJMR1gjXf4cnSUAmvwYMMA69Llzo5yqw+vxMT5\nVFdbYmurIS5OhNq9JIJNEFqhY8eOsX37dlJTU/nxxx/RarXY2Njg2L49X16+bBBamV98weOlpVT0\n7Mn316/za2UlNHGmWa927Zh/9CikpTH78GFyz541HGJctoy4iRMhNhb8/cHCghEpKcyZMZ0PrzZs\nhjXby5MJcXGotWrOFZ4zmBvLLskGpN5YiGcIsaGx+h7Z/eyN1Rdv3PzIyZHmxur5+kqBNXWqYYB5\ne9+7dWMxMYNFkN1HYihSEFqR+v0Wf0lLI/+nn+gRHMyE8eOJiooiPDycmJdf5qCJtWYkJeE6cyaP\nOjsz0MWFL157jcxG21bV6z9rFhlZWVKVYkQEKX5+JFZVUe3qiq2lJXFjxxITFWXYpm9TmDt/Oi7F\nBTioocIK5E72dBjRmVyXXKPe2ADfAUT4RtDHuw/21ve2QvH2ijegXTsIDDQePgwMvKOpQOEBE0OR\ngtAKVVZWcvLkSdLT00lPTyc8PJzXX38duGkh9NChYGND2ZYtFAQHs6ZDB57/+Wcyy8tN3revjQ0Z\nhYXIDhyAM2c4k5qK6upVo4KPTgBnz0JQEFhaEgPcXGqh1WnJv5HPhZILZCmyWJy4mN9H3rx1cSVV\naZeJnXt/emN3WrwxZoxhgPn7wx88lk0wYyLYBMEM/CsxkYSEBJQlJei0WgACAwOJiIgA4LpazcKt\nWxuGGOvmy/KffppVSUl4+vrSz8kJlZUVOSbu737kCLJPP5UqHbp2JcjamqkZGSTOmmVY8NG/P/Ts\nCYCiUkGWIqvhUSL97wXFBarUDUOZFuWmwyrUO5QVI1fc9WdyN8UbI0bcunhDaPtEsAnCfZZy8CAf\nbN5MiUIBVVW889JLBsN5KQcPsuLwYcr9/SE6GoKDeejsWcL+9jeOP/IInY4f51J1NVRWmrx/j6Ii\nzv3rX1icOUPK1avMzskx7Il9+CFx48bB3/8OwcFgb486OpqY/fuNCj5W3bhAxJoIshRZlFSV6J+3\nklkR4BZAoHsgwzoNI9A9UP+IvRTLfvYbtctWdnsHej7o4g2h7RPBJgj3wY0bN1i6dCkHUlM5dfYs\n6vJyqQvi5ka8pydanY6ejz7KL5WVvLRxI/IXXzR4fVFEBBuTkugeEEB/jYYZSiUrcnO5ZuK9qnNz\npSG+Pn2I6d0bZDISN2yg2tYWW+CFF/+HwFB/UhRZZJ35nixFFmeDC3jytCXrFQ1VExPd4OewGgKs\n7ZnYY6JBePm7+je5p2L8lHh+fv8cBX+9on/O83tv4l4xXAhtTsUbQtsmgk0Q7oJarSbpyy/5z9at\nOHfqZLDAGaRd75cuXYqlvT3q0FDo1Am6doWuXcn18GDc2rVo6yd5mtiS6q9ZWRytG4oEOGdnx4mF\nC43mxfo7OaE7eZL88nyyFFnIFVn0UBTphxD/nr4SdZpa/xo3Wze6BXdD+fxgpnybh6vWCmtHN6bN\nnsvmcX+/8w+j1gkuDIRfboB1NahsqdE4s32TEzu3Nl+80asXTJwoijeEe0tURQqCCTef9hw/bhzH\nU1M5d+4cmZmZZGdno1bXhcWOHeDsjOfnnzMqOhrrkBAuVFWRdeMGeWvWwHPPGd3fY8UKFgcEEPTT\nT8wvKCD1o4+Mrol++232DhumD8SZzz7N2N+yjRZCL+mg4/RMmX5nDgA7Kzu6uneVelztAw16X+72\n7nf0WVRWQnFx04+vvnqTa9eMd9KA+Xh4vG1UeSiKN4S7JaoiBcGEmwPrxTFjCAsKIjc3l5ycHHJy\ncugRFsb/btlicIrzmdWrsUhLQ9auHS4BAdi2a4dy5EjptMe6hcgFU6eSlJREe3d3upaXE1lYyK78\nfMpMtMP511+Znp8PgYHMftiPi598zKUX/kf/c9dVy6kdascI31SuFK7nau5VamtLuG4LexvNi010\ngxx3J2aGTTMILx9nH2QWMqP3ramB/Pzmg+rmRxNHqSGT1W8Fafpr5NFHLUlLu41fiiDcYyLYhDbD\nVC9r1JAhaLVaLC0tjTf1ffNNvl25El2jCR4LmQz34cMpfvVVg3tfmzFD+iafNo0rgPWaNfDYY0Zt\nePT8edKGDwcbG7R+vkwpKeWUieFDD00RPaZacKX8G67XXAe5HfzrZ7CyA3UVN7yvkG3nhleNF13a\nd2Hww4PZ3XU3u0Mu0y8D/XqxzHDoXhvGOLtlFP8O2T/A8WZCSqls+vNzc5OqCD08wM9POsiy/s+m\nHq6u0kcSHa1mv3HtCE5OGuMnBeEBEMEmmD1TgWW0SPjgQaYuXEhJQABcuwaFhXy3bh2y6mo+/uor\n/CIimLdpk+E2UkFB6Ly94fJlmDABfHzQdexIcXKyyXY8dOkSR6dO5eGCAsZ072Zyu6nfPS2JWBTA\neV0h5aoc+n4By02U1X8Y0A6fh4IY2mkoXk5edHTwwgkvbFXeyCq9UF13p0QhkwLpNymUKi5Foxx/\nmVOBhu95apUtkVsNn3Nyagighx6SlqU1F1Lt24PVXX4bxMePICfnDaO9D+PijINfEB4EEWzCfXE7\nYQSwZMkS/rtzJ1pra2QqFbPGjGHevHkAXL58mS+3buXDlBQKg4JAoYCiIs5dvEhCbS3BERFcU6ko\nqK0lYc0aSiwsYN066Vu6Y0c0YWFoioqYqVBIi45rbjp0si7k/BcvJuHAAfwKC/FVFDPc0YlLJjb/\ntSgp5PGJRWTbqai9mgOJCyGuoSdm9fG7OPfW4telD73beeFs4U1ZzTV2JW00GD6c7tQZnJZTmBjD\nL3U9qZISqFu+ZsTWVgonZ4d4KnblUPW3hpVqHqmdiZ0WR3RkQ0i5uz/Y8nex96FgbkTxiHDHmgsj\nwPRpyevXkzBqFEH+/hQVFVFUVMTGjRv59tw5al5+WTqfC7BfupRRQ4YQPnEin8XHc27HjoY3traW\nvuFfeEHaYiI2tuFn69ZJvS4bG6nkro7rO++wWqXCteoGb+gsOLH0A6N/H/eEWXh1z6LAERR20PkT\nO3TB4UbDh5bnTxL89kxsVd5YVnkh/7WI3/LOodLZoq2wxLl0LDVlURQXg7qhCBFHUuhOIg5UU4kt\n19zicPGNabYH1fhh32jXqft1DIsgtAbi2Jo/idvtGcXOnM7OjAzpi7+mhjHh4SSt+vSO77VkyRLe\n++47yl5/Ha5eBbkc+y++ICoggO69e3OtuJiTeXlkvvSSURvsZ82i8qYFwVhaQny8tOdRvaQkiI3F\nKjMT3ZYtaKZOlXphjo76hUxd3kngVcsarKoU6GoUvFnhQf5Hnxi9p/vsWVSOu0w7q/ZYyJ1QKgJQ\nxzWEcLsPV+Bj6YeLUxjaSjc0Sjd+/2EcvbVVOLg1VB8qS+Wc0TijLGnYOkomk5p1uwHl4SFtwSjW\nYgnC3RFVkS3kVr0ZuDfDdPX3ublnlFN3InLj+8XOnM6W7FyUH66QStwsLNjy/vswczqr/7sapUbD\n9gMHeGXFCoo6dpQqDJRKDj37LO62tkQ88QQBkZGUqVRsSEujom7vQnbsgE2bqAR2nTvHrp07pTPt\n67Zkulltp074Dx2Km0xGB+DohQtUvvqq0Tf9Q79lMuvdGJSySjapunH14YeN7pWjVfO8ayBqCzeQ\nudIlO5nOJoo0LC54oEj4Df0GUHYH4dx2LBzAUgVOZS9ibR8F9uBoL/WOLqui+ckxle6aLBzKQGEF\nmbaetLcZwt6dhsUTomRdEMxPmwg21+AeDPDvzN6Ubwyev1cBcrv3WrJkCYsPHOD6uw2T6Ivr/rn+\nfrcbRvqeUaN7vfvWW/yQlcVfBw+mrLSUNVu3crlXL/j5Z2mlK5Dz5JPE/ecT9hRf5fS+bzm/91vK\nS0qkb+Bhw6SJnGefRfnWW6xdvZq19QdBrl8PHTpIQ3p2duDkRK2jI1evX+crwF6txrmqiprG42Kj\nR8PAgeDigsO6tTzt74yVtpIdl28gN/F7cqooxCX/GqUWzlzCBau83012X2qUFaysPoi9hRuWmvPY\nLEmmZl7DSc6OHyQzyG0BQQ9HYV8XRsWBfcj8Jp4ujYo0bCrUjFm0gkej0V9nbx+FnV0UMuNKeL3o\n6AD2H3qSUx6J0oLjKlsojmPgkOMMHNj06wRBMA9mG2xHjhxhzpw5qNVqZsyYQVxcXJPXXv94Jd8n\nLOSxmNH6cPsjAfLeokUAvPLKK6jVar7et4+XN23i8rhxoFKBnR2/rPuMF37LpWdQDyqqqqmormbR\npk3ciI6GnTul3SRqa7nu7U3Chg2cs7ahVqcj9XAqRfXDdPv2wY4d5NTUMHrDBmRWVuhqa3EcNRKl\npRXa+p5RnfKuXdmclMTmpKSGJ9PSYMoUfbABXLS1Y42LB1btH0Lr54eFjQ260FApsK5fh759AXAo\nLWXMzrVYqKo5oFBQOOcVePppg26I24IEeh06S+V3F6m2dKH6embDuiwfH+kB2BZXIf/LJzjZOtCj\n4w+Ur0imLL4hjB76JJnY4f9iUJ+GMHr6//ZE885CKt5s6GU5vL0QPxcLzh0Lr3smkJSDLiRu3041\nYAvEzZli4i8oMRxJgW8TE7Gsrkbj6cnwuDgGx9z5/JNU5bePnJy9wCEgUlT53aZDhw4RGRnZ0s1o\nNcTndX+YZbBpNBpiY2M5cOAAPj4+9OvXj2HDhhEUFGT6BcnJVPzlYQ4ePkTngQNx9X+Yy9evo/jf\nlw0uy/H1Zdzkydi5uqDVaNBqNFSVlknB0EjZ66/z2n//y/xp01CtXWvwPgD84x/kPf88byQlQVAX\nsJOBnb006fLee0bNq+zZk+SwEGQaDfz0Q8MPrK2leSMPDyzVKtp7e2FtZYm7vQ0XysoxOgpy+HDs\nz5wh3NcXK1k7ThdcoeT1V/ULhOt55CoIX5OJ1ro72h7hpJfP50b9XwzWrZM2wgWsLl0i6vEUHBws\n+KUonsJGRRf1Qns8wq5/L8fWVupcLVmyhPcWLZLm2Oq4vvsuLz85nnnzutc940fKQSfDMJplHEbL\n3lvCjFdmYjl3FhY2duhqqnDQKVmydJXBdTFRUSZ72jcbHBNzV0F2s8ZVfpmZR+nefZCo8rtN4ov6\nzojP6/4wy2A7ceIEXbp0wd/fH4BJkyaxY8eOpoNtzRoAVBYW5F4twNrXF62rm/F1tbVo1Go0Wi0y\nKyusbWyoqa5B6+xsdKldeTm97WwoiRzM79eKqIqKkqrtrK2lMzKAjmXXidl3ECuZBVYWVnyh01C+\nfr10jY2N9LC2pv3c2SQWgp1je97BFn20RUVJD2Dw5u18vWQ51tbSy7sOHsjFmxvl6YmnjQ0Hv/4a\naKL6MDmZ5QvmGgSBf9f30CQkUPHWW/rnHBYuxFXtwPTp0lCgc8dxzF6/3uheL02ZYpCb9UOqq954\nA421NZYqFTNNDN3eThjFDI9h9furGlX5eZpNlV/9CccLFixgwYIFLd0cQRDugFkGW35+Pn5+fvo/\n+/r6kpGR0fQL9u8HS0tc4l6k7PwvAETHxxsfpDF8OCOUSvYuX65/KmDgQC6a2EHCs7SUY7t2Ndzr\n8ceNrgnp5M+axQn6Px/+bhOXklYbDa15OeiYMjESgHa2k5sMkMabv84aM8Zkz2hmo+rB+uAw6BlN\nMe4ZffzhBzw7fS6WL8RRU6bAJuM07cra8fGny+74XiCF281BdrdihseYRZAJgtB2mGW5/7Zt29i7\ndy+rV68GIDk5mYyMDBITE42utfDxgStXjJ4XBEEQ2pbOnTuTnZ19y+vMssfm4+ODXN5QVyeXy/H1\n9TV5rS4//0E1SxAEQWgFmil6bjl9+/blwoULXLp0idraWjZt2sSYxgt4BUEQBKEJZtljs7KyIikp\nifHjx+vL/ZssHBEEQRCERsxyjk0QBEEQ7pZZDkUKgiAIwt1qtcF25MgRwsLC6NWrl8lqScFQbGws\nHTt25JFHHmnpppg9uVzOkCFDCA4OJjIyknXr1rV0k8xadXU14eHhhISEMGDAAJYtW3brFwloNBpC\nQ0MZPXp0SzfF7Pn7+9OrVy9CQ0Pp37//La9vlUORGo2Gbt26GexMsmHDBjEP14yjR4/i6OjI1KlT\nOXv2bEs3x6wVFBRQUFBASEgIxcXF9OzZk9TUVPHfVzMqKyuxt7enpqaGPn368PXXX9OlS5eWbpZZ\n++CDDzh9+jTl5eXs3LmzpZtj1jp16sTp06dp3779bV3fKntsjXcmsba21u9MIjRt0KBBuLmZ2I1F\nMOLp6UlISAgAHh4e9OvXjytirWSz7Os2x1YqlajVamwe5EmnrVBeXh67d+9m+vTpf+pjt+7EnXxO\nrTLYTO1Mki/Wswn3QXZ2NufPn2fAgAEt3RSzptVq6d27Nx07duTFF180+P+nYGzu3LksXboUWXPH\nTAh6FhYWREVFERoaqt+4ozmt8lO1ECc1Cg+AUqlk0qRJLFu2DAcHh5ZujlmTyWScOXOG7OxsVq5c\nyY8//tjSTTJbu3btokOHDoSGhore2m1KS0vjzJkzfPnllyxatIijR482e32rDLY72ZlEEO6G3V8H\ntAAAAjBJREFUSqViwoQJPPXUU4wdO7alm9Nq+Pv7M2rUKA7Xn/MnGDl27Bg7d+6kU6dOTJ48mYMH\nDzJ16tSWbpZZ8/LyAiAoKIjx48dz4sSJZq9vlcEmdiYR7iedTse0adMIDg5mzpw5Ld0cs1dcXExZ\nmXRKn0KhYM+ePaL6thmLFi1CLpdz8eJFNm7cSFRUFJ9//nlLN8tsVVZWUl5eDkBRURG7d+++5X9f\nZrnzyK2InUnu3OTJkzl8+DAKhQI/Pz8SEhJ47rnnWrpZZiktLY3k5GR9eTHA4sWLeczEKRACXL16\nlWeeeQaNRoOnpycvvfQSQ4cObelmtRpiaqV5165dY/z48QC4u7szd+5cRowY0exrWmW5vyAIgiA0\npVUORQqCIAhCU0SwCYIgCG2KCDZBEAShTRHBJgiCILQpItgEQRCENkUEmyAIgtCmiGATBEEQ2hQR\nbIIgCEKbIoJNEFqZX3/9lUWLFrV0MwTBbIlgE4RWJjU1Vb/VlyAIxkSwCUIrsm/fPtasWUNeXh4F\nBQUt3RxBMEtir0hBaGVGjx7NN99809LNEASzJXpsgtCKFBQU4Onp2dLNEASzJoJNEFqRkydP0r9/\nf06ePEllZWVLN0cQzJIINkFoRTp37oxcLqe6uhp7e/uWbo4gmCUxxyYIgiC0KaLHJgiCILQpItgE\nQRCENkUEmyAIgtCmiGATBEEQ2hQRbIIgCEKbIoJNEARBaFNEsAmCIAhtigg2QRAEoU35/6lBcKR9\n/bANAAAAAElFTkSuQmCC\n", "text": [ "<matplotlib.figure.Figure at 0x1082b7390>" ] } ], "prompt_number": 12 }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks good: the solution is converging on the known result as we increase the number of steps, which tells us that the method is consistent, as any useful method must be. We can check that the order of the approximation is indeed $\\mathcal{O}(h^1)$ by plotting a function of the global error at $t=5$, given by $| \\, y_N - e^5 \\, |$, over a large range of stepsizes $h$." ] }, { "cell_type": "code", "collapsed": false, "input": [ "max_N = 16\n", "N = 2**np.arange(2, max_N) #\u00a0N = 2, 4, 8, ..., 2^max_N\n", "\n", "order_check = 1 # for visual check of the order of accuracy\n", "\n", "y_end = np.zeros(len(N)) # array to fill with the final values\n", "stepsize = np.zeros(len(N)) # array to fill with the stepsizes\n", "\n", "for i, N_i in enumerate(N): # loop over different numbers of steps\n", "\n", " t = np.linspace(0., t_max, N_i+1)\n", " y_end[i] = ode_int_ee(exp, y_0, t, solve_args)[-1]\n", " \n", " stepsize[i] = t_max/N_i\n", " \n", "plt.loglog(stepsize, abs(y_end - np.exp(solve_args['a']*t_max)), \n", " 'b-o', label='Global error')\n", "plt.loglog(stepsize, stepsize**order_check,'k--', label=r'$h^1$')\n", "plt.xlabel(r'$h$')\n", "plt.legend(loc=2)" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "pyout", "prompt_number": 13, "text": [ "<matplotlib.legend.Legend at 0x1121d6610>" ] }, { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEnCAYAAAAafRyJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XlclOX6P/APLoBbJqUIzAgI5QaIimuxmUtqYmoKWKiM\nipaKyOmcTDRA07JcDmKnXCs9RmpZcOQVuQKWC/3Uku9BShES0LSUo7jAMHD//niSRFAZmJlnls/7\n9eKl88zMw+XNcnk/z3Vft5UQQoCIiMgENJE7ACIiovpi0iIiIpPBpEVERCaDSYuIiEwGkxYREZkM\nJi0iIjIZTFpERGQymLSIiMhkNNPnyXNzc5GQkAC1Wo1Ro0Zh3Lhx+vx0RERk5qwM0RFDrVZjypQp\nSEpK0venIiIiM6b15UGVSgV7e3t4enrWOJ6ZmYnevXvDy8sLiYmJ1cdTUlIQGBiIiRMnNj5aIiKy\naFrPtA4fPozWrVtj8uTJyM7OBgBUVlaiS5cu2L9/P5ycnNC3b18kJSWhW7du1e8LCgpCSkqKbqMn\nIiKLovU9LV9fXxQUFNQ4lpWVBXd3d7i4uAAAQkJCkJycjCtXrmD37t0QQmDChAm6iJeIiCyYTgox\niouLoVQqqx8rFAocP34c/v7+8Pf3f+h73d3dkZeXp4swiIjIyLm5ueHcuXMNfr9OSt6trKwa/N68\nvDwIIfTyERsbq7f3Peo1D3q+ruP3H3vY44b+myx1vPQ5ZhwvjhfHS/vxaewkpWlcXFyctm/63//+\nh6SkJLz22msAgBs3biAlJQWvvPIKACA1NRVt27bFs88++8hzxcfHV//97uVFXWroOevzvke95kHP\n13X8/mMPepyeno6AgIBHxtZQ5jZegH7HjOOlHY6XdsxpvNLT0/HJJ58gIyMDDUg7fxENkJ+fLzw8\nPKofV1RUiM6dO4v8/HxRXl4uevbsKXJycup1rgaGYLFiY2PlDsHkcMy0w/HSDsdLO439na/15cHQ\n0FAMGjQIv/zyC5RKJT7++GM0a9YMW7ZswdixY9GnTx+oVKoalYOPEhcXh/T0dG1DsUj6nGWZK46Z\ndjhe2uF41U96enrjZlh/Msji4ocGYGUFmUMgIiIDaezvfKPoPciZFhGRebOImZadnR1KSkoMHBHJ\npV27drh27ZrcYRCRHjV2pqXXhrn1FRcXh4CAgFrXhktKSnjp0II0ZukEERm39PR0nVxRM+qZFu93\nWRZ+vYnMn1nc0yIiIqoPo0haLMQgIjJvFlGIYY6Xi6ZOnYrmzZtj48aN9Xp9kyZN8N1332HQoEEN\n+nwBAQEYOnQoYmJiGvR+QzLHrzcR1cTLg0YkLy8PkyZNQpcuXfDYY4/BxcUFL730Eg4dOlT9Gisr\nK4MWHBj68xER6ZNRJC1zuDx47NgxeHp6olOnTjhw4ABu3LiBY8eOYdy4cbV2bDbl2URFRUW9jjXk\nPERkvnR1edBokpa2rVBSUzMxfPgiBATEYfjwRUhNzWzQ59bVeaZPn45XXnkF7777LhQKBQCgY8eO\nmDRpEjZs2PDA9/36668YM2YM2rdvj06dOmH+/PkoKyur8Zrc3Fw888wzsLOzw+DBg2t0Sf7888/h\n7e2Nxx57DI6Ojpg1axZu375d77ivXr2KadOmQalUon379ggODsaVK1eqn3dxccHq1avxwgsv4Mkn\nn8Tu3bsREBCAhQsX4uWXX4aDgwPWrFmDyspKLFmyBG5ubrCzs8OQIUPw3//+t/o8U6dORXh4OObO\nnQtnZ2dERUXVO0YiMn0BAQE6SVqyd6t9WAgPem7Pngzh5rZQAKL6w81todizJ0Orz62r81y+fFlY\nWVmJgwcPPvK1U6dOFdOnTxdCSI2Ge/ToIaZPny5u374tiouLRa9evcTs2bOrX29lZSW6desm8vLy\nxJ07d8SsWbNE9+7dRWVlpRBCiG+++Ubk5OSIyspKkZ6eLpydncWbb75Z/f6AgACxbNmyOmOpqqoS\nzz77rJgwYYI4c+aMuHDhghg/frx47rnnql/j7OwsOnbsKL788kuh0WjEnTt3hL+/v2jZsqX49NNP\nRVlZmbh9+7ZYvny5cHV1FT///LMoLy8XixcvFg4ODuLGjRtCCCGmTJkimjdvLlauXClu3Lghbt++\nXSseI/h2JCI9a+zPuey/JRqStIYNi6mRaO5+tGu3SPj7i3p/tGtX93mGD1+k1b8hKytLWFlZiZ9/\n/rn6WHJysnj88cdF27Ztha2tbfXxe5PW999/L6ysrEReXl718+vXrxctWrSofmxlZVWji3ReXp6w\nsrISR48erTOWxYsXi379+lU/fljS+uGHH0Tz5s1FcXFx9bHs7GxhZWVVfczFxUVMmzatxvv8/f1r\nJDYhhHjqqafE4sWLqx+XlpYKa2tr8fnnnwshpKTl7u5eZxx3MWkRmb/G/pwbdUeMBykvrzvsqqqm\nWn3eqqq6z1NWpt15nJ2dAQBFRUV4+umnAQBBQUEoKSnB999/D19f3xqvv1sYUVhYiHbt2qFz587V\nz/Xv3x9lZWX4448/8OSTTwIA+vbtW/18586d0a5dOxQVFQEA9u3bhyVLluDnn39GeXk5KisrYW9v\nX6+48/PzUVVVhR49etQ43qJFC1y4cAGOjo4AUGf89x8rKipCv379qh+3bt0aPXr0QGFhYfWx+uyv\nRkTmSVcdMYwmaWnDxkZT5/EBAyqRllb/8wwfrsHevbWP29pWahVPhw4d0L17dyQlJWHw4ME1nhMP\nKbpQKpUoKSlBXl4e3NzcAEgFHba2ttUJCwCysrIwatQoAFKFYklJCRQKBdRqNV588UW8+eab2L9/\nP2xsbBAbG4utW7fWK25nZ2c0bdoUubm5D010zZs3f+QxpVKJ48eP44UXXgAAlJaWIicnB0qlEoCU\n6Jo1M4pvNyKSwd2Jyb0b/zaEURRiaCsychjc3GquO3JzW4i5c4fKch4A2LRpE7Zv344FCxagqKgI\nQgjcvn0bx48fr1FyLv7cchoA+vXrh+7du2PFihW4c+cOLl68iA0bNkClUtU49xdffIHz58+jrKwM\nq1atQrdu3dC/f38IIdCmTRvY2dlBCIF9+/bh448/rhXbgxJn37590b9/f0RHR+PHH3+EEAK///47\nPv/880f+e+8/59SpU/HZZ5/h7NmzUKvVWLlyJezs7KqT7cOSNxFRfZnkf31HjfIDACQmLkZZWVPY\n2lZi7tznq48b+jwAMGDAAPz000+IjY1FYGAgrly5gieeeAK9e/fGwYMHq19377qpZs2aYc+ePYiM\njESnTp1ga2uL8ePH4913361x7tdffx2TJ09GTk4OevbsieTkZFhZWcHGxgbr1q3DypUr8dZbb2Hg\nwIGIiorCunXrarz/Qeu0rKyskJycjMWLF2Ps2LH4448/YG9vj2HDhiEkJOSh/977z/n3v/8d5eXl\nGDZsGK5fv45evXph7969aN26da1/NxFRQ7EjBhkNfr2JzB87YhARkcUwiqRlDh0xiIjowdgwl8wO\nv95E5o+XB4mIyGIwaRERkckwyZJ3IiJjlpqaibVr96K8vBlsbDSIjBzWoKU0VJtek1ZycjJSU1Oh\n0Wgwa9asGm1+iIjMUWpqJubN+xZ5ecuqj+XlSU0MmLgazyCFGFeuXEFsbCw+/PDD2gGwEIP+xK83\nmbo7d4AhQxbhyJG3az03fPhipKUtlSEq49LYn3OtZ1oqlQqpqano0KEDsrOzq49nZmYiKioKGo0G\nM2bMwNy5c6ufW7FiBWbOnNngIImIdK0hl/CEAC5fBs6flz7y8v76+/nzwMWLwIN+rWrbiJvqpnXS\nuruR3+TJk6uPVVZWQqVSYf/+/XByckLfvn0xZMgQdO3aFW+88QZGjhwJb29vnQZORNRQD7uE99xz\nfsjPr5mM7iao/Hzg3j1WrawAhQLo3BkYPlz6c+dODe75/3w1bRtxU920Tlq+vr4oKCiocSwrKwvu\n7u5wcXEBAISEhCA5ORn79+/HwYMHUVpainPnzlnsbOvGjRvYt28fNm3ahG+++UbucIgs3tq1e2sk\nLADIy1uG8eMXo7y85myrVSspGbm7/5WY7n44OwO2tjXP3avXMMybF1Pj/FIj7uf19u+xJDopxCgu\nLq7eggIAFAoFjh8/jsTExBqXCR/k3lXS2uyrZSoee+wxjB8/vlYjWyIyjPJy4NQp4Ngx6SM9ve5f\nfY8/3hSzZ0sJyc1N+rN9e2lGVV+6bMRtDnS1j9ZdOklaje3erYvWHsbg/PnzmDNnDuzt7evcIoSI\n9E8IoKgIOHpUSlBHjwInTwJqtfR8p05Au3YaXL5c+73e3pVYvLjxMYwa5WexSep+909EGruflk6S\nlpOTU40dagsLC6FQKOr9fm13LjZWnTt3Rrdu3Xj/jkhH6lMscecOcOLEXwnq2LG7BRHSpTsfH2De\nPGDAAOnD0RFITeUlPEPT1YyrQSXvBQUFGD16dHX1oEajQZcuXXDgwAE4OjqiX79+SEpKQrdu3R4d\ngJmVvPfp0wcpKSlwcnKq9VxgYCAOHTokQ1SmwRS/3qQ/dRVLuLnFYMGC4WjZ0q86Sf34I6D5czPz\nzp2lxDRwoPRnz55AHRtvV58/MXHfPZfwhnJ2ZACN/jkXWgoJCREODg7C2tpaKBQKsWXLFiGEEOnp\n6cLb21t4eHiIhISEep8PgIiNjRWHDh2q8zlTcu3aNeHs7Cz+9a9/ifXr14upU6cKIYS4deuW2LBh\ng3B1dRUpKSkyR2m8TO3rTfo1bFiMkC723f+xSABCtGwpRECAEG++KURyshC//SZ3xPQwhw4dErGx\nsY3+OTf5Lu8PuqT4oGmotq/XRnJyMpYtW4bvvvsO1tbW1bsWP/74440+tyXgTIsAqaT8wAEgIiIO\nv/0WV+v5p56Kw86dcfDwAJqxEZ3JMfjiYn0wl3ta6enp+Nvf/gZra2uo1Wr8/vvvTFhE9fDbb8Ce\nPUBKCrB/v3SfqmlTTZ2v7dy5ErxtbHpkvaelS+Z0T6tPnz7Ys2cPHBwckJKSgq+//hqJiYlo0aIF\nmjRhQ/1HMbWvNzWcEEB2NvCf/0iJKitLOu7sDIweDQQFAbduZeL11++/p7UQCQmWWz5uDjjTMhIl\nJSW4efMmHBwcAEiXCsePH4/t27djxowZMkdHJD+1GsjI+CtR/fqrdLxfP2DpUilReXreuybKD82b\nc72TueBMy8gcPnwYX331FVavXg0A2LhxI4qLi+Hv74/AwECZozMNpvT1pr88rCz96lXgm2+kJJWW\nBpSWAi1aAEOGSElq1Cjgz//nkYVo7M85kxYZDX69TU9dZemdOsXgueeGIy/PD999B1RVAR07Ai+8\nICWq554DWraUMWiSFS8PEpFs6urhd+HCMnz88WJ4eflh4ULpHpWPD8DbupaNlwfJ7PDrbVp+/RUI\nCIhDQUFcref694/DsWO1jxM19uec//chonq7fh3YvBkICABcXICCgrrL0h9/nNtwkH4YRdKKi4vT\naRdgItKdigogNRUICZHuTU2fLvX2W7IE2LRpGNzcYmq8XurhN1SmaMlYpaen66Q5Oi8PktHg19t4\nCCE1od22DUhKAn7/HbCzA0JDgbAwqUz9bmk6e/iRNlg9SGaDX2/5XbgAbN8uJaszZwBra6mQIiwM\nGDFCekzUGGZRPUhE8rlxA/jiCylR3b1K/+yzwPr1wIQJQLt2soZHVINRJC2WvBPp1/0LgGfPHoZm\nzfywbRvw9ddAWZm0nXx8PPDKK9IWH0S6ZBEl73Z2digpKTFwRCSXdu3a4dq1a3KHYXbqWgDcpEkM\nqqqGw87OD8HBwOTJQP/+2m0rT9QQZn1Pi4ga77nnFuHgwbdrHe/VazGOHVvK+1RkUFynRUR1+n//\nD4iIANLT674L8NhjTZmwyOQYxT0tItKN0lKpRH39euDkSak5bceOGly8WPu1trZcAEymhzMtIjNw\n6hQwaxbg6AjMnCltA5KYKC0C3rCBC4DJfBjFTIvVg0Tau3UL+PxzaVb1ww+ArS0QHCwlrQED/iqq\nuLvQl/tSkZwsonqQiGo7fVpKVP/+t7TGqnt3KVGFhXFNFRk/Li4msgC3bwM7dgAbNgDHjgE2NtLC\n35kzgWeeYak6WQ7OtIiMyP2LgIOChiE3V1oEfP060LXrX7OqJ56QO1oi7XGdFpGZqGsRMBCDpk2H\nIzjYDzNnAr6+nFWRaTPqdVr5+fmYPn06JkyYoM9PQ2QW3n239i7AwDL4+e3D9u2Anx8TFpFek5ar\nqys2bdqkz09BZNKqqoBvv5U6qX/3Xd23mKuqmho4KiLjpXXSUqlUsLe3h6enZ43jmZmZ6N27N7y8\nvJCYmKizAInM0fXrwNq1QLduwPPPA1lZQOfOde8CzEXARH/ROmmFh4cjLS2txrHKykqoVCrs3r0b\nJ06cwObNm3HmzBmdBUlkLnJygNdeA5ycgHnzpI0Vt2+X9rFau5aLgIkeReuSd19fXxQUFNQ4lpWV\nBXd3d7i4uAAAQkJCkJycDHt7eyxcuBA//vgjVqxYgTfeeEMXMROZFI0G2LNH6lBx8KBUrh4aCsye\nDfj4/PU6LgImejSdrNMqLi6GUqmsfqxQKHD8+HHY2dnho48+euT74+Liqv/OzhhkLv74A9i0Cfjw\nQ2kmpVQC77wDTJsGtG9f93tGjfJjkiKzoqtOGHfpJGlZNbKk6d6kRWTqTpwA1q2TGteWlwODBwP/\n/KdUbNGMy/nJwtw/EYmPj2/U+XTyI+Tk5ITCwsLqx4WFhVAoFPV+P3sPkqlTq6Ut69etA44eBVq1\nAlQq6RJgjx5yR0ckP1l7DxYUFGD06NHIzs4GAGg0GnTp0gUHDhyAo6Mj+vXrh6SkJHTr1u3RAXBx\nMZmYe7tWABo4OAzDoUN+uHwZeOopKVFNnQq0bSt3pETGx+C9B0NDQ5GRkYGrV69CqVRiyZIlCA8P\nx5YtWzB27FhoNBrMmDGjXgnrLs60yFQ8qGuFjw/w6ad+GDoUaMINf4hq0dVMS+uklZSUVOdxf39/\nnDp1qtEBERmr8nJgwYK6u1Y88cRiDB/OAgoifWPvQaJH+O034KOPpI/Ll+MAxNV6jb9/HNLTax8n\nopq4NQmRnpw8CSQkSBstqtXAqFHAb79pcOJE7deyawWRYRjF1fe4uDid1vETNZRGA3z5pdRNvU8f\n6e8REcDPP0sLhOPj2bWCqCHS09N1sryJlweJAJSUSAuB162TFgK7uACRkVLZ+v1VgKmpmUhM3HdP\n14qhXBBMVE9msZ9WbGwsqwdJFrm5UuPaTz+VdgcOCJB6Ao4eDTRlc3UinblbPRgfH2/6SYszLTKk\nqipg717pflVamtQLcNIkKVn17Cl3dETmzSxmWkxapGv3b1sfGTkM/v5+2LpVmln9/DPg4CB1XJ85\n88G9AIlIt1g9SHSfuhYAHz8eg4oK4PZtP/TtC/z738CECYC1tYyBEpHWjCJpsSMG6dLatbUXAF+/\nvgz29ouxf78fBgzgtvVEhiZr70Fd4uVB0iW1GujZMw65uXG1nuMCYCL5NfZ3vlGs0yJqrD/+AJYv\nB1xdgdxcbltPZK6YtMiknTkjFVIolUBMDODhAcTFcQEwkbniPS0yOUJIJev//KdUsm5rC4SFSSXr\n0t5VfvDx4bb1RMaE97TI4ty5I1X9/fOfQE4O0LGjtHcVS9aJTAdL3snsXboE/OtfUpf1P/4AevUC\ntm4FJk6UFgYTkeVg0iKjdeoUsGaN1GVdowGCgoD58wE/P5asE1kqJi0yKpWVwH/+IyWrzEygdWvg\n1VeBuXMBd3e5oyMiuTFpkWzubbXUtKkG7u7DcOCAH/LygE6dgJUrgWnTgMcflztSIjIWRpG0WD1o\neepqtXTwYAy6dgV27fLDiy8CzYziu5OIdIHVg2TS+vdfhKyst2sdHz58MdLSlsoQEREZAjtikMnQ\naICdO4GBA4GsrLqnUWVl3MSKiB6MSYv07vp1YNUqwM0NCA4Gfv8d6NqVrZaISHtMWqQ3+flAVBSg\nUACvvy71Bfz6a2kvq5Ur2WqJiLSn13ta5eXlePPNN3Hnzh2MGTMGzz//fO0AeE/LrAgBHDkCrF4t\nJagmTYCQEGl9Ve/eNV+bmpqJxMR997RaGspWS0Rmzqh3Lj548CAuX76M0NBQREREYMOGDbUDYNIy\nCxUVwJdfSsnqhx+Adu2k9kpz5gBOTnJHR2Q4QgikpaXhm2++QUJCAqy4Er4GgxdiqFQq2Nvbw9PT\ns8bxzMxM9O7dG15eXkhMTAQAZGdnw83NDQBw586dBgdJxut//wPefx/o3BkIDZUef/ABUFgIvPMO\nExZZjrKyMmzatAkeHh4YOXIkvvzyS/z2229yh2V2tE5a4eHhSEtLq3GssrISKpUKu3fvxokTJ7B5\n82acOXMGXl5eOH/+PACgZcuWuomYjEJeHhAZKd2v+sc/gKeekjpZ5OYCr70GtGold4REhvOvf/0L\nnTp1wowZM2BtbY1t27YhPz8fDg4OcodmdrRevunr64uCgoIax7KysuDu7g4XFxcAQEhICJKTkzF/\n/nzExMTg+++/x7hx43QRLxnYvV0rbGw0eO65YTh61A/JydLi39BQ6X6Vt7fckRLJR61Wo3///oiO\njkZAQAAvCeqRTnoOFBcXQ6lUVj9WKBQ4fvw4bGxssHLlyke+Py4urvrv7IxhPOrqWrF3bwxatwYW\nLvTDa68Bjo4yBkhkJObNm4eoqCi5wzBKuuqEcZdOklZj/1dxb9Ii47Fq1d4aCUuyDAMGLMbbb7PK\njyxHeXk5Pv/8c3z55Zf46quv0LRpzUXwnFk92P0Tkfj4+EadTyfrtJycnFBYWFj9uLCwEAqFot7v\nj4uL02kmpsY5d07qqp6eXvf/aSoq2LWCLMPVq1exfPlyuLq6YurUqcjPz0dxcbHcYZmk9PR0nUxQ\ndJK0fHx8cPbsWRQUFECtVmPHjh0ICgrSxanJQIQADh8Gxo4Fnn4aWL8ecHBg1wqyXKtWrYJSqURM\nTAx69uyJvXv34vTp0+jUqZPcoVk2oaWQkBDh4OAgrK2thUKhEFu2bBFCCJGeni68vb2Fh4eHSEhI\nqPf5GhAC6ZBaLcT27UL06SMEIMQTTwixaJEQFy8KsWdPhnBzWyiklCZ9uLm9KfbsyZA7bCK927Fj\nh5g2bZr4v//7P7lDMSuN/Z1vFF3eY2NjWYBhYCUlwMaNQGIiUFQEdOkiVQGGhQH3rk5g1woyd0II\n3pMygLsFGfHx8cbbEaNeAbAjhkHl5QEJCcCWLcCtW8DgwUB0NDBihNRyichSlJSUYOPGjdi+fTuO\nHDmCVlxcaBDcmoQeSQjgu++AceOkRcAffQSMHw+cOgUcOACMGsWERZbj/PnziIyMhFKpxBtvvIEn\nnngCV65ckTssqiej2BuWOxfrR0UF8MUXwJo1Uj9AOzvgzTeB2bO5voos0zvvvIOYmBg0a9YMoaGh\nmD9/Pry5Mt4guHMxVbu/a4VKNQwXLvhh7VrpftXTT0v3qyZPrnm/isjSZGRk4Ntvv8Xs2bPhxMaY\nsjDqLu/1CoCFGI1SV9cKK6sYCDEcgYF+iI4GRo7k5T+yLGq1GtbW1nKHQfdgIQYBAIYPX4S9e9+u\ndXzgwMU4cmSpDBERyefXX3/F2rVrsW3bNmRnZ8Pe3l7ukOg+jf2dbxT3tEh7Go20f9WRI3V/Ca2t\n2bWCLEdWVhZWrVqFL7/8EgAwceJElJWVyRwV6QMvGpmY69eBVasANzdpR+CqKnatIMu2bNky9O/f\nH99++y2io6ORn5+Pzz77DM7OznKHRnpgFDMtVg8+Wn4+sHYtsHkzUFoK+PtLC4OtrIZh/vyYGve0\n3NwWYu7c52WMlshwXnzxRbRu3RoqlQpt2rSROxx6AFYPWohjx6SZ1e7dUjFFcLBUCdinz1+vYdcK\nsgTXrl2DnZ2d3GFQI5lF9SCTVk0aDfD118Dq1cDRo8DjjwMzZwJz5kg7BRNZkpMnT2L16tXYtWsX\nfvrpJ3Tt2lXukKgRWIhhRm7ckNorJSQABQVA587SJcCpU4HWreWOjshwqqqqkJqaitWrVyM9PR1t\n2rTB7Nmz8dhjj8kdGsnMKJKWpd/TunBBul+1caOUuJ59VpplBQUBTVkESBbo/fffx4IFC6BUKrFy\n5UpMnz4dbdu2lTssagTe0zIx93etiIwchvbt/bBmDbBrl/SaCROk+1X9+skbK5HcLl68iMzMTIwf\nPx7NmzeXOxzSId7TMgF1da2wtY1BWdlwPPaYHyIipJ2CubccWZqff/4ZTz/9NLcGsSDs8m4C1q7d\nWyNhAUBZ2TJ07boPRUXA++8zYZHlEEIgLS0NQ4cORdeuXXH48GG5QyITwqSlZ0VFQE5O3bcO7e2b\ngstKyFKUlZVh06ZN8PDwwIgRI5CTk4N33nkHnp6ecodGJoRJS09OngReeQVwdQWKiti1gmjTpk2Y\nMWMGrK2tsXXrVuTn52PBggVo166d3KGRCWH1oA5VVQF79kiVfxkZQJs20r0qD49hWL6cXSvIsk2e\nPBndu3dHYGAg72FZIFYPGpFbt4CtW6XNFs+eBZRKYN48YPp04G6VLrtWkCUQQiA9PR2+vr5o1swo\n/k9MRobVgzK6eBH44ANp+/pr14C+fYG//U3ayp4/r2RJ1Go1kpKSsHr1apw+fRpffPEFxo8fL3dY\nZITYEUMGP/0kXQJMSpJaLr34IhAdDTzzDMCrHmRJrl69ivXr12PdunW4dOkSPDw8sGXLFowaNUru\n0MhMMWnVU1UVkJYmJasDB4BWrYBZs6TLgG5uckdHJI+9e/ciJiYGw4cPxyeffIKhQ4fyfhXplV4v\nD+bn52PZsmW4fv06dt1t+3B/AEZ0ebCurhWDB/th2zbpflVuLuDoCERGAhERAIueyNJVVFTgl19+\nQY8ePeQOhUyESdzTmjBhgtEnrbq6VrRrFwONZjhKS/3Qq5d0v2rCBMDaWsZAiQysoqICu3btQlBQ\nEFqzczM1kkE6YqhUKtjb29daBJiZmYnevXvDy8sLiYmJDQ7CGNTVtaKkZBlsbfchPR04cQJ4+WUm\nLLIcJSUtCyNKAAAVQElEQVQleO+99+Dq6oqXX375gf/xJDKkeiWt8PBwpKWl1ThWWVkJlUqF3bt3\n48SJE9i8eTPOnDmDbdu2Yf78+bh48aJeAtYHIYBLl+q+vde9e1P4+7PAgizHhQsXMG/ePCiVSrzx\nxhvo2rUrUlNTMWXKFLlDI6pf0vL19a21aj0rKwvu7u5wcXFB8+bNERISguTkZISFhWHNmjVwdHTE\ntWvXMGvWLPz4449YsWKFXv4BjVFWJu1f5ekJZGezawURAJw7dw4ffvghxo8fj1OnTmH//v0YOXIk\nmjRhAx2SX4OrB4uLi6FUKqsfKxQKHD9+vMZr7Ozs8NFHHz3yXHFxcdV/N0RnjN9/l9ZWrVsHXLkC\neHkBUVHDkJISg/Pn2bWCLFtgYCAuXLiAjh07yh0KmQFddcK4q8FJS5dlrfcmLX3KzZWqALdulWZZ\nI0ZIxRWDBwNWVn4YMgRITFx8T9eK59m1gszSjRs3sHnzZoSEhMDBwaHGc1ZWVkxYpDP3T0Ti4+Mb\ndb4GJy0nJycUFhZWPy4sLIRCoWjQufTZe1AI4NAhaX1VaipgYwNMngxERQHdu9d87ahRfkxSZNZ+\n/fVXrF27Fhs3bkRpaSlatWqFiIgIucMiC6CrGVeDk5aPjw/Onj2LgoICODo6YseOHUhKSmp0QLqi\nVgM7dkjJ6scfgfbtgbg44NVXgQ4d5I6OyLB++eUXvPXWW/jiiy8AAMHBwZg/fz58fHxkjoxIS6Ie\nQkJChIODg7C2thYKhUJs2bJFCCFEenq68Pb2Fh4eHiIhIaE+p6qlniHU29WrQixfLoSDgxCAEN27\nC7FpkxB37uj00xCZlOzsbNG2bVvx+uuviwsXLsgdDlmwxv7ON4qGubGxsY2+PHj2LPDPfwKffALc\nvg0MHSr1Axw+nOXqRABw584dtGjRQu4wyELdvTwYHx9v/B0xHhqAlquj72+1NHjwMBw96oeUFKB5\nc2kB8Pz5Uhk7kSUpKirCunXrqvetIjJGFtXlva5WS3v3xqB1ayAmxg+zZwMseiJLc/LkSaxevRo7\nduxAVVUVXF1dmbTIbBnFasG4uLh6VZXU1WoJWIYBA/Zh6VImLLIsOTk5CAwMRJ8+fZCcnIw5c+bg\n3LlzmDlzptyhEdWSnp6uk+VNJnV5MCAgDhkZcbWO+/vHIT299nEic3bx4kX4+flh1qxZmDFjBtre\n3SabyIiZxeXB+q7TsrFhqyWiuxwdHXH27FnuX0UmQVfrtExqplXXPS03t4VISGDnCjJPp0+fxpo1\nazB16lT4+/vLHQ5Ro5nFTKu+7iYmtloicyaEwLfffotVq1Zh//79aNmyJQYNGsSkRQQTm2kRmbv/\n/ve/mDhxInJycuDo6Ii5c+ciIiICdnZ2codGpBNmMdPSZ+9BIlPSqVMntG/fHlu3bkVwcDCsueso\nmQmLvKdFRESmrbG/841inRaRpRBC4MCBAxg5ciR27twpdzhEJodJi8gAysvL8emnn8Lb2xtDhgzB\nyZMncefOHbnDIjI5vKdFpGe5ubkYPHgwLl26BA8PD2zZsgWTJk2CjY2N3KERGQzvaRGZCI1GA5VK\nhbCwMAwZMoSLgcmiNfZ3PpMWkY4IIVBVVYWmTZvKHQqR0WIhBpHMKioqsH37dvj4+OCDDz6QOxwi\ns8akRdRAJSUleO+99+Dq6opXXnkFt27dgoODg9xhEZk1oyjEIDI158+fh5eXF27duoXBgwdj/fr1\nGDFiBJo04f8DifTJKJIWqwfJ1Li6uiIqKgovvfQSvL295Q6HyOixepDIADQaDcrLy9GqVSu5QyEy\nCyzEINKD69evY/Xq1XBzc8Py5cvlDoeI/mQUlweJjMWvv/6KhIQEbNq0CaWlpfDz84Ovr6/cYRHR\nn5i0iP508eJFuLu7QwiB4OBgzJ8/Hz4+PnKHRUT30Ps9reTkZKSmpkKj0WDWrFno169fzQB4T4uM\nyMaNG/H8889DqVTKHQqRWTKZjhhXrlxBbGwsPvzww5oBMGmRgd28eRO3b99Ghw4d5A6FyOIYrBBD\npVLB3t4enp6eNY5nZmaid+/e8PLyQmJi4gPfv2LFCsycObPBgRI1VlFRERYsWAClUomFCxfKHQ4R\nNUC9k1Z4eDjS0tJqHKusrIRKpcLu3btx4sQJbN68GWfOnMG2bdswf/58XLx4EUII/OMf/8DIkSO5\nnoVkcfLkSYSFhcHV1RXvv/8+hg4dihkzZsgdFhE1QL0LMXx9fVFQUFDjWFZWFtzd3eHi4gIACAkJ\nQXJyMhYsWICwsDAAwNq1a3Hw4EGUlpbi3LlznG2RQV29ehUDBw6EjY0N5s6di8jIyOrvVyIyPY2q\nHiwuLq5xw1qhUOD48eM1XhMZGYnIyMiHnicuLq767+yMQbr0xBNP4Ouvv8agQYPQtm1bucMhsji6\n6oRxV6OSlq72Bbo3aRE1xG+//YbS0lI89dRTtZ4bMWKEDBEREVB7IhIfH9+o8zWqI4aTkxMKCwur\nHxcWFkKhUGh9nri4OJ1mYrIcp0+fRnh4OJydnREdHS13OET0AOnp6TqZoDQqafn4+ODs2bMoKCiA\nWq3Gjh07EBQU1OigiB5GCIG0tDQMHToUPXv2xM6dOzFjxgysWbNG7tCISM/qvU4rNDQUGRkZuHr1\nKjp06IAlS5YgPDwcGRkZiIqKgkajwYwZMx55/6pWAFynRVq6efMmlEolWrZsiblz5yIiIgJ2dnZy\nh0VE9WAyi4sfGICVFWJjY1mAQVo5deoUevToAWtra7lDIaJ6uFuQER8fb/pJizMtqktOTg5KS0vR\nv39/uUMhIh3h1iRkVoQQ2L9/P0aOHIkePXrg9ddflzskIjIiRpG0WD1IGo0Gn376Kby9vTF06FCc\nOHECS5YswVdffSV3aESkA7qqHuTlQTIKGo0G7u7uaNOmDaKjoxEaGgpbW1u5wyIiHWvs73yj2E8r\nLi6OhRgWrlmzZjh8+DAUCoXOFq0TkfHQVWcMzrTIYIQQyMjIQGlpKUaPHi13OEQkAxZikNGrqKjA\n9u3b4ePjg8DAQCxdulTukIjIRDFpkd5oNBq89957cHV1xSuvvILbt29jw4YNyMjIkDs0IjJRvKdF\netO0aVPs3LkTXbt2xYYNG/D888+jSRP+P4nIEvGeFpmEmzdvonXr1nKHQURGgve0SFYajQY7d+7E\nhg0b6nyeCYuIdIlJixrkxo0bWLNmDdzd3REcHIzNmzdzxkxEemcUSYsdMUxHVVUVXn/9dSgUCkRH\nR8PZ2Rlff/01jhw5wvVVRPRA7IhBshkzZgxatWqF6Oho+Pj4yB0OEZkQs9iahEnLtFRVVbEKkIga\nhIUYpHM3b95EYmIi3nrrrTqfZ8IiIrkYxTotMg5FRUVITEzEhg0b8L///Q8BAQGcVRGRUeFvI4IQ\nAtOmTYOrqytWrlyJoUOH4siRIzh06BATFhEZFaOYabEjhrysrKzQunVrzJkzB5GRkXB1dZU7JCIy\nM+yIQUREJoeFGFQvly5dwqJFizB58mS5QyEiajAmLTN3+vRphIeHw8XFBcuXL8etW7dQUVEhd1hE\nRA1iFPe0SD+Cg4Oxc+dOtGzZEhEREZg3bx7c3d3lDouIqMH0mrRyc3ORkJAAtVqNUaNGYdy4cfr8\ndHSffv36oVevXoiIiICdnZ3c4RARNZpBCjHUajWmTJmCpKSk2gGwEKPRuJaKiEyFQQoxVCoV7O3t\n4enpWeN4ZmYmevfuDS8vLyQmJtb53pSUFAQGBmLixIkNDpLqdubMGURERMDPz4+Jn4gsQr2SVnh4\nONLS0mocq6yshEqlwu7du3HixAls3rwZZ86cwbZt2zB//nxcvHgRABAUFITvv/8eH3/8se6jt0BC\nCBw4cAAjR45E9+7dsW3bNnh4eKCsrEzu0IiI9K5e97R8fX1RUFBQ41hWVhbc3d3h4uICAAgJCUFy\ncjIWLFiAsLAwAEBGRgZ2794NIQQmTJig08At1dixY5GcnAx7e3ssXboUs2bNwpNPPil3WEREBtHg\nQozi4mIolcrqxwqFAsePH6/xGn9/f/j7+z/yXPfuscLOGA8XEhKCMWPGYNKkSbCxsZE7HCKih9JV\nJ4y7Gpy0dLnhny42BjM3t27dQqtWrWodDwkJkSEaIqKGuX8iEh8f36jzNbjkzMnJCYWFhdWPCwsL\noVAoGnQu7lwsEUIgIyMDY8aMgaenJzQajdwhERHphK52Lm5w0vLx8cHZs2dRUFAAtVqNHTt2ICgo\nqNEBWaKKigps374dPj4+CAgIwJEjRxAWFga1Wi13aERERqVe67RCQ0ORkZGBq1evokOHDliyZAnC\nw8ORkZGBqKgoaDQazJgxA5GRkdoHwHVaGDNmDFJSUtClSxdER0cjLCwMLVq0kDssIiKda+zvfKPo\n8h4bG2vRBRiHDh3C7du3MWLECC4SJiKzdLcgIz4+3vSTliXMtIQQKCoqqlFxSURkabg1iZHTaDTY\nuXMnBgwYAC8vL9y8eVPukIiITJZRJC1zrB68fv06Vq9eDTc3NwQHB+PatWtYtmwZmjZtKndoREQG\np6vqQV4e1JOJEydi165d8PPzQ3R0NF544QUmLCKyeCzEMFKnT5+GWq2Gj4+P3KEQEcmOhRhGoLKy\nEj/++CP69OkjdyhERCaBhRgyuHnzJhITE/H0009j4MCBuHTpktwhERFZBCYtLRQVFWHBggVQKpWI\njIxEx44dkZSUhA4dOsgdGhGRRWhww1xdiouLM4l7WvHx8diyZQvGjx+P6OhoDBgwQO6QiIhMgq66\nvfOelhYuXLiAqqqq6j3EiIhIO7ynpWO3b9/Gf/7znzqf69SpExMWEZGMmLT+dOnSJSxatAhKpRJB\nQUE4e/as3CEREdF9jCJpydkRIzs7G+Hh4XBxccHy5cvh5+eHw4cPw93dXZZ4iIjMETti6MiiRYuw\nZs0ahIeHIyoqismKiEiPzKIjhpwhlJSUQAgBOzs72WIgIrIULMSohytXruCDDz6oc6DatWvHhEVE\nZCLMOmnl5ORgxowZ6NSpE+bMmYPs7Gy5QyIiokYwy6T1/fffY8SIEejRowf+/e9/Y+rUqThz5gy8\nvLzkDo2IiBrBLDtinDhxAqdOncLSpUsxa9YsPPnkkzo5LxERNQw7YjxEWVkZAMDW1lan5yUiosax\n2EKMX375BW+88QYqKipqPWdra8uERURkhkx2pjVu3Dikpqbi8OHD6Nevnx4iIyIiXbPYdVrnz59H\nq1atYG9vr4eoiIhIH4z+8uCtW7fQt29fpKam6vS8nTt3ZsIiIrIwek9a7733HoKDg/X9aSyGXD0a\nTRnHTDscL+1wvAyrXklLpVLB3t4enp6eNY5nZmaid+/e8PLyQmJiYq337du3D927d0f79u11Ey3x\nB6QBOGba4Xhph+NlWPVKWuHh4UhLS6txrLKyEiqVCrt378aJEyewefNmnDlzBtu2bcP8+fNx8eJF\nZGRk4NixY/jss8+wceNGg/cYbOg3U33e96jXPOj5uo7ff+xRj/WF46Udjpd2OF7a4XjVrV5Jy9fX\nF+3atatxLCsrC+7u7nBxcUHz5s0REhKC5ORkhIWFYc2aNXB0dMTbb7+NNWvWYNKkSYiIiICVlZVO\ng38UftG1w/HSDsdLOxwv7XC8HkDUU35+vvDw8Kh+vGvXLjF9+vTqx9u2bRNz5syp7+mqubm5CQD8\n4Ac/+MEPC/hwc3PTOk/cq8FtnHQ1azp37pxOzkNEROavwdWDTk5OKCwsrH5cWFgIhUKhk6CIiIjq\n0uCk5ePjg7Nnz6KgoABqtRo7duxAUFCQLmMjIiKqoV5JKzQ0FIMGDcIvv/wCpVKJjz/+GM2aNcOW\nLVswduxY9OnTByqVCt26ddN3vEREZMFkb+NERERUX0bd5V1fLaDMUW5uLl599VVMmzYNu3fvljsc\no5ecnIyIiAioVCpkZWXJHY7Ry8/Px/Tp0zFhwgS5QzFq5eXliI6OxquvvlprbSvV1pDvK6OeacXG\nxqJNmzbo1q0bRo0aJXc4JkGtVmPKlClISkqSOxSTcOXKFcTGxuLDDz+UOxSTMGHCBOzatUvuMIzW\nwYMHcfnyZYSGhiIiIgIbNmyQOySToM33ld5nWmwBpZ2GjhcApKSkIDAwEBMnTjREqEahMeMFACtW\nrMDMmTP1HabRaOx4WSJtxiw7Oxtubm4AgDt37hg8VmOg9++xRq3yqofMzExx8uTJGguTNRqNcHNz\nE/n5+UKtVouePXuKnJwcsXXrVhEVFSWKi4tFTEyMiIqKEsOGDRNjxowRVVVV+g7VKDR0vO41evRo\nQ4ctm4aOV1VVlfj73/8u9u/fL2P0htfY76+XXnpJjrBlpc2YHTx4UCQlJQkhhIiIiJArZFlpM153\nafN9pfekJUTtbhpHjhwRw4cPr378zjvviHfeeafO937yySciNTVV7zEak4aMV3p6uoiMjBRz584V\nW7duNVisxqAh45WQkCD69OkjZs2aJT766CODxWoMGjJeV69eFTNnzhTu7u7i3XffNVisxqK+Y1ZW\nVib+9re/iTlz5oi0tDQ5QjUK9R2vhnxfNbgjRmMUFxdDqVRWP1YoFDh+/Hidr50yZYqhwjJa9Rkv\nf39/+Pv7Gzo0o1Sf8YqMjERkZKShQzNK9RkvOzs7fPTRR4YOzWg9aMxsbGywcuVKGSMzTg8ar4Z8\nX8lSPWjoxrmmjuOlHY6Xdjhe2uOYaUeX4yVL0mILKO1wvLTD8dIOx0t7HDPt6HK8ZElabAGlHY6X\ndjhe2uF4aY9jph2djpfO78DdJyQkRDg4OAhra2uhUCjEli1bhBBS4YC3t7fw8PAQCQkJ+g7DZHC8\ntMPx0g7HS3scM+3oe7yMenExERHRvYy6jRMREdG9mLSIiMhkMGkREZHJYNIiIiKTwaRFREQmg0mL\niIhMBpMWERGZDCYtIiIyGUxaRERkMpi0iAxs586deOaZZ+QOg8gkMWkRGZinpyf69esndxhEJolJ\ni8jAjh49Ch8fH7nDIDJJbJhLZGARERFwdnaGQqFA8+bNMWnSJLlDIjIZnGkRGVhubi5mzpyJ0aNH\nIysrS+5wiEwKkxaRAd28eRN2dnZ48skncezYMXh7e8sdEpFJYdIiMqAffvgBAwcOBACkpKRg0KBB\nOHnypMxREZkOJi0iA8rNzUVgYCAAQKlU4rvvvoOXl5fMURGZDhZiEBGRyeBMi4iITAaTFhERmQwm\nLSIiMhlMWkREZDKYtIiIyGQwaRERkclg0iIiIpPBpEVERCbj/wOhn/POuGD07wAAAABJRU5ErkJg\ngg==\n", "text": [ "<matplotlib.figure.Figure at 0x1121deb50>" ] } ], "prompt_number": 13 }, { "cell_type": "markdown", "metadata": {}, "source": [ "We see that the slope is constant below $h \\approx 0.1$, which tells us that the order is constant. We could find that order formally by fitting this logarithm, but it's good enough for now to compare the global error with the function $h^1$ (black dashed line). That the slopes are the same indicates visually that the method is $\\mathcal{O}(h^1)$. (Try comparing with $h^2$ by changing `order_check` in the code above to $2$.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Next\n", "\n", "In [Part 2][part2] we'll look at a method with a higher order of accuracy: the **Runge-Kutta method**; we'll go beyond first-order to **higher-order ODEs** and solve a more interesting physical problem; and we'll look at **multistep methods** where more than one point is used to get to the next step.\n", "\n", "[part2]:./7_The-Runge-Kutta-Method-Higher-Order-ODEs-and-Multistep-Methods.ipynb" ] } ], "metadata": {} } ] }