{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Solving Ordinary Differential Equations with the Runge-Kutta Methods " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## List of Problems \n", "\n", "\n", "\n", "- [Problem midpoint](#problem_midpoint)\n", "\n", "- [Problem tableau](#problem_tableau)\n", "\n", "- [Problem Runge Kutta4](#problemrk4)\n", "\n", "- [Problem embedded](#problem_embedded)\n", "\n", "- [Problem coding A](#prob_a)\n", "\n", "- [Problem coding B](#prob_b)\n", "\n", "- [Problem coding C](#prob_c)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Assignment handin -- upload a single, fresh notebook that contains your answers\n", "\n", "* Assignment: A409 do Problems 2 (tableau), 3 (Runge Kutta4), 4 (embedded), Coding A and Coding B\n", "\n", "* Assignment: Grads do Problems 2 (tableau), 3 (Runge Kutta4), 4 (embedded), Coding A, Coding B and Coding C\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Objectives\n", "In this lab, you will explore Runge-Kutta methods for solving ordinary\n", "differential equations. The goal is to gain a better understanding of\n", "some of the more popular Runge-Kutta methods and the corresponding\n", "numerical code.\n", "\n", "Specifically you will be able to:\n", "\n", "- describe the mid-point method\n", "\n", "- construct a Runge-Kutta tableau from equations or equations from a\n", " tableau\n", "\n", "- describe how a Runge-Kutta method estimates truncation error\n", "\n", "- edit a working Octave code to use a different method or solve a\n", " different problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Readings\n", "\n", "\n", "There is no required reading for this lab, beyond the contents of the\n", "lab itself. However, if you would like additional background on any of\n", "the following topics, then refer to the sections indicated below.\n", "\n", "**Runge-Kutta Methods:**\n", "\n", " - Newman, Chapter 8\n", "\n", " - Press, et al.  Section 16.1\n", "\n", " - Burden & Faires  Section 5.4\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "Ordinary differential equations (ODEs) arise in many physical\n", "situations. For example, there is the first-order Newton cooling\n", "equation discussed in , and perhaps the most famous equation of all, the\n", "second-order Newton’s Second Law of Mechanics $F=ma$ .\n", "\n", "In general, higher-order equations, such as Newton’s force equation, can\n", "be rewritten as a system of first-order equations . So the generic\n", "problem in ODEs is a set of N coupled first-order differential equations\n", "of the form, \n", "\n", "$$\n", " \\frac{d{\\bf y}}{dt} = f({\\bf y},t)\n", "$$ \n", " \n", "where ${\\bf y}$ is a vector of\n", "variables.\n", "\n", "For a complete specification of the solution, boundary conditions for\n", "the problem must be given. Typically, the problems are broken up into\n", "two classes:\n", "\n", "- **Initial Value Problem (IVP)**: the initial values of\n", " ${\\bf y}$ are specified.\n", "\n", "- **Boundary Value Problem (BVP)**: ${\\bf y}$ is\n", " specified at the initial and final times.\n", "\n", "For this lab, we are concerned with the IVP’s. BVP’s tend to be much\n", "more difficult to solve and involve techniques which will not be dealt\n", "with in this set of labs.\n", "\n", "Now as was pointed out in , in general, it will not be possible to find\n", "exact, analytic solutions to the ODE. However, it is possible to find an\n", "approximate solution with a finite difference scheme such as the forward\n", "Euler method . This is a simple first-order, one-step scheme which is\n", "easy to implement. However, this method is rarely used in practice as it\n", "is neither very stable nor accurate.\n", "\n", "The higher-order Taylor methods discussed in are one alternative but\n", "involve higher-order derivatives that must be calculated by hand or\n", "worked out numerically in a multi-step scheme. Like the forward Euler\n", "method, stability is a concern.\n", "\n", "The Runge-Kutta methods are higher-order, one-step schemes that makes\n", "use of information at different *stages* between the\n", "beginning and end of a step. They are more stable and accurate than the\n", "forward Euler method and are still relatively simple compared to schemes\n", "such as the multi-step predictor-corrector methods or the Bulirsch-Stoer\n", "method. Though they lack the accuracy and efficiency of these more\n", "sophisticated schemes, they are still powerful methods that almost\n", "always succeed for non-stiff IVPs." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Runge-Kutta methods" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Midpoint Method: A Two-Stage Runge-Kutta Method \n", "\n", "The forward Euler method takes the solution at time $t_n$ and advances\n", "it to time $t_{n+1}$ using the value of the derivative $f(y_n,t_n)$ at\n", "time $t_n$ \n", "\n", "$$y_{n+1} = y_n + h f(y_n,t_n)$$ \n", "\n", "where $h \\equiv \\Delta t$." ] }, { "cell_type": "markdown", "metadata": { "title": "[markdown" }, "source": [ "![fig1](images/euler.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Figure euler: The forward Euler method is essentially a straight-line approximation\n", "to the solution, over the interval of one step, using the derivative at\n", "the starting point as the slope. \n", "\n", "The idea of the Runge-Kutta schemes is to take advantage of derivative\n", "information at the times between $t_n$ and $t_{n+1}$ to increase the\n", "order of accuracy.\n", "\n", "For example, in the midpoint method, the derivative at the initial time\n", "is used to approximate the derivative at the midpoint of the interval,\n", "$f(y_n+\\frac{1}{2}hf(y_n,t_n), t_n+\\frac{1}{2}h)$. The derivative at the\n", "midpoint is then used to advance the solution to the next step. The\n", "method can be written in two *stages* $k_i$,\n", "\n", "
eq:midpoint
\n", "$$\n", "\\begin{aligned}\n", " \\begin{array}{l}\n", " k_1 = h f(y_n,t_n)\\\\\n", " k_2 = h f(y_n+\\frac{1}{2}k_1, t_n+\\frac{1}{2}h)\\\\\n", " y_{n+1} = y_n + k_2\n", " \\end{array}\n", "\\end{aligned}\n", "$$ \n", "\n", "The midpoint method is known\n", "as a 2-stage Runge-Kutta formula.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![fig2](images/midpoint.png)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Figure midpoint: The midpoint method again uses the derivative at the starting point to\n", "approximate the solution at the midpoint. The derivative at the midpoint\n", "is then used as the slope of the straight-line approximation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Second-Order Runge-Kutta Methods\n", "\n", "As was shown in lab 2 , the error in the forward Euler method is\n", "proportional to $h$. In other words, the forward Euler method has an\n", "accuracy which is *first order* in $h$.\n", "\n", "The advantage of the midpoint method is that the extra derivative\n", "information at the midpoint results in the first order error term\n", "cancelling out, making the method *second order* accurate.\n", "This can be shown by a Taylor expansion of equation\n", "[eq:midpoint](#eq:midpoint)\n", "\n", "
Problem midpoint
\n", "\n", "Even though the midpoint method is second-order\n", "accurate, it may still be less accurate than the forward Euler method.\n", "In the demo below, compare the accuracy of the two methods on the\n", "initial value problem \n", "\n", "
eq:linexp
\n", "\\begin{equation}\n", "\\frac{dy}{dt} = -y +t +1, \\;\\;\\;\\; y(0) =1\n", "\\end{equation}\n", "\n", "which has the exact\n", "solution \n", "\\begin{equation}\n", "y(t) = t + e^{-t}\n", "\\end{equation}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Why is it possible that the midpoint method may be less accurate\n", " than the forward Euler method, even though it is a higher order\n", " method?\n", "\n", "2. Based on the numerical solutions of [eq:linexp](#eq:linexp), which method\n", " appears more accurate?\n", "\n", "3. Cut the stepsize in half and check the error at a given time. Repeat\n", " a couple of more times. How does the error drop relative to the\n", " change in stepsize?\n", "\n", "4. How do the numerical solutions compare to $y(t) = t + e^{-t}$ when\n", " you change the initial time? Why?" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "******************************\n", "context imported. Front of path:\n", "/Users/phil/repos/numeric\n", "back of path: /Users/phil/.ipython\n", "******************************\n", "\n", "through /Users/phil/repos/numeric/notebooks/lab4/context.py\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZzN1R/H8dfHkF3WlCwzSmUJ2bcYWaIUKiGyE0VUokiU0EKlUpIQ2SJr1rEM2fetKHv2nQZjmZnz++Nc8xvTLHfGXWc+z8djHs3c7/d+v5+5mfc993zP9xwxxqCUUirlSuPtApRSSrmXBr1SSqVwGvRKKZXCadArpVQKp0GvlFIpnAa9UkqlcBr0yutE5A8RCfZ2HckhIpdFpLC361AqIRr0yuuMMcWNMaHO7Csih0SktptLiu/coSLSIeZjxpgsxpgDbjxnfxExCf3OItJVRDaJyHURGeeuWpT/SuvtApTyFBERQIwxUd6uxRki8gDwAnAikV2PAx8BTwIZ3V2X8j/aoldeF7OVLiIDROQXERkvImGObp1yjm0TgILAXEeXSS/H45VEZI2IXBSR7TG7gRyt8EEishq4ChQWkbYisttx/AMi8kqsehqKyDYR+VdE9otIPREZBDwOfOM49zeOfY2IPOio4aSIBMQ4TmMR2eH4Po2IvOM43jnH75gzkZfmG6A3cCOhnYwxM4wxs4Bzib7YKlXSoFe+6FlgCpAdmIMNPIwxLwP/AM84ukw+FZH7gXnYFm1OoCfwq4jkiXG8l4FOQFbgMHAaaABkA9oCX4hIGQARqQCMB952nL86cMgY0xf4HejqOHfXmAUbY9YBV4AnYjz8EjDJ8f3rQCOgBpAPuACMiO8FEJEmwA1jzHxnXjClEqJBr3zRKmPMfGNMJDABKJXAvi2B+Y79o4wxIcAm4KkY+4wzxvxhjIkwxtw0xswzxuw31gpgMba1DtAeGGOMCXEc75gxZo+TdU8GmgOISFZHDZMd214B+hpjjhpjrgMDgBdE5D/dpyKSBRgM9HDyvEolSINe+aKTMb6/CmSIKxAdCgFNHN02F0XkIlANuC/GPkdiPkFE6ovIOhE579j/KSC3Y3MBYH8y654EPCci6YHngC3GmMMx6pwZo8bdQCSQN47jfABMMMYcTGYdSt1Gg175m9jTrR7BhmL2GF+ZjTEfx/UcRwj/CgwF8hpjsgPzAYlxvAecPPftG435E9s1VJ/bu21uHbd+rDozGGOOxXGoWsDrjj7/k9g3n19EpHdC51cqPhr0yt+cAmKOW/8ZeEZEnhSRABHJICLBIpI/nuffBaQHzgARIlIfqBtj+49AWxGp5biAer+IPBLPueMyCdsfXx2YFuPxkcAgESkEICJ5RKRhPMeoBZQASju+jmO7fuLs0xeRtCKSAQgAbr0GOqJORdOgV/5mCPCeowukpzHmCNAQ6IMN7yPYC6lx/ts2xoRhg/gX7AXRl7AXfG9t34DjAi1wCViB7XYBGI7tV78gIl/FU99kIBhYZow5G+Px4Y7zLBaRMGAdUDGeGs8ZY07e+sJ28VwwxlwGEJE+IrIgxlPeA8KBd7DXLMIdjykF2DHF3q5BKaWUG2mLXimlUjgNeqWUSuE06JVSKoXToFdKqRTOJ4dg5c6d2wQGBnq7DKWU8hubN28+a4zJE9c2nwz6wMBANm3a5O0ylFLKb4jI4fi2adeNUkqlcBr0SimVwmnQK6VUCueTffRxuXnzJkePHuXatWveLsXvZMiQgfz585MuXTpvl6KU8gK/CfqjR4+SNWtWAgMDsSvCKWcYYzh37hxHjx4lKCjI2+UopbzAb7purl27Rq5cuTTkk0hEyJUrl34SUioV85ugBzTkk0lfN6VSN78KeqWUSqnWHV3HsDXD3HJsDXoPGzduHF27dk10n+PHj0f/3KFDB/788093l6aU8pJ1R9dRd0JdRm4eSdj1MJcfX4PeB8UO+tGjR1OsWDEvVqSUcpdbIZ83S16Wt15O1vRZXX4ODfokaNSoEWXLlqV48eKMGjUKgCxZstC3b19KlSpFpUqVOHXqFABz586lYsWKPPbYY9SuXTv68VvCwsIICgri5s2bAPz7778EBgYybdo0Nm3aRIsWLShdujTh4eEEBwdHTwmxcOFCypQpQ6lSpahVq5YHf3ullKvFDvn82eJbAfPO+M3wyph69IBt21x7zNKl4csvE95nzJgx5MyZk/DwcMqXL8/zzz/PlStXqFSpEoMGDaJXr1788MMPvPfee1SrVo1169YhIowePZpPP/2UYcP+3/+WNWtWgoODmTdvHo0aNWLKlCk8//zzNGnShBEjRjB06FDKlSt32/nPnDlDx44dWblyJUFBQZw/f961L4JSymNih/y2lfkZthQ+/xxcPX7CL4PeW7766itmzpwJwJEjR9i7dy933XUXDRo0AKBs2bKEhIQAdtx/06ZNOXHiBDdu3IhzDHuHDh349NNPadSoEWPHjuWHH35I8Pzr1q2jevXq0cfKmTOnK389pZSHxA75rSvy8/zzUKoUXLkCWbK49nx+GfSJtbzdITQ0lCVLlrB27VoyZcpEcHAw165dI126dNHDFwMCAoiIiACgW7duvPnmmzz77LOEhoYyYMCA/xyzatWqHDp0iBUrVhAZGUmJEiUSrMEYo0MllfJzsUN+S2h+XnjB9iosXuz6kAfto3fapUuXyJEjB5kyZWLPnj2sW7cu0f3vv/9+AH766ad492vVqhXNmzenbdu20Y9lzZqVsLD/XnmvXLkyK1as4ODBgwDadaOUn4kd8puX25B/7DEb8tmzu+e8GvROqlevHhEREZQsWZJ+/fpRqVKlBPcfMGAATZo04fHHHyd37tzx7teiRQsuXLhA8+bNox9r06YNnTt3jr4Ye0uePHkYNWoUzz33HKVKlaJp06Z3/osppTxi3dF1PPnzk9yT+R6Wt17OpmU25MuUcW/IA4gxxn1HT6Zy5cqZ2AuP7N69m6JFi3qpIveZPn06s2fPZsKECW49T0p9/ZTyB7dCPk+mPIS2CWXj0vy8+CKUKwcLF8Ldd9/5OURkszGmXFzb/LKPPqXo1q0bCxYsYP78+d4uRSnlJrFDfsOS/DRtakN+0SLIls39NWjQe9HXX3/t7RKUUm4UO+TXh+SnWTMoX9625D0R8qB99Eop5RaxQ37dYtuSr1DBsyEPGvRKKeVysUN+7SLbkq9UyfMhDxr0SinlUuuPrr8t5NcszE/z5lC5MixYAFldP5VNojTolVLKRdYfXU/dn+tGh/zqBfl56SWoUgXmz/dOyIMGvdsMGDCAoUOHersMpZSHxAz55a2Xs2q+DfmqVZ0MeWPg6lW31KZB7yNuTZ2glPI//w35ArRoAdWqwbx5TkxrYAy88w7UqAGXL7u8Pg36JPj555+pUKECpUuX5pVXXiEyMpIsMf4PTp8+nTZt2vznefv376devXqULVuWxx9/nD179gD2Dtg333yTmjVr0rt3b0/9GkopF4od8r/PK0DLllC9um3JOzV3zcCB8Omndtxl5swur9E/x9F7YZ7i3bt3M3XqVFavXk26dOl49dVXmThxolOH7tSpEyNHjqRIkSKsX7+eV199lWXLlgHw999/s2TJEgICAlzyayilPCd2yK+YW4DWrW3DfO5cJzN72DDo3x9at4ZvvnH9HMX4a9B7wdKlS9m8eTPly5cHIDw8nHvuuSfR512+fJk1a9bQpEmT6MeuX78e/X2TJk005JXyQwmF/G+/QaZMThxk5Ejo2ROaNIHRoyGNezpZEg16ERkDNABOG2P+M4+uiLQAbvU7XAa6GGO2O7bVA4YDAcBoY8zHLqnaC/MUG2No3bo1Q4YMue3xmIuJXLt27T/Pi4qKInv27GyL5xNIZjd8TFNKuVfskF8+uwBt2kDNmrYl71TIjx8PXbpAgwbw88+Q1n3tbmfePsYB9RLYfhCoYYwpCQwERgGISAAwAqgPFAOai4jfLnxaq1Ytpk+fzunTpwE7RfDhw4fJmzcvu3fvJioqKnpRkpiyZctGUFAQ06ZNA+wbxvbt2z1au1LKdWKH/LJZNuSfeCIJIT99OrRtC7VqwbRpcNddbq050aA3xqwE4p343BizxhhzwfHjOuDWoocVgH3GmAPGmBvAFKDhHdbrNcWKFeOjjz6ibt26lCxZkjp16nDixAk+/vhjGjRowBNPPMF9990X53MnTpzIjz/+SKlSpShevDizZ8/2cPVKKVfYcGwDdX+uS+5MuVneejlLZxaIzmunQ37ePKLvoJo9GzJkcHvdGGMS/QICgV1O7NcT20UD8MKt7x0/vwx848z5ypYta2L7888///OYcp6+fkrdmfVH15tsQ7KZwsMLm38u/mPGjDFGxJg6dYy5etXJgyxZYkz69MaULWvMxYsurQ/YZOLJVJd1ColITaA9UO3WQ3G9ryTw/E5AJ4CCBQu6qiyllLpjG45toM6EOuTOlJvQ1qGE/FqADh2gTh2YNQsyZnTiIKtXw7PPQpEidn5iV0xC7ySXXOIVkZLAaKChMeac4+GjQIEYu+UHjsd3DGPMKGNMOWNMuTx58riiLKWUumOxQ37x9AK0bw9169qeF6dCfvNmeOopuP9+CAmBXLncXndMdxz0IlIQmAG8bIz5O8amjUAREQkSkbuAZsCcOz2fUkp5SuyQXzTNtuTr1bMteae613ftsu8KOXLA0qVw771urzs2Z4ZXTgaCgdwichToD6QDMMaMBN4HcgHfih3oH+FomUeISFdgEXZ45RhjzB9u+S2UUsrFYof8gqkFeOUVqF8fZsxwMuT37oXate3OS5dCgQKJP8cNEg16Y0zzRLZ3ADrEs20+oOvkKaX8SuyQnz+lAJ07296XGTMgfXonDnL4sB2OExkJy5fDAw+4ve746J2xSikVw4ZjG6g7oW50yM+bXIAuXeDpp+HXX50M+ePHbciHhdmQL1rU7XUnRCc1c6E5c+bw8cdx3/ybxamZjeL21FNPcfHixQT3GTduHMePx3utWynlhFshnytTLpa3Xs5vkwpE37zqdMifOWO7a06dsstJlS7t9roToy16F3r22Wd59tlnXX7c+fMT7/0aN24cJUqUIF++fC4/v1KpQeyQn/tzQbp2hWeesTevOhXyFy7YC68HD9qQr1jR7XU7Q1v0Tjp06BCPPPIIHTp0oESJErRo0YIlS5ZQtWpVihQpwoYNGxg3bhxdu3YF4ODBg1SuXJny5cvTr1+/6OOEhoZSvXp1GjduTLFixejcuTNRUVEATJ48mUcffZQSJUrcNm1xYGAgZ8+e5dChQxQtWpSOHTtSvHhx6tatS3h4ONOnT2fTpk20aNGC0qVLEx4e7tkXRyk/F1/IP/usna3AqZAPC7Od+H/8ATNn2tnNfIRftuh7LOzBtpOunaa49L2l+bJewpOl7du3j2nTpjFq1CjKly/PpEmTWLVqFXPmzGHw4ME0atQoet/u3bvTpUsXWrVqxYgRI247zoYNG/jzzz8pVKgQ9erVY8aMGVSpUoXevXuzefNmcuTIQd26dZk1a9ZtxwTYu3cvkydP5ocffuDFF1/k119/pWXLlnzzzTcMHTqUcuXKue5FUSoViB3ycyYUpFs3aNgQfvnFyWlowsPtu8LGjfadoV5C04N5nrbokyAoKIhHH32UNGnSULx4cWrVqoWI8Oijj3Lo0KHb9l29ejXNm9sBSy+//PJt2ypUqEDhwoUJCAigefPmrFq1io0bNxIcHEyePHlImzYtLVq0YOXKlXHWUNrR51e2bNn/nFcp5bzYIT/rJxvyjRolIeSvX4fnnoMVK+yMlLEaZ77AL1v0ibW83SV9jM9vadKkif45TZo0cS4FKPEsIBD7cRG5NR9QkmoICAjQbhqlkmnjsY23hfzMcQXp0QMaN4apUyFdOicOEhFhJyhbuBB++AFeesntdSeHtujdpGrVqkyZMgXgPytRbdiwgYMHDxIVFcXUqVOpVq0aFStWZMWKFZw9e5bIyEgmT55MjST08WXNmpWwsDCX/g5KpVQbj22kzoQ65MyYk+WtlzNjrA35555LQshHRkKbNrY/fvhw6BDn7UQ+QYPeTYYPH86IESMoX748ly5dum1b5cqVeeeddyhRogRBQUE0btyY++67jyFDhlCzZk1KlSpFmTJlaNjQ+Vmd27RpQ+fOnfVirFKJiBnyoW1C+XVMQd54A55/HqZMcTLkjbGLhkycCIMHw+uvu73uOxLftJbe/ErJ0xQvX77cPP300x4/b0p5/ZS6ExuObjB3D7nbBH0ZZA5fPGyGDjUGjGnSxJgbN5w8SFSUMd272yf26ePWepOCBKYp1ha9UipViN2S/+WHgtHLtU6c6GRLHqBfP9tV0707fPSRW2t2Fb+8GOvPgoODCQ4O9nYZSqUqsUN+6qiC9OoFL75oQ97p5VqHDIFBg2x//BdfQDwDLnyNX7XojZMjU9Tt9HVTqVnskJ/yvQ35Zs2SGPJffQV9+tiRNSNH+k3Igx8FfYYMGTh37pyGVhIZYzh37hwZPLEupVI+JnbIT/quIL172xGREyYkIeR//NF21TRuDD/9BAEBbq3b1fym6yZ//vwcPXqUM2fOeLsUv5MhQwby58+f+I5KpSCbjm+6LeQnflswukH+009JCPlJk6BjR3u36+TJSXii7/CbitOlS0dQUJC3y1BK+YFNxzdRe3zt6HHyP48oSN++0KIFjBuXhKyeNQtatYLq1ZMwfaXv8ZuuG6WUckbskJ/wTSH69oWWLZPYkl+0CJo2hXLlYO5cyJTJrXW7kwa9UirFiB3y478uRL9+8PLLtiXvdNf6ihV2zppixWDBAsia1Z1lu50GvVIqRYgd8j99VYj337c9L2PHJiHk16+3K40EBcHixXZRbz+nQa+U8nuxQ37sl4Xo399ORTNmTBJCfts2e9H1nntgyRLIk8edZXuMBr1Syq/FFfIffABt28Lo0UkI+d27oU4d202zdCmkoNXa/GbUjVJKxRYz5Je1Ws6YLwrx4Yf/D/k0zjZl9++3i3kHBNiWfGCgO8v2OA16pZRfujVOPkfGHNEhP3AgtG8Po0YlIeSPHLEhf/26vQj70ENurdsbNOiVUn7nVshnz5Cd5a1CGT2sUPQUNN9/n4SQP3nShvyFC7BsGZQo4da6vUWDXinlV2KH/A9DCzF4sL15deTIJIT8uXO2T/7YMTu6pmxZt9btTRr0Sim/ETvkR31WiCFDoFMn+O67JIT8pUvw5JOwdy/MmwdVq7q1bm/ToFdK+YXYIT/yk0J88gl07gwjRiQh5K9cgaefhu3b7RQHtWq5tW5foEGvlPJ5sUP+u48L8emndjW/b75JQshfu2bveF271q4b+PTTbq3bV2jQK6V82ubjm28L+W+HFOKzz+DVV23IOz0t/I0bdjmpJUvsfAhNmrizbJ+iN0wppXzW5uObqT2hdnTIjxhsQ/6115IY8pGRdlaz336Db7+F1q3dWrev0Ra9UsonxQz5Za2W8/VHhfj8c+ja1S725HTIR0XZwfXTpsHQoba/J5XRFr1Syuf8J+QHBvL559CtWxJD3hj7zvDTTzBgALz1ljvL9lka9EopnxI75L/6MJAvvrAr+Q0fnsSQ79XLjrt8+214/3231u3LNOiVUj4jdsgP/yCQL7+EHj3giy+SuB73hx/arppXX4VPPvGrxbxdTYNeKeUTYof8lwMCGT4c3ngDPv88iTk9dKjtqmnTBr7+OlWHPDgR9CIyRkROi8iueLY/IiJrReS6iPSMte2QiOwUkW0isslVRSulUpbYIf9F/0C++grefBOGDUtiTn/7re2qefHFJE5hmXI58wqMA+olsP088DowNJ7tNY0xpY0x5ZJYm1IqFdhyYgt1JtTh7vR3s6zVcob1C+Trr6FnT9swT1LI//STHXvZoAFMmJCEyehTtkSD3hizEhvm8W0/bYzZCNx0ZWFKqZRvy4kt1B5fm2zps7G8dShD3wtkxAjbIP/00ySG/C+/QLt2ULu2HUp5111uq9vfuPszjQEWi8hmEemU0I4i0klENonIpjNnzri5LKWUt8UM+WWtQvmsbyDffmsHyiT52ulvv0GLFlC5sp2/JkMGt9Xtj9wd9FWNMWWA+sBrIlI9vh2NMaOMMeWMMeXypJB1GpVScYsr5L/7Dt55Bz7+OIkhv2QJvPAClC5tZ6LMnNltdfsrtwa9Mea447+ngZlABXeeTynl+2KH/CfvBjJyJLz7LgwenMSQX7UKGjaEIkVg4UK4+2631e3P3Bb0IpJZRLLe+h6oC8Q5ckcplTrEDvmP3wlk1Cjo0wcGDUpiyG/aZGefzJ/ftupz5XJb3f4u0bluRGQyEAzkFpGjQH8gHYAxZqSI3AtsArIBUSLSAygG5AZmiv0/lxaYZIxZ6I5fQinl+2KH/JDegYweDX37wsCBSQz5nTvtwiE5ctiQz5vXbXWnBIkGvTGmeSLbTwL549j0L1AqmXUppVKQ2CE/uFcgP/4I/frBBx8kMeT//tsuAZghg13ntUABt9WdUujslUopt4od8oPeDmTMGDv1zIABSQz5Q4fsilBRURAaCoULu6foFEaDXinlNrFD/qOegYwdC/3725BPkmPHbMhfvmxD/pFH3FBxyqRBr5Ryi5ghv/TlUAa+Fci4cTbg+/dP4sFOn7Y3Qp0+bfvkS2mvcFJo0CulXC52yH/4ZiDjx9v++CTPFnzhAtStC4cP2yGUFSu6peaUTINeKeVSsUP+gzcCmTDBzhrcr18SDxYWBvXrw+7dMGcOVI/3nkuVAA16pZTLxA75AT0C+fln+OgjO4wySa5etZOTbdoE06fb4ZQqWTTolVIuETvk+3cPZOJEeyNUnz5JPNj16/Dcc/D77zBxIjRq5JaaUwsNeqXUHdt6Yiu1x9cma/qsLGkZyvuvBzJpkp3S4N13k3iwmzehWTNYtMjOJ988wVt5lBM06JVSd2Tria3UGl/rtpCfPBmGDLGTlCVJZCS0bm1noBw+HNq3d0vNqY0GvVIq2f4T8t2CmDLFTjPcq1cSDxYVBa+8ApMn248Cr7/ulppTIw16pVSyxA75fl2DmDrVLhjy9ttJPJgxdnHYH3+0V22T3N+jEqJBr5RKstgh/95rQfzyC3z2mV0CMMneew+++gp69LAznCmX0qBXSiVJzJAPaRFK31eDmDbNLuL95pvJOODgwfarY0f4/PMkTn6jnKFBr5RyWuyQ79MliF9/tfn8xhvJOOCXX9qumhYt4LvvNOTdRINeKeWUrSe2UntC7f+E/Bdf2B6XJPvhB/vu0LgxjBsHAQGuLlk5aNArpRJ1K+Sz3JWFxS8t593OQcyYYRvk3bsn44ATJ9oRNvXq2VE2aTWK3ElfXaVUgmKG/KLmy3m3c2FmzrTD3JM1AnLmTDtWvkYNmDED0qd3ec3qdhr0Sql4xQ75d14pzOzZdoBMt27JOODChdC0KZQvbycpy5jR5TWr/9KgV0rFKb6Q/+YbeO21ZBwwNNT2xxcvDgsWQNasri5ZxUODXin1HzFDfmHz5fTqWJi5c2HECHj11WQccN06OxNlUBAsXgzZs7u8ZhU/DXql1G22ndx2e8h3KMxvv8G330KXLsk44Nat9qLrvffa1aHy5HF5zSphGvRKqWjbTm6j1vhaZE6XmYXNl/N2+8LMm2eHuHfunIwD/vmnXR0qWzZYuhTy5XN5zSpxGvRKKeD2kF/UPJSe7Qozfz58/z106pSMA+7bZ9d5DQiwIV+okMtrVs7RoFdK/Sfk32pXmAUL7iDk//kHatWCGzfsRdgiRVxdskoCDXqlUrnYIf9m28IsXGhvXO3QIRkHPHnStuQvXoRly6BECZfXrJJGg16pVCyukL+1sFOy1vw4e9aG/LFjEBICZcu6vGaVdBr0SqVSt114bRZKj9aFCQmxId+uXTIOeOmSXcB73z6YNw+qVHF5zSp5NOiVSoXiCvklS+y6H23bJuOAly/DU0/Bzp12ioNatVxes0q+NN4uQCnlWVN3TaXmTzX/E/JjxiQz5MPDoWFDe1PUpEnw9NMur1ndGQ16pVKJs1fP0nR6U5r92oyHcj3EwqYr6N7KhvzYsdCmTTIOeuMGNGliL7qOHQsvvODqspULaNeNUqnAnL/m0GluJ86Hn2fQE4Po+lgvGjdMy/Lldir4Vq2ScdCICLtgyK07qpJ1EOUJGvRKpWCXrl2ix6IejNs2jlJ5S7Go5SLyUoqGz8DKlTB+PLRsmYwDR0XZK7bTp9s1BJN126zyFO26USqFCtkfQonvSjB++3j6Pt6X9R028MeyUhQvbrvTJ0xIZsgbY6evnDABPvggmQvFKk/SFr1SKczlG5fpFdKL7zZ9x8O5HmZNuzUUSluRZk1g1iyoVMl21zz8cDIObgy8/TaMHAm9ekG/fq4uX7mBtuiVSkFW/bOK0iNLM3LTSN6o9AZbOm3lwO8Vo6eA/+wzWLUqmSEPMGCA7ap57TX4+GNdzNtPJBr0IjJGRE6LyK54tj8iImtF5LqI9Iy1rZ6I/CUi+0TkHVcVrZS63bWIa/Rc3JPqY6sTZaIIbRNK79Kf06JpRl56CR56CLZtg54972AN7k8/hQ8/tMNzvvpKQ96PONOiHwfUS2D7eeB1YGjMB0UkABgB1AeKAc1FpFjyylRKxWfjsY2U+b4Mw9YO45Wyr7C98w6OralOsWK2Ff/pp7YV/8gjd3CSESOgd2+7DODo0ZBGOwP8SaJ99MaYlSISmMD208BpEYl9l0QFYJ8x5gCAiEwBGgJ/JrtapVS0G5E3+GjlRwz+fTD3ZrmXhS0WUjrrk7R5ya65XbGi7Yu/o4AHOz6+a1d45hl7ATbZHwmUt7jzbfl+4EiMn486HouTiHQSkU0isunMmTNuLEsp/7fz1E4qjq7IwJUDaVGyBTu77OLCpicpXtwOa//kE1i92gUhP3WqncKyTh345RdIl84l9SvPcueom7g68Ex8OxtjRgGjAMqVKxfvfkqlZhFREQxdM5T3l79Pjow5mNl0JpVzNKJDS9uKr1DBtuKLFnXByebOteMvq1Sx89dkyOCCgypvcGfQHwUKxPg5P3DcjedTKkX7+9zftJ7VmjRTq/gAAB6iSURBVHVH1/F80ef59qnvWD4vD8Vfs3OKffKJHdKe1hV/1UuW2OkMSpe2HxEyZ3bBQZW3uDPoNwJFRCQIOAY0A15y4/mUSpGiTBRfr/+ad5e+S4a0GZj03CSeuKcZr7YRfv3VtuLHjoVirhrqsGqVnaTs4Ydh0SK73qvya4kGvYhMBoKB3CJyFOgPpAMwxowUkXuBTUA2IEpEegDFjDH/ikhXYBEQAIwxxvzhnl9DqZTp4IWDtJvTjtBDoTxV5Cl+eOYHVi3IR/GaEBZmh7K/9ZaLWvEAGzfa6Ybz57cLh+TM6aIDK29yZtRN80S2n8R2y8S1bT4wP3mlKZV6GWMYvWU0by5+E0EY/cxoGtzfjq7thOnToXx52xfvslY8wI4dduGQXLnsYt5587rw4MqbdAoEpXzMsX+P0XFuRxbsW0DNwJqMbTiW9YsLUeJJ+PdfGDLE3vjkslY8wF9/2ZE1mTLZkM8fZ9tN+SkNeqV8hDGGSTsn0XVBV65HXOfr+l/TJPBVunZME92KHzsWihd38YkPHrQrQhljL8IWLuziEyhv09vblPIBp6+c5oVpL9ByZkuK5i7K9s7byXu4KyWKp2HOHNuKX7PGDSF/7JgN+atXbZ/8HQ+8V75IW/RKedmM3TPo/FtnLl2/xCe1P6HVg2/xepcApk2DcuVsX7zLAx7g9GmoXRvOnLHdNaVKueEkyhdo0CvlJRfCL/D6wtf5ecfPPHbvYyxrvIw9K0tQsjFcugSDB9sZgV3aF3/L+fO2T/7wYVi40I7RVCmWBr1SXrBw30Laz2nPqcun6F+jP50e6csb3dLxyy9QtqxtxZco4aaTHzsGjRrBnj327tfq1d10IuUrtI9eKQ8Kux7GK3Nfof7E+mTPkJ31Hdbz6LkBlC6ZjpkzYdAgu/qT20J+/nx7t+uff8K0aVC3rptOpHyJtuiV8pDQQ6G0nd2WwxcP06tKL7o9+gFvdc8Q3YpftsyNAX/jBvTpYxcNKVnSTlamF15TDQ16pdzs6s2r9Fnah+Hrh/Ngzgf5ve3vnNxYlTIl4eJF+Ogjuyqf2yaGPHAAmjWzd7126WLDPmNGN51M+SINeqXcaN3RdbSe1Zq/z/1N1/Jd6fnYx/R+IzNTp3qgFQ92auGOHe1qUNOnw/PPu/FkyldpH71SbnA94jp9lvah6piqhN8MZ8nLS6h57WsqlM7MjBkwcCCsXevGkA8Ph86d7YpQRYvC1q0a8qmYtuiVcrFtJ7fRamYrdp7eSbvS7ehb7nP69rybKVOgTBl7X1LJkm4sYPduG/A7d9rxmYMG6YIhqZwGvVIuEhEVwcerPuaDFR+QO1Nu5jafy80/GlC5DFy4YFvxvXu7MXONsXMkdOtm56yZPx/q13fTyZQ/0aBXygV2n9lN61mt2Xh8I81KNOPDit/Qv1cuJk+Gxx7zQCs+LMx21UyaBDVrws8/Q758bjyh8ica9ErdgcioSIavH06fpX3IclcWfnnhF9L+3YRqZW0r/sMP4Z133NxzsmWL7ao5cAA++AD69tUFvNVtNOiVSqb95/fTdnZbfv/nd559+FmGVPmeQe/ey6RJHmrFGwNff2374XPntkN4atRw4wmVv9KgVyqJjDF8v/l7ei7uSUCaAMY1HEe2g614ooJw7pxtVL/7rptb8efPQ7t2MHs2PP20nTMhd243nlD5Mw16pZLgyKUjtJ/TnpADIdQpXIehj//IJ30LMGmSnVlg0SIPTAK5ejU0bw4nT9qbn954w46TVyoeGvRKOcEYw/jt43l94etERkXy3dPfce/RV6hbyYOt+Kgou0js++9DoUI28MuXd+MJVUqhQa9UIk5dPsUrv73C7L9mU61gNb4MHscX/R5g4kQPtuJPnoSXX7YrQDVtCt9/D3ff7eaTqpRCg16pBEz7Yxpd5nXh8o3LDKs7jKDT3WlQJYCzZ2HAANuKv+suNxcREgItW9oFY0eNgg4dtKtGJYlOgaBUHM5dPUfzX5vz4vQXCcoRxPJmW9n67Zs81yiAvHnt/GD9+7s55CMi7IyTTz4JuXLBhg3/n7dGqSTQFr1Ssfz29290nNuRs1fPMrDmQEpcfIfnHk/L2bM23Pv08UAr/p9/7AXXNWugfXsYPhwyZ3bzSVVKpS16pRz+vf4v7We355nJz5AnUx6WvLiRv354j8YN03LPPbYVP2CAB0J+1izb+b9jh73TdfRoDXl1R7RFrxSw9MBS2s1px9F/j/JutXcpd7k/zWum58wZO8ilb18PBPz16/bmp6+/trOfTZ0KDz7o5pOq1EBb9CpVu3LjCt3md6P2hNpkSJuBBS+s5uhPg3m+UXry5LHd4h984IGQ37sXKle2Id+9u+2y0ZBXLqItepVqrTmyhtazWrPv/D66V+xOtRuDaVM7k2db8QATJ9oJydKls3e6PvusB06qUhNt0atU51rENXqF9KLamGpEREUw57nlXJj8JU0aZiJ3bg+24q9csdMYtGxpB+Jv364hr9xCW/QqVdl8fDOtZrXizzN/0qlMJ2qboXSul5VTp6BfP3jvPQ+14nfutDc+7dljPzoMGABp9c9RuYe26FWqcDPyJgNCB1Dpx0pcvHaRX55dwPVfv+fFhlmjh6h/+KEHQt4Ye1drhQp2YrLFi+3q4Bryyo30X5dK8Xad3kXrWa3ZcmILLUu2pEGar+jRIAenTtkWfL9+HmrFX7pkb3iaNg3q1oXx4yFvXg+cWKV22qJXKVZkVCSfrv6UsqPKcuTSEcY//Stp50ygWaMc5MwJ69fb5f08EvIbNthJ6mfMsBOTLVigIa88Rlv0KkXae24vrWe1Zu3RtTxX9DmeT/8dbze8h1OnbJd4v36QPr0HComKgi++sMtM5csHK1dClSoeOLFS/6ctepWiRJkovtnwDaVGlmL32d18/+TPZF0wnRaN7iFnTli3znaJeyTkz5yBZ56Bnj2hQQPYulVDXnlFoi16ERkDNABOG2NKxLFdgOHAU8BVoI0xZotj2yEgDIgEIowx5VxXulK3O3zxMO3mtGPZwWXUe7AeL2UZzbsv3M/Jkx5uxQOsWAEvvQRnz9qboF57TScjU17jTIt+HFAvge31gSKOr07Ad7G21zTGlNaQV+5ijGHM1jE8+t2jbDi2gS9rjeLepfNp1eh+smf3cCs+MtIOwn/iCTs/zbp10LWrhrzyqkRb9MaYlSISmMAuDYHxxhgDrBOR7CJynzHmhItqVCpex8OO03FuR+bvnU9wYDCt7x7Le80COXnSzjL5/vsebMUfPw4tWkBoqF0kZMQIyJrVQydXKn6u6KO/HzgS4+ejjscADLBYRDaLSKeEDiIinURkk4hsOnPmjAvKUimZMYbJOydT4tsSLDu4jI9rDCdo5VLaNg7k7rttQ3rQIA+G/IIF9u7WDRvsQt3jx2vIK5/hiqCP6zOpcfy3qjGmDLZ75zURqR7fQYwxo4wx5Ywx5fLkyeOCslRKdebKGV6c/iIvzXiJh3M/zNfFtvHNy6/z07g0vPsubNkC5TzVUXjjhr3Y+tRTdlTN5s3QurWHTq6Uc1wxvPIoUCDGz/mB4wDGmFv/PS0iM4EKwEoXnFOlUrP3zKbTb524EH6B/lWHcHhyTzqOSUuxYnaIukfXyj54EJo1s634Ll1g2DDImNGDBSjlHFe06OcArcSqBFwyxpwQkcwikhVARDIDdYFdLjifSoUuXrtI61mtaTS1Efmy5mN48U382PYdxo9Lyzvv2Ia0R0N++nR7A9Rff9k7Xb/9VkNe+SxnhldOBoKB3CJyFOgPpAMwxowE5mOHVu7DDq9s63hqXmCmHX1JWmCSMWahi+tXqcDi/YtpN7sdJy+f5O2K/Tgz/T1e/fEuihaFtWvttDEeEx4Ob74JI0dCxYoweTIEBXmwAKWSzplRN80T2W6A1+J4/ABQKvmlqdTu8o3LvL34bUZuHknR3EXpVWgmn3Usz/Hj9kbT/v0hQwYPFrR7t51xcudOuxLUoEF2DnmlfJxOgaB80srDK2kzqw2HLh6iW9mehM0eSPfRGbzTijfGjqTp2hUyZYL586F+fQ8WoNSd0SkQlE8JvxnOW4veInhcMCLCZ8VXMOvVzxg/JgO9e9sRNR4N+bAwOya+XTt74u3bNeSV39EWvfIZG45toPWs1uw5u4f2JbsQseBTenbPwiOP2CVUK1b0cEFbt9qumv377d2ufftCQICHi1DqzmmLXnndjcgbvLfsPar8WIXLNy4zuOhiFnf/lgk/ZqFXL5u3Hg15Y+z8NJUq2eX+li2zt9hqyCs/pS165VU7Tu2g1cxWbD+1nRbF2hCw5Av6vJmdRx6B1att1nrU+fPQvj3MmgVPP2375nPn9nARSrmWtuiVV0RERTDk9yGUG1WOk5dP8sEjs/m951h+Hp2dt9+2rXiPh/yaNXZs/Lx59uanuXM15FWKoC165XF7zu6h9azWbDi2gUZFmpDt92/p/3ZuHn7YS634qCj45BM7j3GhQrYIj959pZR7adArj4kyUXy1/iveXfoumdJlou9DU5jwTlOOHrXD0j/4wAs3l546ZUfVhITAiy/CqFFw990eLkIp99KgVx5x8MJB2s5uy4rDK6hXuAG5145iUO/7eOghWLUKKlf2QlFLlkDLlnbR7u+/twt367zxKgXSPnrlVsYYRm0exaPfPcqWE1t468Ex/PHeHCaOvI+ePWHbNi+EfESEHSpZty7kzGknJevUSUNepVjaolcud/LySZYcWELIgRCWHFjC8bDjBBesxf2bxjCsT0Eeesh2g3ulFf/PP3aJv9Wr7eia4cPtSlBKpWAa9OqOXblxhZWHVxJyIISQAyHsOm0nKc2VMRe1CtfigZsNmfhuM1b8k4a33oKBA7000eOcOdCmDdy8CZMmQfMEp3FSKsXQoFdJFhkVyeYTmwnZb4N9zZE13Iy6SfqA9FQrWI0XqrUg18U6HN30GCs/T8Mva4nui69SxQsFX78OvXrBV19BmTIwdSo8+KAXClHKOzToVaKMMRy4cCC6xb7s4DIuXrsIQOl7S9PlsR7kCavN+W3VWD0iEwM32zWy06a1Kz0NHAhvveWlVvzevXZxkC1boHt3O4zSY+sLKuUbNOhVnM5dPceyg8uiw/3QxUMAFMhWgAaFnyPv1dpc3VWLjaPv4Zstdih6unR23q/evSE42PbBZ8nixV9i0iR45RVb2OzZ8OyzXixGKe/RoFcAXIu4xpoja6K7Y7ac2ILBkC19Nqrmq0ndzD25trs2O2Y/xMRtgjFw11325qa+faFGDRvsmTJ5+zfBzk/z+uswZgxUrWoXBylQIPHnKZVCadCnUlEmih2ndkSPjvn98O+ER4STNk1ayuatxIv3DCBqX232LK3Awh1pMcYu8lGpkl3wo0YNO9GYz62et2uXvfFpzx77DjRggO1DUioV07+AVOTIpSO3DXs8c/UMAA/lKEpwto6kOViHA8trsH5bVtZjQ7xKFXvHanCw7Zbx2e5tY+CHH2w//N13w+LFULu2t6tSyido0Kdg/17/l+UHl0eH+1/n/gIgT8a8PJyuLiUu1uHIytr8veV+/sYOJ69aFVoOssFerpztnvF5ly7ZG55++cXeBDV+POTN6+2qlPIZGvQpyM3Im6w/tj462NcfXU+kiSRDQEYeSl+DSpc6cWJVHQ5vLMEZhCxZ4PHHocPHtiumbFk/XAJ140Y7qubwYRgyxA6jTKM3fCsVkwa9HzPGsOfsnuiumNBDoYTdCEMQgtKXo8Sl3pxZX4fj6yqzIzI92bLZYO/6mQ32xx7z4+7rqCj44gu7Sni+fLBypZcG6Svl+/z1zzzVOnX5FEsOLGHJwSWE7A/hWNgxAO5J+wD3XXiJjJvqcHpDTQ6E5yR7dqheHWp8artiSpVKIYsknT1r73CdNw8aNYIff7Rz1iil4qRB7+Ou3rzKysMro7tjdpzaAUCmNDnI/W8tcm6vw/mNdTh9MYiInLalXmOI/e+jj6aQYI9p5Uo7V82ZM3a5v9de08nIlEqEBr2PiYyKZMuJLdHdMauPrOZG5A3Sche5rlYl267B/LutDldPPMbVXAHUrAHBA22wFy+egrunIyNh0CA7BOiBB2DdOtv3pJRKlAa9Dzhw4QAh+0NYcnAJSw8s5cK1CwBkv1aKdH9148aOOkT88zjkzES9GlCjr+2KKVo0lTRmjx+HFi0gNNQuEjJiBGTN6u2qlPIbGvRecD78PMsOLovujjlw4QAAGW/mx+xrBH/WhgO1yJgtL/WDocbbtsX+8MOpJNhjWrAAWrWCq1ftQt2tW3u7IqX8jga9B1yPuG6nF3B0x2w6vgmDIW1kVtIcrgl7esCBOuTM8DDBNYTg122wP/hgKgz2W27etHe2fvaZvdgwdar9CKOUSjINejcwxrDz9E5C9oew+EAIKw+t5FpkOGICCDhRCfNXfzhQm3wBFQiuno4ar9iumKCgVBzsMR08aOeKX78eunSBYcN8cK4FpfyHBr2LHPv3GCEHQli8L4TF+5Zy7vopANKcf4SovR3gQG0KEUzNKtmo0da22AMDvVuzT5o+HTp0sN9PmwYvvODdepRKATTokynsehihh0JZtC+EebuXcOjKbgDkyj2Y/bXhQB0KRdWidvkCBL9kg10nUExAeDi8+SaMHGlnS5s82X7EUUrdMQ16J92MvMnG4xtZuDeE2TtD2HVxPVFEwM2McLg67G9Poaja1C31KMHPp6FGDbj/fm9X7Sf27IGmTWHHDnj7bTuM0u/mYlDKd2nQx8MYw9/n/mbB3yH8ujWEjWeXc50wMALHy8KBtykYUYcni1WmVoMMVK8O993n7ar90E8/wauv2ons58+H+vW9XZFSKY4GfQynr5xm4d9LmbIhhDUnl3CJI3bDhSDY35wCN+tQt8gT1Kubk+rV4Z57vFuvXwsLs3e1Tphgr0RPnGjnrFFKuVyqDvqrN6+ydO8qJq4NIfRICKdku90QngMOPkH+632pFVSHhjUK8/h7kDu3d+tNMbZts4uD7N9v73Tt2zcFztWglO9IVUEfZaJYe3ArY1eGsPRQCIejVmMCrkNkOvinKvmuDaJGgTq8ULUMwb0CdJ4sVzPG3tX61lv2XXPZMnuVWinlVokGvYiMARoAp40xJeLYLsBw4CngKtDGGLPFsa2eY1sAMNoY87ELa3fKrmMHGb1sCYv2hrAvcikRd523G06VJO+V16iWrw7NKj9O7Tczkz27p6tLRS5cgPbtYeZMePppe5erfkRSyiOcadGPA74BxsezvT5QxPFVEfgOqCgiAcAIoA5wFNgoInOMMX/eadEJ+efMBUYtXs68PSHsvh7C9cz77YYr+cj97zNUzFOH5hVr8UzNe8mWzZ2VqGhr19rFQU6csDc/vfGG3hmmlAclGvTGmJUiEpjALg2B8cYYA6wTkewich8QCOwzxhwAEJEpjn3dEvQXwsIJ7pmHXfdeISoNZIpMQ9UTmXjieF6ePJOF0jfvIm2ajcBGmD3YHSWo+Pz1FxQsCKtXQ/ny3q5GqVTHFX3098Ot4SmAbb3fH8/jFeM7iIh0AjoBFCxYMMlF5MiakXvD81Dlz4LUv5aXOmlzkjFtGsiJ/VLe8+ST0L+/XbRbKeVxrgj6uD6DmwQej5MxZhQwCqBcuXLx7peQReMPJudpSimVorki6I8CMW/uzw8cB+6K53GllFIe5Ir1iOYArcSqBFwyxpwANgJFRCRIRO4Cmjn2VUop5UHODK+cDAQDuUXkKNAfSAdgjBkJzMcOrdyHHV7Z1rEtQkS6AouwwyvHGGP+cMPvoJRSKgHOjLppnsh2A7wWz7b52DcCpZRSXpJSl5JWSinloEGvlFIpnAa9UkqlcBr0SimVwom9lupbROQMcDiZT88NnHVhOa7m6/WB1ugKvl4f+H6Nvl4f+FaNhYwxeeLa4JNBfydEZJMxppy364iPr9cHWqMr+Hp94Ps1+np94B81gnbdKKVUiqdBr5RSKVxKDPpR3i4gEb5eH2iNruDr9YHv1+jr9YF/1Jjy+uiVUkrdLiW26JVSSsWgQa+UUimcXwa9iNQTkb9EZJ+IvBPHdhGRrxzbd4hIGR+s8RERWSsi10Wkp6frc7LGFo7Xb4eIrBGRUj5WX0NHbdtEZJOIVPNkfc7UGGO/8iISKSIv+FJ9IhIsIpccr+E2EXnfk/U5U2OMOreJyB8issKX6hORt2O8frsc/599a107Y4xffWGnPN4PFMYubrIdKBZrn6eABdhVrioB632wxnuA8sAgoKePvo5VgByO7+t78nV0sr4s/P86U0lgj6+9hjH2W4adyfUFX6oPOwX5b57+95fEGrNj15ou6Pj5Hl+qL9b+zwDLvPV6xvfljy36CjgWHTfG3ABuLToeU/SC5caYdcCtBct9pkZjzGljzEbgpgfrismZGtcYYy44flyHXSXMl+q7bBx/XUBmEliq0ls1OnQDfgVOe7I4nK/Pm5yp8SVghjHmH7B/Oz5WX0zNgckeqSwJ/DHo41uMPKn7uJO3z++MpNbYHvspyVOcqk9EGovIHmAe0M5Dtd2SaI0icj/QGBjpwbpucfb/cWUR2S4iC0SkuGdKi+ZMjQ8BOUQkVEQ2i0grj1WXhL8TEckE1MO+qfsUV6wZ62nOLDqepIXJ3cDb53eG0zWKSE1s0HuyD9yp+owxM4GZIlIdGAjUdndhMThT45dAb2NMpEhcu7uVM/Vtwc6RcllEngJmAUXcXtn/OVNjWqAsUAvICKwVkXXGmL/dXRxJ+1t+BlhtjDnvxnqSxR+DPr7FyJO6jzt5+/zOcKpGESkJjAbqG2POeag2SOJraIxZKSIPiEhuY4ynJplypsZywBRHyOcGnhKRCGPMLF+ozxjzb4zv54vItz74Gh4FzhpjrgBXRGQlUArwRNAn5d9hM3yw2wbwy4uxaYEDQBD/vzhSPNY+T3P7xdgNvlZjjH0H4J2Lsc68jgWxawFX8dH6HuT/F2PLAMdu/ewrNcbafxyevRjrzGt4b4zXsALwj6+9hkBRYKlj30zALqCEr9Tn2O9u4DyQ2VOvXVK+/K5Fb+JZdFxEOju2x7tguS/VKCL3ApuAbECUiPTAXs3/N94De7hG4H0gF/Cto0UaYTw0U5+T9T0PtBKRm0A40NQ4/up8qEavcbK+F4AuIhKBfQ2b+dpraIzZLSILgR1AFDDaGLPLV+pz7NoYWGzspw6fo1MgKKVUCuePo26UUkolgQa9UkqlcBr0SimVwmnQK6VUCqdBr5RSKZwGvfJZIpJdRF6N8XM+EZnupnM18sbMjfFx3O4f71BWERkqIk94siblvzTolS/LDkQHvTHmuDHGXdP89gK+ddOx3eFrIN5pkZWKSYNe+bKPgQcc83x/JiKBIrILQETaiMgsEZkrIgdFpKuIvCkiW0Vk3a35wB3TIix0TIb1u4g8EvskIvIQcN04bvsXkSaOecW3O263R0QCHDVsdMyB/0qM5/cSkZ2O/T92PFbaUccOEZkpIjkcj4eKyCciskFE/haRxx2PZxSRKY79p2LndLl13nGOenaKyBsAxpjDQC7HjXdKJcjv7oxVqco72FvdSwOISGCs7SWAx4AM2LugextjHhORL4BW2AnFRgGdjTF7RaQittUeu8ujKnZyr1veB540xhwTkeyOx9oDl4wx5UUkPbBaRBYDjwCNgIrGmKsxFpwYD3QzxqwQkQ+B/kAPx7a0xpgKjknE+mMnYusCXDXGlHTML3SrntLA/caYEo7X4FY9OPapig/Olqh8iwa98mfLjTFhQJiIXALmOh7fCZQUkSzYxVOmxZg5Mn0cx7kPOBPj59XAOBH5BZjheKyu45i3uo7uxs7yWBsYa4y5CmCMOS8idwPZjTG3VkL6CZgW4/i3jrkZCHR8Xx34ynGMHSKyw/H4AaCwiHyNnYp5cYzjnAbyxfXCKBWTBr3yZ9djfB8V4+co7L/tNMDFW58IEhCODW4AjDGdHa3/p4FtIlIaO0FeN2PMophPFJF6JH0K6lt1RnL732Bc0zBfELuE45PAa8CL/H/e/QyO2pVKkPbRK18WBmRN7pMdE8QdFJEmEL2WcFzr3u7GzoSJY78HjDHrjTHvA2ex09Quwk7+lc6xz0Mikhnbwm4ndtEJRCSnMeYScOFW/zvwMpDYOqcrgRaOY5TALo2IiOQG0hhjfgX6YWfpvOUh7EyOSiVIW/TKZxljzonIascF2AXAiGQcpgXwnYi8B6TDLgW3PdY+K4FhIiKOmRs/E5Ei2Fb8Usf+O7DdLFvE9gOdARoZYxY6WvybROQGdubUPkBrYKTjDeAAic+g+h0w1tFlsw3Y4Hj8fsfjtxpl7wI43nAexM6AqlSCdPZKpQARGQ7MNcYs8XYtzhCRxkAZY0w/b9eifJ923ShlDcYuauEv0gLDvF2E8g/aoldKqRROW/RKKZXCadArpVQKp0GvlFIpnAa9UkqlcBr0SimVwv0P7CoNcsDqln0AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import context\n", "from numlabs.lab4.lab4_functions import initinter41,eulerinter41,midpointinter41\n", "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "initialVals={'yinitial': 1,'t_beg':0.,'t_end':1.,'dt':0.25,'c1':-1.,'c2':1.,'c3':1.}\n", "coeff = initinter41(initialVals)\n", "timeVec=np.arange(coeff.t_beg,coeff.t_end,coeff.dt)\n", "nsteps=len(timeVec)\n", "ye=[]\n", "ym=[]\n", "y=coeff.yinitial\n", "ye.append(coeff.yinitial)\n", "ym.append(coeff.yinitial)\n", "for i in np.arange(1,nsteps):\n", " ynew=eulerinter41(coeff,y,timeVec[i-1])\n", " ye.append(ynew)\n", " ynew=midpointinter41(coeff,y,timeVec[i-1])\n", " ym.append(ynew)\n", " y=ynew\n", "analytic=timeVec + np.exp(-timeVec)\n", "theFig,theAx=plt.subplots(1,1)\n", "l1=theAx.plot(timeVec,analytic,'b-',label='analytic')\n", "theAx.set_xlabel('time (seconds)')\n", "l2=theAx.plot(timeVec,ye,'r-',label='euler')\n", "l3=theAx.plot(timeVec,ym,'g-',label='midpoint')\n", "theAx.legend(loc='best')\n", "theAx.set_title('interactive 4.1');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In general, an *explicit* 2-stage Runge-Kutta method can be\n", "written as\n", "\n", "\n", "\n", "
eq:explicitrk1
\n", "\n", "\\begin{align}\n", "k_1 =& h f(y_n,t_n)\\\\\n", "k_2 =& h f(y_n+b_{21}k_1, t_n+a_2h) \\\\\n", "y_{n+1} =& y_n + c_1k_1 +c_2k_2\n", "\\end{align}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", " The scheme is said to be\n", "*explicit* since a given stage does not depend\n", "*implicitly* on itself, as in the backward Euler method ,\n", "or on a later stage.\n", "\n", "Other explicit second-order schemes can be derived by comparing the\n", "formula [eq: explicitrk2](#eq:explicitrk2) to the second-order Taylor method and\n", "matching terms to determine the coefficients $a_2$, $b_{21}$, $c_1$ and\n", "$c_2$.\n", "\n", "See [Appendix midpoint](#app_midpoint) for the derivation of the midpoint\n", "method." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Runge-Kutta Tableau \n", "\n", "A general s-stage Runge-Kutta method can be written as,\n", "\n", "$$\n", "\\begin{array}{l}\n", " k_i = h f(y_n+ {\\displaystyle \\sum_{j=1}^{s} } b_{ij}k_j, t_n+a_ih), \n", " \\;\\;\\; i=1,..., s\\\\\n", " y_{n+1} = y_n + {\\displaystyle \\sum_{j=1}^{s}} c_jk_j \n", "\\end{array}\n", "$$\n", "\n", "\n", "\n", "\n", "An *explicit* Runge-Kutta method has $b_{ij}=0$ for\n", "$i\\leq j$, i.e. a given stage $k_i$ does not depend on itself or a later\n", "stage $k_j$.\n", "\n", "The coefficients can be expressed in a tabular form known as the\n", "Runge-Kutta tableau. \n", "\n", "$$\n", "\\begin{array}{|c|c|cccc|c|} \\hline\n", "i & a_i &{b_{ij}} & & && c_i \\\\ \\hline\n", "1 & a_1 & b_{11} & b_{12} & ... & b_{1s} & c_1\\\\\n", "2 & a_2 & b_{21} & b_{22} & ... & b_{2s} & c_2\\\\ \n", "\\vdots & \\vdots & \\vdots & \\vdots & & \\vdots & \\vdots\\\\\n", "s &a_s & b_{s1} & b_{s2} & ... & b_{ss} & c_s\\\\\\hline\n", "{j=} & & 1 \\ 2 & ... & s & \\\\ \\hline\n", "\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "An explicit scheme will be strictly lower-triangular.\n", "\n", "For example, a general 2-stage Runge-Kutta method, \n", "\n", "\n", "$$\n", " \\begin{array}{l}\n", " k_1 = h f(y_n+b_{11}k_1+b_{12}k_2,t_n+a_1h)\\\\\n", " k_2 = h f(y_n+b_{21}k_1+b_{22}k_2, t_n+a_2h)\\\\\n", " y_{n+1} = y_n + c_1k_1 +c_2k_2\n", " \\end{array}\n", "$$\n", " \n", " \n", " has the coefficients,\n", "\n", "$$\n", "\\begin{array}{|c|c|cc|c|} \\hline\n", "i & a_i & {b_{ij}} & & c_i \\\\ \\hline\n", "1 & a_1 & b_{11} & b_{12} & c_1\\\\\n", "2 & a_2 & b_{21} & b_{22} & c_2\\\\ \\hline\n", "{j=} & & 1 & 2 & \\\\ \\hline\n", "\\end{array}\n", "$$\n", "\n", "\n", "\n", "In particular, the midpoint method is given by the tableau,\n", "\n", "$$\n", "\\begin{array}{|c|c|cc|c|} \\hline\n", "i & a_i & {b_{ij}} & & c_i \\\\ \\hline\n", "1 & 0 & 0 & 0 & 0\\\\\n", "2 & \\frac{1}{2} & \\frac{1}{2} & 0 & 1\\\\ \\hline\n", "{j=} & & 1 & 2 & \\\\ \\hline\n", "\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "##### Problem tableau \n", "\n", "Write out the tableau for\n", "\n", "1. [Heun’s method](#eq:heuns)\n", "\n", "2. the fourth-order Runge-Kutta method ([eq:rk4](#eq:rk4)) discussed in the\n", " next section." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "### Explicit Fourth-Order Runge-Kutta Method \n", "\n", "\n", "\n", "\n", "Explicit Runge-Kutta methods are popular as each stage can be calculated\n", "with one function evaluation. In contrast, implicit Runge-Kutta methods\n", "usually involves solving a non-linear system of equations in order to\n", "evaluate the stages. As a result, explicit schemes are much less\n", "expensive to implement than implicit schemes.\n", "\n", "However, there are cases in which implicit schemes are necessary and\n", "that is in the case of *stiff* sets of equations. See\n", "section 16.6 of Press et al. for a discussion. For these labs, we will\n", "focus on non-stiff equations and on explicit Runge-Kutta methods.\n", "\n", "The higher-order Runge-Kutta methods can be derived by in manner similar\n", "to the midpoint formula. An s-stage method is compared to a Taylor\n", "method and the terms are matched up to the desired order.\n", "\n", "Methods of order $M > 4$ require $M+1$ or $M+2$ function evaluations or\n", "stages, in the case of explicit Runge-Kutta methods. As a result,\n", "fourth-order Runge-Kutta methods have achieved great popularity over the\n", "years as they require only four function evaluations per step. In\n", "particular, there is the classic fourth-order Runge-Kutta formula:\n", "\n", "$$\n", " \\begin{array}{l}\n", " k_1 = h f(y_n,t_n)\\\\\n", " k_2 = h f(y_n+\\frac{k_1}{2}, t_n+\\frac{h}{2})\\\\\n", " k_3 = h f(y_n+\\frac{k_2}{2}, t_n+\\frac{h}{2})\\\\\n", " k_4 = h f(y_n+k_3, t_n+h)\\\\\n", " y_{n+1} = y_n + \\frac{k_1}{6}+ \\frac{k_2}{3}+ \\frac{k_3}{3} + \\frac{k_4}{6}\n", " \\end{array}\n", "$$\n", "\n", "\n", "\n", "
Problem RK4\n", " \n", "In the cell below, compare compare solutions to the test\n", "problem\n", "\n", "
eq:test
\n", "$$\n", "\\frac{dy}{dt} = -y +t +1, \\;\\;\\;\\; y(0) =1\n", "$$ \n", "\n", "generated with the\n", "fourth-order Runge-Kutta method to solutions generated by the forward\n", "Euler and midpoint methods.\n", "\n", "1. Based on the numerical solutions of ([eq:test](#eq:test)), which of the\n", " three methods appears more accurate?\n", "\n", "2. Again determine how the error changes relative to the change in\n", " stepsize, as the stepsize is halved." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3wU1d7H8c+PJBBKIEDoJQmI9NBCLwYRBESKil5EEBURFa9e733UR0W5VlQsICAihN4RBOmC9B4UEKSHACFAEmqABFLO88cuPBFTFtjNZDe/9+uVV3bnzM58Zwm/TM7OnCPGGJRSSnmufFYHUEop5Vpa6JVSysNpoVdKKQ+nhV4ppTycFnqllPJwWuiVUsrDaaFXlhORvSISZnWOOyEil0WkitU5lMqKFnplOWNMbWPMGkfWFZEoEXnAxZEy2/caEemffpkxpogxJtKF+3xfRExmxywiBURkvIgcE5EEEfldRDq5Ko9yT1roVZ4hNm7zMy8iVYHHgFNZrOYNnADuA4oBg4HZIhLk6nzKfbjND73yXOnP0kVkiIjMFpHJ9jPUvSISam+bAlQGfrZ3mbxhX95MRDaJyAUR2ZW+G8h+Fv6xiGwErgJVROQZEdln336kiLxwS55uIrJTRC6JyBER6SgiHwOtgZH2fY+0r2tE5B57htMi4pVuOz1EZLf9cT4Recu+vbP2YyyRzVszEngTuJ7ZCsaYK8aYIcaYKGNMmjFmEXAUaOTIe6/yBi30KjfqCswE/IGF2Aoexpg+wHHgYXuXyeciUgFYDHwElAD+A/woIqXSba8PMADwA44BsUAXoCjwDPC1iDQEEJEmwGTgf+z7bwNEGWPeAdYDg+z7HpQ+sDFmC3AFuD/d4ieB6fbH/wS6YzvzLg+cB0Zl9gaISE/gujFmiSNvWLrXlQHuBfbezuuUZ9NCr3KjDcaYJcaYVGAKUC+LdZ8CltjXTzPG/AJEAJ3TrTPRGLPXGJNijEk2xiw2xhwxNmuBFdjO1gGeA8KNMb/Yt3fSGLPfwdwzgF4AIuJnzzDD3vYC8I4xJtoYcw0YAjwmIt63bkREigCfAK85uN8br/MBpgGTbiOzygO00Kvc6HS6x1cB34wKol0g0NPebXNBRC4ArYBy6dY5kf4FItJJRLaIyDn7+p2BAHtzJeDIHeaeDjwiIgWAR4DfjDHH0uWcny7jPiAVKJPBdv4LTDHGHHV0x/bPHqZg6+YZlM3qKo/RQq/cza3DrZ7AVhT9030VNsYMzeg19iL8IzAMKGOM8QeWAJJue1Ud3PdfG435E1vXUCf+2m1zY7udbsnpa4w5mcGm2gH/tPf5n8b2y2e2iLyZ0X5FRIDx2H5pPGqMSc4qp8p7tNArd3MGSH/d+lTgYRF5UES8RMRXRMJEpGImr88PFADigBT7pYgd0rWPB54RkXb2D1AriEiNTPadkenY+uPbAHPSLR8DfCwigQAiUkpEumWyjXZAHaC+/SsGW9dPZn363wE1sX12kZhNPpUHaaFX7uZT4F17F8h/jDEngG7A29iK9wlsH6Rm+LNtjEnAVohnY/tA9ElsH/jeaN+G/QNa4CKwFlu3C8BwbP3q50VkRCb5ZgBhwK/GmPh0y4fb97NCRBKALUDTTDKeNcacvvGFrYvnvDHmMoCIvC0iS+2PA7H9EqgPnLZfEXRZRHpnkk/lQaITjyillGfTM3qllPJwWuiVUsrDaaFXSikPl22hF5FwEYkVkT2ZtHcTkd32W8YjRKRVurYoEfnjRpszgyullHJMth/Gikgb4DIw2RhTJ4P2IsAVY4wRkRBgtjGmhr0tCgi95eqDbAUEBJigoKDbeYlSSuVpO3bsiDfGlMqoLbO7DW8yxqzLaiS8G5d82RUmm5tKHBEUFEREhP4BoJRSjhKRY5m1OaWP3j5K335sg0s9m67JYLtueIeIDMhmGwPsXT8RcXFxzoillFIKJxV6Y8x8e3dNd+DDdE0tjTENsd0S/rK9GyizbYw1xoQaY0JLlcrwrw+llFJ3wKlX3Rhj1gFVRSTA/jzG/j0WmA80ceb+lFJKZS/bPvrsiMg9wBH7h7ENsY0lclZECgP5jDEJ9scdgA/udD/JyclER0eTlJR0t5HzHF9fXypWrIiPj4/VUZRSFsi20IvIjbE7AkQkGngf8AEwxowBHgX6ikgykAg8YS/6ZbANy3pjP9ONMcvuNGh0dDR+fn4EBQVh36ZygDGGs2fPEh0dTXBwsNVxlFIWcOSqm17ZtH8GfJbB8kiynjDitiQlJWmRvwMiQsmSJdEPuJXKu9zqzlgt8ndG3zel8ja3KvRKKeWppn60jE/vH4krRhTWQp/DJk6cyKBBWc/0NnHiRGJiYm4+79+/P3/++aeroymlLLJ40jaKfnydKr/5c+rAKadvXwt9LnRroR83bhy1atWyMJFSylW2rDxI4itRpOVLIfBf1ylfo7zT96GF/jZ0796dRo0aUbt2bcaOHQtAkSJFeOedd6hXrx7NmjXjzJkzAPz88880bdqUBg0a8MADD9xcfkNCQgLBwcEkJ9um97x06RJBQUHMmTOHiIgIevfuTf369UlMTCQsLOzmkBDLli2jYcOG1KtXj3bt2uXg0SulnO3wn7EcfmIzha4VpuCTh2j2/rPZv+gO3PV19FZ47TXYudO526xfH775Jut1wsPDKVGiBImJiTRu3JhHH32UK1eu0KxZMz7++GPeeOMNfvjhB959911atWrFli1bEBHGjRvH559/zpdffnlzW35+foSFhbF48WK6d+/OzJkzefTRR+nZsyejRo1i2LBhhIaG/mX/cXFxPP/886xbt47g4GDOnTvn3DdBKZVj4s5cZk27+VS+GExC59U8+sOnLtuXntHfhhEjRtw8cz9x4gSHDh0if/78dOnSBYBGjRoRFRUF2K77f/DBB6lbty5ffPEFe/fu/dv2+vfvz4QJEwCYMGECzzzzTJb737JlC23atLl5PXyJEiWceHRKqZySeDWZmU3DqXK6GqebLuLRBZ+4dH9ueUaf3Zm3K6xZs4aVK1eyefNmChUqRFhYGElJSfj4+Ny8fNHLy4uUlBQAXnnlFV5//XW6du3KmjVrGDJkyN+22bJlS6Kioli7di2pqanUqfO3UaD/whijl0oq5eZSU9MY1ehbQo815FDIXJ5fNxxc/P9az+gddPHiRYoXL06hQoXYv38/W7ZsyXb9ChUqADBp0qRM1+vbty+9evX6y9m8n58fCQkJf1u3efPmrF27lqNHjwJo141Sbujz5t8Qur8h++9dwPNbh4GXl8v3qYXeQR07diQlJYWQkBAGDx5Ms2bNslx/yJAh9OzZk9atWxMQEJDper179+b8+fP06vX/NyD369ePgQMH3vww9oZSpUoxduxYHnnkEerVq8cTTzxx9wemlMoxQx8cSfPtDdkf9AvPbxoMvr45st9sZ5iyQmhoqLl14pF9+/ZRs2ZNixK5zty5c1mwYAFTpkxx6X489f1Tyl1802cSIVMrcaTCVp5a8yQF7wl06vZFZIcxJjSjNrfso/cUr7zyCkuXLmXJkiVWR1FKudD4txZQa0Y5jpf+k0dmd3J6kc+OFnoLffvtt1ZHUEq52I8j11Hua2/ii56kzYjqlGxRP8czaKFXSikXWfPTH+R7K47EAinUeseXKk+0tySHfhirlFIusHvbCU73+wPvVG9KPnea+v/OcsR3l9JCr5RSThYddZ4dXZbjf6UEpusOwr5+1dI8WuiVUsqJLl1M4ufW06kcH8y5NivoOvO/VkfSQu8qQ4YMYdiwYVbHUErloOTkVMY3+o6a0bU53mA+T/4yzOV3vTpCC30ucWPoBKWUe0pNTeOr0OE0ONKAAzXn8czmryBf7iixuSOFm5g6dSpNmjShfv36vPDCC6SmplKkSJGb7XPnzqVfv35/e92RI0fo2LEjjRo1onXr1uzfvx+w3QH7+uuv07ZtW958882cOgyllAt83vZbmu5uyP6qS+i/6UPIn9/qSDe55+WVFoxTvG/fPmbNmsXGjRvx8fHhpZdeYtq0aQ5tesCAAYwZM4Zq1aqxdetWXnrpJX799VcADh48yMqVK/HKgfEulFKu8dlD39F8fT0OVFrLc2v/hZd/Masj/UW2hV5EwoEuQKwx5m/DK4pIN+BDIA1IAV4zxmywt3UEhgNewDhjzFAnZs9Rq1atYseOHTRu3BiAxMRESpcune3rLl++zKZNm+jZs+fNZdeuXbv5uGfPnlrklXJjXz05gcZLqnO4wlZ6L30SnwrlrI70N46c0U8ERgKTM2lfBSw0xhgRCQFmAzVExAsYBbQHooHtIrLQGHP3k59aME6xMYann36aTz/96+QA6ScTSUpK+tvr0tLS8Pf3Z2cmf4EULlzYuUGVUjlm9EuzCJlViWNl9tBjZnuK1K5mdaQMZdtHb4xZB2Q6Hq4x5rL5/5HRCgM3HjcBDhtjIo0x14GZQLe7zGuZdu3aMXfuXGJjYwHbEMHHjh2jTJky7Nu3j7S0NObPn/+31xUtWpTg4GDmzJkD2H5h7Nq1K0ezK6WcL/ydRVT9wZ9TJY7wwNgQSrZqaHWkTDnlw1gR6SEi+4HFwI1JDysAJ9KtFm1fltk2BohIhIhExMXFOSOWU9WqVYuPPvqIDh06EBISQvv27Tl16hRDhw6lS5cu3H///ZQrl/GfbNOmTWP8+PHUq1eP2rVrs2DBghxOr5Rypllfr6bsMOGc3ymaflGeSl3DrI6UJYeGKRaRIGBRRn30t6zXBnjPGPOAiPQEHjTG9Le39QGaGGNeyW5/eWmY4pyi759SzrF40nZSX4zhunci9wxOof7/PGV1JCDrYYqdenmlvZunqogEYDuDr5SuuSIQ48z9KaVUTlq3+E+SBkWRJqmUf/l8riny2bnrQi8i94h9IlMRaQjkB84C24FqIhIsIvmBfwAL73Z/Sillhd82RBHTeyf5Uwrg99QhWnz6otWRHObI5ZUzgDAgQESigfcBHwBjzBjgUaCviCQDicAT9g9nU0RkELAc2+WV4caYvS45CqWUcqGDe86wt9taAq6WJrX7Jtp9/6HVkW5LtoXeGJPl2JrGmM+AzzJpWwLo9ElKKbd18vgF1rdbSKWLgVx6YAWPzXa/Max0CASllMrEufir/NxiOkFxVYhrvojHln5hdaQ7ooVeKaUykHg1mamNx1HjZC1ONPyR3mu+zhUjUd4JLfROtHDhQoYOzXiUh/SDn92uzp07c+HChSzXmThxIjExelGTUs6QnJzK6AYjCYkK4VCdOfTbPBzceKgS9xzULJfq2rUrXbt2dfp2lyzJ/mOOiRMnUqdOHcqXL+/0/SuVl6SmpvFVo+E0PdiQA9Xn8cK2YeDjY3Wsu6Jn9A6KioqiRo0a9O/fnzp16tC7d29WrlxJy5YtqVatGtu2bWPixIkMGjQIgKNHj9K8eXMaN27M4MGDb25nzZo1tGnThh49elCrVi0GDhxIWloaADNmzKBu3brUqVPnL8MWBwUFER8fT1RUFDVr1uT555+ndu3adOjQgcTERObOnUtERAS9e/emfv36JCYm5uybo5QH+bzFcJr+0ZD9VRfRf/MHULCg1ZHumlue0b+27DV2nnbuMMX1y9bnm45ZD5Z2+PBh5syZw9ixY2ncuDHTp09nw4YNLFy4kE8++YTu3bvfXPfVV1/lxRdfpG/fvowaNeov29m2bRt//vkngYGBdOzYkXnz5tGiRQvefPNNduzYQfHixenQoQM//fTTX7YJcOjQIWbMmMEPP/zA448/zo8//shTTz3FyJEjGTZsGKGhGd4Yp5RywKf3j6T5tgbsD1rJ8xvfwKu4v9WRnELP6G9DcHAwdevWJV++fNSuXZt27dohItStW5eoqKi/rLtx40Z69bJdmdqnT5+/tDVp0oQqVarg5eVFr1692LBhA9u3bycsLIxSpUrh7e1N7969WbduXYYZ6tevD0CjRo3+tl+l1J35vNtYmq+uw8FKG3jml/54lcl+GHJ34ZZn9NmdebtKgQIFbj7Oly/fzef58uXLcCpAyeQT+luXiwiOjDl0awYvLy/tplHKCYY9EU7owns4Uj6CJ37qQYF7gqyO5FR6Ru8iLVu2ZObMmQB/m4lq27ZtHD16lLS0NGbNmkWrVq1o2rQpa9euJT4+ntTUVGbMmMF9993n8P78/PxISEhw6jEolRd83WcyDeYEElV2F92n30exhrWtjuR0WuhdZPjw4YwaNYrGjRtz8eLFv7Q1b96ct956izp16hAcHEyPHj0oV64cn376KW3btqVevXo0bNiQbt0cH76/X79+DBw4UD+MVeo2DH9+OiHTKnCi9J90mdSEkvc1tjqSSzg0THFO8+RhitesWcOwYcNYtGhRju7XU94/pZxl9KA5VPvOn1MlD3H/2FpU7B5mdaS7kmPDFCullDv4/t/zuWdMMU6XOEqbb6u6fZHPjlt+GOvOwsLCCAsLszqGUnlW+ODFBI4oSJz/CVoMK0fQEw9aHcnl9IxeKZVnTP54OeU+y8c5v9M0+qgYVZ9+2OpIOULP6JVSecKML1dT8sMULhU6R8jg/NQY+JjVkXKMntErpTze3JEb8HvnMlcLXKT6mynU+deTVkfKUVrolVIebcH4reR/4yxJ+a8Q9GoC9f/3Gasj5Tgt9Hcoq2GHL126RIUKFW4OcKaUssay6b9j/hlDSr7rVHgxjsYfvGB1JEtoob8DxpibI05mZPDgwbd1V6tSyvlW/ribqwMiAUOp547T/LNXrI5kGS30DroxRPBLL71Ew4YNb959Gh8fT/PmzVm8eDEAO3bs4MyZM3To0MHKuErlaeuW7OPCM/vxTvOmWO8DtB7+b6sjWSrbq25EJBzoAsQaY+pk0N4buDF4+mXgRWPMLntbFJAApAIpmd21dbsOvXaIyzsvO2NTNxWpX4Rq31TLcp0DBw4wYcIERo8eTZEiRThz5gxdu3blo48+on379qSlpfHvf/+bKVOmsGrVKqfmU0o5ZvPKw5zqtYtCyYXwfXwnbce+Z3UkyzlyeeVEYCQwOZP2o8B9xpjzItIJGAs0Tdfe1hgTf1cpc4nAwECaNWsGQHJyMu3atWPUqFE3u2lGjx5N586dqVSpkpUxlcqzdmyIIqrnVopeK0a+7ptpP+ljqyPlCtkWemPMOhEJyqJ9U7qnW4CKdx8ra9mdebtK4cKFbz729vamUaNGLF++/Gah37x5M+vXr2f06NFcvnyZ69evU6RIkUznkVVKOc/u7dHs77qeElcCSHloNQ/N/NzqSLmGs/vonwOWpntugBUiskNEBmT1QhEZICIRIhIRFxfn5FjOJyKEh4ezf//+m4V82rRpHD9+nKioKIYNG0bfvn21yCuVA3Zvi2ZXx18pmVCKa+1/4eF5n1kdKVdxWqEXkbbYCv2b6Ra3NMY0BDoBL4tIm8xeb4wZa4wJNcaElipVylmxXMrLy4uZM2eyevVqRo8ebXUcpfKknVtO8EfH1ZS6VJqrbZfSfdGXkMmkP3mVU4ZAEJEQYBzQyRhz9sZyY0yM/XusiMwHmgB/nx/PDQQFBbFnz56bzy9ftn0YnD9/fpYvX/639fv160e/fv1yKp5SedKODVHs77qBkpcDSHpgGY8s+UaLfAbu+oxeRCoD84A+xpiD6ZYXFhG/G4+BDsCejLeilFK3Z9vaoxx4eCMlLpckqf0vdF/ytRb5TDhyeeUMIAwIEJFo4H3AB8AYMwZ4DygJjLbPhXrjMsoywHz7Mm9gujFmmQuOQSmVx2xdHUlkjy34J/qT/OAqui/U7pqsOHLVTa9s2vsD/TNYHgnUu/NoGe4r0wm3VeZy4yxiSt2pzSsPc+yx7RRLLIbptJau87/QIp8Nt7kz1tfXl7Nnz2rRuk3GGM6ePYuvr6/VUZS6a+uXHeD4oxH4JflhuqzjofmfaZF3gNuMR1+xYkWio6Nxh0svcxtfX18qVnT57Q1KudS6Jfs53Wsnha8VxqvLBjrO1UsoHeU2hd7Hx4fg4GCrYyilLLB6wV7i++ylYHIhfLpt4cFZen/K7XCbQq+UyptW/bSHc33/xDe5AL49ttF+ug5rcLvcpo9eKZX3rJiziwt99lMgpQC+PXfQfvpHVkdyS1rolVK50rIZO0l49hA+qd4Ufvx32k/+wOpIbku7bpRSuc7iKb9x7cUovNO8KPqPPwgLH2J1JLemZ/RKqVzl58kRXB8YRT4DxZ7aS1j4YKsjuT09o1dK5RoLJmzDDIoBSaNEn4O0+e5dqyN5BD2jV0rlCvN/2AIvx5AmyZR+5ghtvnvb6kgeQwu9Uspy079cg/c/Y0n1Sqbsc8do+e2b2b9IOUy7bpRSlho3eAnlP4fLvgkED4ij8Rf/sTqSx9EzeqWUZUa++iOVhnpxoXAsNd64SuMvXrM6kkfSM3qllCW+enoaIVPLcKrEEZr9txjVXupndSSPpYVeKZXjPntkPKE/BXGi9J+0+zKYSr27WB3Jo2mhV0rlqE8e/I4WK2pypNx2uo1rTEDnMKsjeTwt9EqpHJGamsbQ+0bRcmNdDlVczxMzOlG0VajVsfIELfRKKZdLTU3jsybf0vK3ehwIXEnfn3pRsH5tq2PlGVrolVIudS0phW8ajaTFn/XZH7yY51a8iM89VayOladooVdKuczlhGt832gMTQ/V58A983h+zVt4VShvdaw8J9vr6EUkXERiRWRPJu29RWS3/WuTiNRL19ZRRA6IyGERecuZwZVSudu5+KuEh/xAo0P1OFR9Ji9sHaJF3iKO3DA1EeiYRftR4D5jTAjwITAWQES8gFFAJ6AW0EtEat1VWqWUW4g5fpGZ9ScTElWHyNrTeH77F1CihNWx8qxsC70xZh1wLov2TcaY8/anW4Abs1A3AQ4bYyKNMdeBmUC3u8yrlMrljhyIZ0mT2VSPqcbxelN5dttw8POzOlae5uwhEJ4DltofVwBOpGuLti9TSnmo3duj2dBqEYFxQcSGzqTvttFQqJDVsfI8p30YKyJtsRX6VjcWZbCayeL1A4ABAJUrV3ZWLKVUDtmyOpIjj26iTEJZLrX8kV6rx4CXl9WxFE46oxeREGAc0M0Yc9a+OBqolG61ikBMZtswxow1xoQaY0JLlSrljFhKqRyyeNrvnOgWgf9Vf1LaLeHRNaO1yOcid13oRaQyMA/oY4w5mK5pO1BNRIJFJD/wD2Dh3e5PKZW7TPp0FWn9T+CT6k3+ruvpsnQ45NOBcXOTbLtuRGQGEAYEiEg08D7gA2CMGQO8B5QERosIQIr9zDxFRAYBywEvINwYs9clR6GUssTwl+ZQ44eiXChykSpPn6LxN59ZHUllQIzJtNvcMqGhoSYiIsLqGEqpLAztMZ7GC4I4GXCIZm/5cu/r/ayOlKeJyA5jTIaDB+mdsUqp22IbnGw0LTfW4Uj5bTw8IoTSj2Z1q42ymhZ6pZTDkpKS+SZ0NC331uNA4C/0nvkIRZo1sjqWyoYWeqWUQ+LPXGFqswk0i6rHgarz6L/iX3hVCbY6lnKAFnqlVLYO7o1lzQMLCDldi8iaU3lh/WdQsqTVsZSDtNArpbK0eVUkRx7fSNDFQE43mMqzG0bq3a5uRgu9UipTCyfv4PrLkRRPLk5i63k8+ctY8Nay4W70rgalVIYmfLSSfANi8DJCoa7r6Pbrd1rk3ZQWeqXU33wzcA7l/5vK5YLnqdw/krazPwfJaPgq5Q600Cul/uLTbuOpO7YEZ4ofIXRwGo2+ecPqSOou6d9hSinAfiNU69G03FyHI+W38vCIenojlIfQQq+U4nLCNUY1/Z6W+0I4GLiC3jMepXBzvRHKU2ihVyqPO7wvjhXtf6TpyRAO3jOX55b/W2+E8jBa6JXKw1Yv2k9M3wiqXazKibqTGLDmS53b1QPph7FK5VETP1nFpcf3UzTRj6Q2c+gT8b0WeQ+lhV6pPOiLp6ZR/r1UkvInENB7Bw//+j3kz291LOUiWuiVykNSU9P4qOUoGk+rQEzAPhq/n0bzcR/oNfIeTvvolcojzsVfJbzFOFodCuFg5V/oFd4Rv3atrY6lcoAWeqXygL2/n2Jj50WEng7hcLUZPLf0TbyqVrE6lsohWuiV8nBLZu7i0sB9BF6pzKn6E+i/ZgQUK2Z1LJWDtI9eKQ/2/TtLSet3HN/k/Jh2P9Fr+3gt8nmQFnqlPNSnj06k6lAfLhWKp/Jz++m4dLSOPplHZVvoRSRcRGJFZE8m7TVEZLOIXBOR/9zSFiUif4jIThGJcFZopVTmkpNT+bjxtzSfF8TxMr/TZmghGo54W6+sycMc+fU+ERgJTM6k/RzwT6B7Ju1tjTHxtx9NKXW7zsRcYnqrybQ8WpcDQYvpO+UxCrZqanUsZbFsC70xZp2IBGXRHgvEishDTsyllLpNOzYdZ1ePldSLrUVkjSn0Xz4Er8qVrI6lcgFX99EbYIWI7BCRAVmtKCIDRCRCRCLi4uJcHEspzzLt6/UcfXAb5c6X42zjyTy7fbgWeXWTqwt9S2NMQ6AT8LKItMlsRWPMWGNMqDEmtFSpUi6OpZTn+LTnJEr9TyImXzK+Dy2n5+bxUKSI1bFULuLSQm+MibF/jwXmA01cuT+l8pKEhCQ+qTeC5nMDiS71Bw3fvkTb+d+Al5fV0VQu47JrrUSkMJDPGJNgf9wB+MBV+1MqL9m1LZot3ZfS4lQIB6vMo8/EJynYupnVsVQulW2hF5EZQBgQICLRwPuAD4AxZoyIlAUigKJAmoi8BtQCAoD5YrukyxuYboxZ5oqDUCovmTZ8PQXePUNgUiVi6oYzYOUwKF3a6lgqF3Pkqpte2bSfBipm0HQJqHeHuZRSGRj6xGQa/liWi4WTKdBpEU/+OB58fKyOpXI5vU1OKTeQkJDEt61+oMXuukSW3cF9gwpR9Z2RVsdSbkILvVK53O6Ik2zqtpgWMXU5GDyfPhOeoOB9LayOpdyIFnqlcrEZIzbg884pgpIqE1N3PAN+GQZlylgdS7kZHdRMqVxq6D+mUPL1qyCp5H/wZ56MGK9FXt0RPaNXKpe5nHCNEa3H0mJXXY6W/Y02L+Wn6uBRVsdSbkwLvVK5yO6Ik2y098cfCv6Jp8Y/RlBw0lAAABYfSURBVMG2rayOpdycFnqlcokZozbh/b8nqZIYyMk643n+ly+gbFmrYykPoH30SlksNTWNjx8Op+SrlxFS8Wn/E713jNcir5xGz+iVstDRQ2eZ12kmLY/U5kj5LbQdWJgqg7+zOpbyMFrolbLIjJGbkXePUf9SDSKrTebp6a/gHdrA6ljKA2mhVyqHpaam8WnXCTRdVplLhYWU1lN5dtEI8POzOpryUFrolcpBRw7Es6DzLFpF1uZI+c3cP7AIwe+G63yuyqW00CuVQ6aP2ES+944TklCdo9Um02/mK3g11K4a5Xpa6JVysdTUND7tEk7TFYFcLAxpbabxzELtqlE5Rwu9Ui50eF8sPz80l1ZHa3Gk/CYeeLEoge9oV43KWVrolXKRqcM34PPeSepevpeoapPoN+tVvBrUtzqWyoO00CvlZKmpaXzy0ASa/xLIhcIG2kyjn3bVKAtpoVfKiQ7ujWVJlx9pHVWTIxU20v5Ffyq/rV01ylpa6JVykslfrcf3v6eoffkeoqpNoN+sf2lXjcoVtNArdZcSEpL4puMEWmy+l3N+qch9M+i3YKR21ahcI9tBzUQkXERiRWRPJu01RGSziFwTkf/c0tZRRA6IyGERectZoZXKLRbP3M2ce2bQelNNIiuuosXb53lg1Xgt8ipXceSMfiIwEpicSfs54J9A9/QLRcQLGAW0B6KB7SKy0Bjz5x2nVSqXSE5OZWj3STRZUYEAnxLE1RnD8z99DFWrWh1Nqb/JttAbY9aJSFAW7bFArIg8dEtTE+CwMSYSQERmAt0ALfTKrW1ZE0lEn19oHV2dyPJbaPVIEvd+PRW8tSdU5U6u/MmsAJxI9zwaaJrZyiIyABgAULlyZRfGUurOGGP4/JnZ1JpZiKpplTlRfRxPT3kNr8ahVkdTKkuuLPQZXU9mMlvZGDMWGAsQGhqa6XpKWeHgnjMsfGQuTQ/V5njpPdS67yCdJv0ABQtaHU2pbLmy0EcDldI9rwjEuHB/SrnEqLeWUnbUVepfrU5k1Sn0HvEEBToPsjqWUg5zZaHfDlQTkWDgJPAP4EkX7k8ppzoTc4lxnSfTclcdThU/g3fYIp798Wvw97c6mlK3JdtCLyIzgDAgQESigfcBHwBjzBgRKQtEAEWBNBF5DahljLkkIoOA5YAXEG6M2euaw1DKuaaO2IgMOUbL83U4HDSPx95pgn//CVbHUuqOOHLVTa9s2k9j65bJqG0JsOTOoimV8y4nXOOrThNouakaF4p4caXRGPov+BIqVLA6mlJ3TK8HU8puyew/OPVKBG1ia3C40i90fC6Aiu/N0HFqlNvTQq/yvIsXEvmm+2SabwimtE8J4mp9R/95H0D16lZHU8optNCrPG3isHV4fXqc+85V50jF1bR5OIVqw6eCj4/V0ZRyGi30Kk+KOnyOaY/OoPnumpz38+FC3VE8N2MI1K5tdTSlnE4LvcpTjDF8/cpCAiel0PxyTSIDf+Lh56tQ5n9nQb5sx/hTyi1poVd5xtb1UWzst5SGkTU5WeIItPyR/jM/h4oZXjSmlMfQQq883vXrKXz2xHQaLi1J7dQqHK0ykSf+25lCvafpFTUqT9BCrzzagqm/c/p/fqP16aocLbuDek0P8eDEb/TuVpWnaKFXHulc/FVGdp9C881VKFcggFP3fkff71/CK+zfVkdTKsdpoVce54ePVlPwqxjanK/OoUoraPewUOWryVCggNXRlLKEFnrlMQ7vi2V2zzm02FubuKJCQsi3PD/9A71kUuV5WuiV20tKSmZY71nUWVaEpok1OBI0m24v1CTgjTl6yaRSaKFXbm7MkJX4joqhVXxljpXZRYnQrTw3/QsdhEypdLTQK7e0auE+dv1rDQ0jaxJX1If4Gt/y1Ncv4NXxVaujKZXraKFXbuVY5Dkm9ZpNsx1VqJkviKPBE3lkYAOKvT5LJ+dWKhP6P0O5hWtJKXzRdxa1Fxei1dV7iay8nPtaX6XTiK+gRAmr4ymVq2mhV7neDx+twnt4NK3iAzlW+g+Kh0yjf/iHULOm1dGUcgta6FWutXrxAX579VcaHalJvF8B4u79lqeGPYvXw69YHU0pt6KFXuU6J6LOE95rNs23B1M7XxBHgyfR47na+L8xU8eJV+oOaKFXuca1pBSGPTOHWgt9ue9qdQ5XWk7r5hfoOOoLCAiwOp5SbivbQi8i4UAXINYYUyeDdgGGA52Bq0A/Y8xv9rYoIAFIBVKMMaHOi648RWpqGiP+vQT/aedoGV+Z46X3ULz2VPqPHwJ161odTym358gZ/URgJDA5k/ZOQDX7V1PgO/v3G9oaY+LvIqPyUKmpaYx6axmFJ8fSIDaI2GIXias2gt5D++LVY64OIayUk2Rb6I0x60QkKItVugGTjTEG2CIi/iJSzhhzykkZlYdJTU1jzOCV+ISfJORMMHFFC3C66rd0f+V+fF+cCfnzWx1RKY/ijD76CsCJdM+j7ctOAQZYISIG+N4YM9YJ+1NuyhjD9//9FcZGUftUVc76FeJU1ZF0G9iKQoOmga+v1RGV8kjOKPQZ/X1t7N9bGmNiRKQ08IuI7DfGrMtwIyIDgAEAlStXdkIslZuM/2QNSaMPU/vkPZwrUpSTVUbT7dnGFHl9ChQsaHU8pTyaMwp9NFAp3fOKQAyAMebG91gRmQ80ATIs9Paz/bEAoaGhJqN1lPuZNGw9l0bsp+6JapwvXJwTwWPo9lQIRd8IhyJFrI6nVJ7gjEK/EBgkIjOxfQh70RhzSkQKA/mMMQn2xx2AD5ywP+UGpo7YxNkv91Dv+L1cLBjA8eCxdP1HdfzfGgtFi1odT6k8xZHLK2cAYUCAiEQD7wM+AMaYMcASbJdWHsZ2eeUz9peWAebbrr7EG5hujFnm5Pwql5k1Zisxn+2kQVR1ihYsw7GgcTzUI5CA90brPK1KWcSRq256ZdNugJczWB4J1LvzaMpdpKamEf7pWi6HH6TB0eoU8i3HsaBwOnUpS+khw6FkSasjKpWn6Z2x6o5dvJDIqFcWUXb5VarFBXKpYFmOBk2iUwd/yn44DEqXtjqiUgot9OoO7N0Zw+zXl1Fvmz8trpTiVPEoTgV/S+du91LsjU+gXDmrIyql0tFCrxy2YOrv7PlsK6H7g2mbUoXIclvxrbiWx//ZFa9+U6BQIasjKqUyoIVeZSk5OZUx7/6CmXGCkBPVCPUOJrr8CuqWP8Kz770AD/6PTsCtVC6nhV5l6EzMJb5/dRHBqwx1z1fgfOHiHAsaz/2h+Xnwvdd0sDGl3IgWevUXW9dHsfStlTT8rTRtksoTHXCQ+OCv6NSrAYVf/VI/YFXKDWmhVyQlJfPDf1dxdV4UjQ7fQ5u0KhytsIFCRTfS61+P49Vnto5Do5Qb00Kfh/00+Td2jt5B3T0lqXulBAm+5Tle8ScaB53k/sEvQ7vBOlSwUh5AC30eszsimnn/XU2lbWlUjQ2kVb5gjpfdineZDbTrXIsiLw+CGjWsjqmUciIt9HnAufirjBu8Au9lsdQ7VpUwU4nogAOcDBxFq7o+PPBSH2j/Jnjrj4NSnkj/Z3uo1NQ0pn6zkeNT/qT+/go0uebPucLJHK88mzoljvPUwB7w+A86/oxSeYAWeg+zbtkhVn2xgXt/8yXwQjnKeAcSXW49Pvm30aVnE/I/+x+oVs3qmEqpHKSF3s2lpqaxaPpOfpv6B2X3pFIzpgptCSaqzE5ig6bTtklJOg58Gu77UG9sUiqP0kLvhuLOXGbGl+u5sPIE1Y/4U+ZSadoSyMmShzgWGE7jCucIG/gE9Jihk3sopbTQu4tNvx5h5ehtFPntMrVOVCIkpSBJPpU4Weo38J9Gg0qphD3eAbp/BDoVo1IqHS30uVRSUjIzRmwmasFBgg4WIDi+Em0oR1zRGE6WX0J5rz20ahGE3yNdoP0/wc/P6shKqVxKC30usu+PUywYvgmz8Sy1jpUlOLEolSWYE6X/ILryYmoVi6N1t5Z4de0DjRppn7tSyiFa6C1y8UIii6b8zuFVkfgcvErgKT8qXChHM0pyqaAXcSW3YPJF0KxmEdo91hk6fwXly1sdWynlhrTQ54DU1DTWLT/I1nl7ubbrHGWjfQiOK0+F1PxUoCLnC5/lfLF9nCi2iEDvE3TqEIJP1y4Q9p6OMaOUumta6F3gyIE4lk35ndjNJyl6NIWqZ0rhf9WfZpTkmncRTpc4SHSFnwgwB6heCtq0rEu+Fs2g2VMQGKjjyyilnEoL/R26fj2FHRuPsWfzcWL/jOP68cv4xqVQMbYIlc6Vpzb5qU0wp/2PcbbkFq4F7CMw/1maNaiMb6tm0Kwf1KsHBQpYfShKKQ+XbaEXkXCgCxBrjKmTQbsAw4HOwFWgnzHmN3tbR3ubFzDOGDPUidld7uiReLavjuT4ztNcjryAnL5G0XNCqUuFKX2pJPlT81MNoRqlSZWSnC8SzyW/IxyvvJTyJpKQID/C2oRC05bQ9HUdy10pZQlHzugnAiOByZm0dwKq2b+aAt8BTUXECxgFtAeige0istAY8+fdhnbU9espxMZcJvZUAudiL3Mh7goJ5xO5ej6RpEtJXLt0jZQryaRdTcEkpiCJafheNJS4kJ8yl4pTLLEopYHSFAWKcrlAAhf8TnGt0AFOFj1F4bRTFE+Lo5yfcE9QSXzvrQIhIdCsm20ESC+vnDpUpZTKVLaF3hizTkSCslilGzDZGGOALSLiLyLlgCDgsDEmEkBEZtrXdVmhn1xqHAWvFyF/SgEKXC+Ib8r/f5DpDQTYvyC//euvkryTuFT4LFcKxhAfsIsEOYV/SgylfK4QXLowpe+tBFWrQpUqUPUB2/eyZfUyR6VUruaMPvoKwIl0z6PtyzJa3jSzjYjIAGAAQOU7vLPzkl8slzgDXtcxXsnky5eMl1cy3l4p+HilUMArBV+vVAp6p1LYKw0/nzSKeqfh75NGUZ80vPKJrXulalWo0sb2PThYhxFQSrk1ZxT6jC4RMVksz5AxZiwwFiA0NDTT9bIyKPLtO3mZUkp5NGcU+migUrrnFYEYbH0jGS1XSimVg5zRubwQ6Cs2zYCLxphTwHagmogEi0h+4B/2dZVSSuUgRy6vnAGEAQEiEg28D/gAGGPGAEuwXVp5GNvllc/Y21JEZBCwHNvlleHGmL0uOAallFJZcOSqm17ZtBvg5UzalmD7RaCUUsoiel2gUkp5OC30Sinl4bTQK6WUh9NCr5RSHk5sn6XmLiISBxy7w5cHAPFOjONu9Pj1+PX486ZAY0ypjBpyZaG/GyISYYwJtTqHVfT49fj1+PPu8WdGu26UUsrDaaFXSikP54mFfqzVASymx5+36fGrv/G4PnqllFJ/5Yln9EoppdLRQq+UUh7OLQu9iHQUkQMiclhE3sqgXURkhL19t4g0tCKnqzhw/L3tx71bRDaJSD0rcrpSdu9BuvUai0iqiDyWk/lczZHjF5EwEdkpIntFZG1OZ3QlB/4PFBORn0Vkl/34n7EiZ65hjHGrL2xDHh8BqmCb3GQXUOuWdToDS7HNctUM2Gp17hw+/hZAcfvjTp50/I6+B+nW+xXbCKqPWZ07h38G/LHNz1zZ/ry01blz+PjfBj6zPy4FnAPyW53dqi93PKNvgn3ScWPMdeDGpOPp3Zyw3BizBbgxYbknyPb4jTGbjDHn7U+3YJvdy5M48jMA8ArwIxCbk+FygCPH/yQwzxhzHMAY40nvgSPHbwA/ERGgCLZCn5KzMXMPdyz0mU1GfrvruKvbPbbnsP1140myfQ9EpALQAxiTg7lyiiM/A/cCxUVkjYjsEJG+OZbO9Rw5/pFATWzTl/4BvGqMScuZeLmPM+aMzWmOTDp+WxOTuxmHj01E2mIr9K1cmijnOfIefAO8aYxJtZ3UeRRHjt8baAS0AwoCm0VkizHmoKvD5QBHjv9BYCdwP1AV+EVE1htjLrk6XG7kjoU+s8nIb3cdd+XQsYlICDAO6GSMOZtD2XKKI+9BKDDTXuQDgM4ikmKM+SlnIrqUo/8H4o0xV4ArIrIOqAd4QqF35PifAYYaWyf9YRE5CtQAtuVMxNzFHbtuHJl0PLMJyz1BtscvIpWBeUAfDzmDu1W274ExJtgYE2SMCQLmAi95SJEHx/4PLABai4i3iBQCmgL7cjinqzhy/Mex/TWDiJQBqgOROZoyF3G7M3qTyaTjIjLQ3p7phOWewMHjfw8oCYy2n9GmGA8a0c/B98BjOXL8xph9IrIM2A2kAeOMMXusS+08Dv77fwhMFJE/sHX1vGmMyavDF+sQCEop5encsetGKaXUbdBCr5RSHk4LvVJKeTgt9Eop5eG00CullIfTQq9yLRHxF5GX0j0vLyJzXbSv7iLyniu2fSfsQxdkekmsiAwTkftzMpNyX1roVW7mD9ws9MaYGGOMq4YbfgMY7aJtu8K3QKbDMyuVnhZ6lZsNBarax1T/QkSCRGQPgIj0E5Gf7GOOHxWRQSLyuoj8LiJbRKSEfb2qIrLMPrDXehGpcetORORe4NqNG2pEpKeI7LGPZb7OvszLnmG7fZz/F9K9/g0R+cO+/lD7svr2HLtFZL6IFLcvXyMin4nINhE5KCKt7csLishM+/qzsI1Pc2O/E+15/hCRfwEYY44BJUWkrKvefOU53O7OWJWnvAXUMcbUBxCRoFva6wANAF9sd0G/aYxpICJfA32xDWw2FhhojDkkIk2xnbXf2uXREvgt3fP3gAeNMSdFxN++7DlsQ2k0FpECwEYRWYFt/JTuQFNjzNUbv2CAycArxpi1IvIB8D7wmr3N2xjTREQ625c/ALwIXDXGhNjHKbqRpz5QwRhTx/4e3MiDfZ2W2IZiVipTWuiVO1ttjEkAEkTkIvCzffkfQIiIFME2CcucdCNYFshgO+WAuHTPN2K7fX42tjGDADrYt3mj66gYUA1bkZ5gjLkKYIw5JyLFAH9jzI1ZnSYBc9Jt/8Y2dwBB9sdtgBH2bewWkd325ZFAFRH5FlgMrEi3nVigfEZvjFLpaaFX7uxausdp6Z6nYfvZzgdcuPEXQRYSsRVuAIwxA+1n/w8BO0WkPrbxUl4xxixP/0IR6cjtD4F9I2cqf/0/+LftGGPOi20qyAeBl4HHgWftzb727EplSfvoVW6WAPjd6YvtY48fFZGecHMu4Yzmz90H3HPjiYhUNcZsNca8B8RjGxJ3OfCiiPjY17lXRApjO8N+1j5CJCJSwhhzETh/o/8d6ANkN2frOqC3fRt1gBD74wAgnzHmR2AwkH7+43sBjxioTLmWntGrXMsYc1ZENto/gF0KjLqDzfQGvhORdwEfbNPO7bplnXXAlyIi9vHLvxCRatjO4lfZ19+NrZvlN7H1A8UB3Y0xy+xn/BEich3byKlvA08DY+y/ACLJfgTV74AJ9i6bnfz/uOkV7MtvnJT9L4D9F849QMTtviEq79HRK5UCRGQ48LMxZqXVWRwhIj2AhsaYwVZnUbmfdt0oZfMJUMjqELfBG/jS6hDKPegZvVJKeTg9o1dKKQ+nhV4ppTycFnqllPJwWuiVUsrDaaFXSikP9389kQaaD2K7MgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from numlabs.lab4.lab4_functions import initinter41,eulerinter41,midpointinter41,\\\n", " rk4ODEinter41\n", "initialVals={'yinitial': 1,'t_beg':0.,'t_end':1.,'dt':0.05,'c1':-1.,'c2':1.,'c3':1.}\n", "coeff = initinter41(initialVals)\n", "timeVec=np.arange(coeff.t_beg,coeff.t_end,coeff.dt)\n", "nsteps=len(timeVec)\n", "ye=[]\n", "ym=[]\n", "yrk=[]\n", "y=coeff.yinitial\n", "ye.append(coeff.yinitial)\n", "ym.append(coeff.yinitial)\n", "yrk.append(coeff.yinitial)\n", "for i in np.arange(1,nsteps):\n", " ynew=eulerinter41(coeff,y,timeVec[i-1])\n", " ye.append(ynew)\n", " ynew=midpointinter41(coeff,y,timeVec[i-1])\n", " ym.append(ynew)\n", " ynew=rk4ODEinter41(coeff,y,timeVec[i-1])\n", " yrk.append(ynew)\n", " y=ynew\n", "analytic=timeVec + np.exp(-timeVec)\n", "theFig=plt.figure(0)\n", "theFig.clf()\n", "theAx=theFig.add_subplot(111)\n", "l1=theAx.plot(timeVec,analytic,'b-',label='analytic')\n", "theAx.set_xlabel('time (seconds)')\n", "l2=theAx.plot(timeVec,ye,'r-',label='euler')\n", "l3=theAx.plot(timeVec,ym,'g-',label='midpoint')\n", "l4=theAx.plot(timeVec,yrk,'m-',label='rk4')\n", "theAx.legend(loc='best')\n", "theAx.set_title('interactive 4.2');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Embedded Runge-Kutta Methods: Estimate of the Truncation Error \n", "\n", "\n", "\n", "It is possible to find two methods of different order which share the\n", "same stages $k_i$ and differ only in the way they are combined, i.e. the\n", "coefficients $c_i$. For example, the original so-called embedded\n", "Runge-Kutta scheme was discovered by Fehlberg and consisted of a\n", "fourth-order scheme and fifth-order scheme which shared the same six\n", "stages.\n", "\n", "In general, a fourth-order scheme embedded in a fifth-order scheme will\n", "share the stages \n", "\n", "$$\n", " \\begin{array}{l}\n", " k_1 = h f(y_n,t_n)\\\\\n", " k_2 = h f(y_n+b_{21}k_1, t_n+a_2h)\\\\\n", " \\vdots \\\\\n", " k_6 = h f(y_n+b_{51}k_1+ ...+b_{56}k_6, t_n+a_6h)\n", " \\end{array}\n", "$$\n", "\n", " \n", "\n", "\n", "\n", "\n", "The fifth-order formula takes the step: \n", "\n", "$$\n", " y_{n+1}=y_n+c_1k_1+c_2k_2+c_3k_3+c_4k_4+c_5k_5+c_6k_6\n", "$$ \n", "\n", "while the\n", "embedded fourth-order formula takes a different step:\n", "\n", "\n", "\n", "$$\n", " y_{n+1}^*=y_n+c^*_1k_1+c^*_2k_2+c^*_3k_3+c^*_4k_4+c^*_5k_5+c^*_6k_6\n", "$$\n", "\n", "If we now take the difference between the two numerical estimates, we\n", "get an estimate $\\Delta_{\\rm spec}$ of the truncation error for the\n", "fourth-order method, \n", "\n", "\n", "$$\n", " \\Delta_{\\rm est}(i)=y_{n+1}(i) - y_{n+1}^{*}(i) \n", "= \\sum^{6}_{i=1}(c_i-c_{i}^{*})k_i\n", "$$ \n", "\n", "This will prove to be very useful\n", "in the next lab where we provide the Runge-Kutta algorithms with\n", "adaptive stepsize control. The error estimate is used as a guide to an\n", "appropriate choice of stepsize.\n", "\n", "An example of an embedded Runge-Kutta scheme was found by Cash and Karp\n", "and has the tableau: " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\n", "\\begin{array}{|c|c|cccccc|c|c|} \\hline\n", "i & a_i & {b_{ij}} & & & & & & c_i & c^*_i \\\\ \\hline\n", "1 & & & & & & & & \\frac{37}{378} & \\frac{2825}{27648}\\\\\n", "2 & \\frac{1}{5} & \\frac{1}{5}& & & & & & 0 &0 \\\\\n", "3 & \\frac{3}{10} & \\frac{3}{40}&\\frac{9}{40}& & & & &\\frac{250}{621}&\\frac{18575}{48384}\\\\\n", "4 & \\frac{3}{5}&\\frac{3}{10}& -\\frac{9}{10}&\\frac{6}{5}& & & &\\frac{125}{594}& \\frac{13525}{55296}\\\\\n", "5 & 1 & -\\frac{11}{54}&\\frac{5}{2}&-\\frac{70}{27}&\\frac{35}{27}& & & 0 & \\frac{277}{14336}\\\\\n", "6 & \\frac{7}{8}& \\frac{1631}{55296}& \\frac{175}{512}&\\frac{575}{13824}& \\frac{44275}{110592}& \\frac{253}{4096}& & \\frac{512}{1771} & \\frac{1}{4}\\\\\\hline\n", "{j=} & & 1 & 2 & 3 & 4 & 5 & 6 & & \\\\ \\hline\n", "\\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "##### Problem embedded\n", "\n", "Though the error estimate is for the embedded\n", "fourth-order Runge-Kutta method, the fifth-order method can be used in\n", "practice for calculating the solution, the assumption being the\n", "fifth-order method should be at least as accurate as the fourth-order\n", "method. In the demo below, compare solutions of the test problem\n", "[eq:test2](#eq:test2]) \n", "\n", "
eq:test2
\n", "$$\\frac{dy}{dt} = -y +t +1, \\;\\;\\;\\; y(0) =1$$\n", "\n", "generated by the fifth-order method with solutions generated by the\n", "standard fourth-order Runge-Kutta method. Which method\n", "is more accurate? Again, determine how the error decreases as you halve\n", "the stepsizes. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dd3hUVf7H8fcXCCBFIAQBaQEiiiAgRJqiUYqICkhRWWyoq65l11XXdW3gqqvr6trLIgKiP7EC6lIElSIiShCkIyG00ENJ6Ibk/P6YGzZiOpncmcnn9Tx5mLn35N7PTIZvTm45x5xziIhI5CrndwAREQkuFXoRkQinQi8iEuFU6EVEIpwKvYhIhFOhFxGJcCr04jszW25mCX7nKA4z229mzfzOIZIfFXrxnXOulXNuVmHamtl6M+sR5Eh57XuWmd2cc5lzrppzLjmI+xxuZi6/12xm75rZVjNLN7Ofj88ookIvZYYFhM1n3syaA4OArQU0fQqIdc6dDPQFnjCzDsHOJ+EjbD70Erly9tLNbISZfWhm48xsn3dYJ95b9w7QGPjcO2Ryv7e8s5nNM7O9ZvZTzsNAXi/8STP7FjgINDOzYWa20tt+spndelyefma22OshrzWz3mb2JNANeMXb9yteW2dmcV6GbWZWPsd2rjCzJd7jcmb2gLe9Xd5rjC7grXkF+CvwS36NnHPLnXNHsp96X80L2LaUJc45fenL1y9gPdDDezwCOAz0AcoT6K3Oz62t97wBsMtrXw7o6T2v462fBWwEWgEVgCjgUgKF0IALCPwCaO+17wikedsp523/jBzbuvm47A6I8x6vBXrmWPcR8ID3+G5gPtAQqAT8Bxifz3syGPg0t9ecR/vXvNfhgB+Ban7/XPUVOl/q0Usomuucm+KcywTeAdrm0/YaYIrXPss5NwNIJFD4s411gV7vUedchnNusnNurQuYDUwn0FsHuAkY7Zyb4W1vs3NuVSFzjweGAJhZdS/DeG/drcBDzrkUF+h9jwAGmVmF4zdiZtWAfxD45VAozrnbgere65gAHMn/O6QsUaGXULQtx+ODQOXcCqKnCTDYO2yz18z2AucB9XO02ZTzG8zsEjObb2a7vfZ9gBhvdSMCPfPieA8YYGaVgAHAj865DTlyTsyRcSWQCdTNZTuPAe8459YVZefOuUzn3FwCfzX8oZivQSKQCr2Em+OHW91EoCjWzPFV1Tn3dG7f4xXhT4BngbrOuZrAFAKHcbK3l9fx7XyHenXOrQA2AJcAvyNQ+HPmvOS4nJWdc5tz2VR34I/eMf9tBH75fGhmf81v/zlUyOc1SBmkQi/hZjuQ87r1d4HLzexiMytvZpXNLMHMGubx/RUJHCPfCRw1s0uAXjnWvwUMM7Pu3gnUBmZ2Rh77zs17wB+B8wkco8/2BvCkmTUBMLM6ZtYvj210B1oD7byvLQQO/bx6fEMzO8XMrjazat7rv5jA4aOvC8gpZYgKvYSbp4CHvUMg9znnNgH9gAcJFO9NwF/I47PtnNtHoBB/COwh0PP+LMf6H4BhwPMETsrOJnDYBeBFAsfV95jZS3nkGw8kAF8751JzLH/R2890M9tH4MRspzwy7nLObcv+InCIZ49zbj+AmT1oZlOzmxM4TJPivZ5ngbudc5/mkU/KIHNOE4+IiEQy9ehFRCKcCr2ISIRToRcRiXAq9CIiES6vm1B8FRMT42JjY/2OISISNhYuXJjqnKuT27qQLPSxsbEkJib6HUNEJGyY2Ya81unQjYhIhFOhFxGJcCr0IiIRLiSP0ecmIyODlJQUDh8+7HeUkFa5cmUaNmxIVFSU31FEJESETaFPSUmhevXqxMbGYmYFf0MZ5Jxj165dpKSk0LRpU7/jiEiIKPDQjZmNNrMdZrYsj/VDzWyJ9zXPzNrmWLfezJZ607Kd0GU0hw8fpnbt2iry+TAzateurb96RORXCnOMfizQO5/164ALnHNtgMeBkcetv9A51845F1+8iP+jIl8wvUcicrwCC71zbg6wO5/185xze7yn2XNiiohIEbz7/Dweu/TdoGy7pK+6uQmYmuO5IzD+9kIzuyW/bzSzW8ws0cwSd+7cWcKxQsfYsWO58847C2yzZcuWY89vvvlmVqxYEexoIuKTj0YmUvNve4ibF0XK+r0lvv0SOxlrZhcSKPTn5Vh8rnNui5mdAswws1XeXwi/4ZwbiXfYJz4+vkwPkj927Fhat27NqaeeCsCoUaN8TiQiwfLf95ZQ8e4tHKj4C00/OoeGsTVLfB8l0qM3szbAKKCfc25X9nLn3Bbv3x3ARKBjSezPL/3796dDhw60atWKkSMDpyKqVavGQw89RNu2bencuTPbt28H4PPPP6dTp06cffbZ9OjR49jybPv27aNp06ZkZGQAkJ6eTmxsLB999BGJiYkMHTqUdu3acejQIRISEo4NCTFt2jTat29P27Zt6d69eym+ehEpaV9/torDt6zlaLlMose1omuPgmaqLJ4T7tGbWWNgAnCtc+7nHMurAuWcc/u8x72Av5/o/gDuvhsWLy6JLf1Pu3bwwgv5txk9ejTR0dEcOnSIc845h4EDB3LgwAE6d+7Mk08+yf3338+bb77Jww8/zHnnncf8+fMxM0aNGsUzzzzDc889d2xb1atXJyEhgcmTJ9O/f3/ef/99Bg4cyODBg3n11Vd59tlniY//9fnrnTt38vvf/545c+bQtGlTdu/O89SJiIS4eV8ms/OaZVTOqki5kU3o2b9l0PZVYKE3s+w5MGPMLAUYDkQBOOfeAB4FagOveVd8HPWusKkLTPSWVQDec85NC8JrKDUvvfQSEydOBGDTpk2sWbOGihUrctlllwHQoUMHZsyYAQSu+7/qqqvYunUrv/zyS67Xtd98880888wz9O/fnzFjxvDmm2/mu//58+dz/vnnH9tWdHR0Sb48ESkli+ensG7wAmr8UpWDL9TlymvaFvxNJ6DAQu+cG1LA+puBm3NZngwEJX1BPe9gmDVrFl9++SXfffcdVapUISEhgcOHDxMVFXXsksby5ctz9OhRAO666y7uuece+vbty6xZsxgxYsRvtnnuueeyfv16Zs+eTWZmJq1bt843g3NOl0+KhLlVS7az+LI51D1Qix1PVef6284J+j411k0hpaWlUatWLapUqcKqVauYP39+ge0bNGgAwNtvv51nu+uuu44hQ4YwbNiwY8uqV6/Ovn37ftO2S5cuzJ49m3Xr1gHo0I1ImFmftJu5PWdQP602Gx+qxPX3nlfwN5UAFfpC6t27N0ePHqVNmzY88sgjdO7cOd/2I0aMYPDgwXTr1o2YmJg82w0dOpQ9e/YwZMj//nC64YYbuO22246djM1Wp04dRo4cyYABA2jbti1XXXXVib8wESkV2zalM/X8z2mSWpfVdztuHX5Rqe3bnAu9Kxnj4+Pd8ROPrFy5kpYtg3eywi8ff/wxn376Ke+8806JbTNS3yuRcLU79SDvnv0BrTY3ZtFtB7jvtb4lvg8zW5jXCARhM6hZJLrrrruYOnUqU6ZM8TuKiARJetphxnQcT4eU5sy7LpUHXxtU6hlU6H308ssv+x1BRILoyKGjvH7Ou3RaF8ecAVt59O18r20JGh2jFxEJgoyMTJ47Zyyd1sQxq88mHv3EnyIPKvQiIiUuMzOLpzuOoevyOGYmrGfE5Gt9zaNCLyJSgjIzs3ji3LF0WxzHzM5refTL6/yOpEIvIlKSHu8xjgu+b8acs5N4dO4wypf3v8z6nyBMVatWrdBtcw5KJiKR67E+75IwK5ZvW63hb9+HRpEHFfpicc6RlZXldwwRCSGPD3qfC6Y2ZH6LJO5ZcANRUeX9jnSMCn0hrV+/npYtW3L77bfTvn37Y3espqam0qVLFyZPngzAM888w1lnnUXbtm154IEHfrWNrKwsrr/+eh5++OFSzy8iwfPUtZ/Q7ZN6LGy6ltt/uIaTToryO9KvhOV19HdPu5vF20p2nOJ29drxQu/8R0tbvXo1Y8aM4bXXXqNatWps376dvn378sQTT9CzZ0+mTp3KpEmT+P7776lSpcqvxqI5evQoQ4cOpXXr1jz00EMlml1E/PPc7Z/T8f9q8lOjZIb9MISTa1T2O9JvqEdfBE2aNDk2xk1GRgbdu3fnmWeeoWfPngB8+eWXDBs2jCpVqgC/Hkb41ltvVZEXiTAv/2Uabf9zEqvrbWLwtwOJjqnid6RchWWPvqCed7BUrVr12OMKFSrQoUMHvvjiCy644AIg/2GEu3btysyZM7n33nupXDn0fuOLSNGMfOxrTn++HOtittJnTl9ObVTD70h5Uo++mMyM0aNHs2rVKp5++mkAevXqxejRozl48CDw62GEb7rpJvr06cPgwYOPjVkvIuHp7Wfn0ujJI2yumcp5M3oSGxfakwCp0J+A8uXL8/777zNz5kxee+01evfuTd++fYmPj6ddu3Y8++yzv2p/zz330L59e6699lpdtSMSpj58YwG1H05jV9V0zp58Pi3b1PM7UoE0THEE0nslEhyfv/sTWbds4FDUERpP6EDX7sGZzLs4NEyxiMgJ+urTlfxyWzKuHNR556yQKvIFUaEXESnA3BlJ7Lp2BRWzKlDhraZ073uG35GKRMfoRUTy8eN3m9gw+EdO+qUSGS824LIhbfyOVGQq9CIieVj50zaWXjaX6IPV2fNULQb/PtdD4CFPhV5EJBfJP6fybc8vqZcezaZHKnPdn8/1O1KxqdCLiBxny8Y0vkiYTJNddfn5z45bHrnQ70gnRIW+mIoyTLGIhI/U7Qf45NwJtNjWkCW3HeauZ3r7HemEqdAXg4YpFolM6WmHGdf5A1qlNGHBdWnc++rlfkcqESr0hVTYYYpFJDwdOpjBG/Hv0n59M74dvIMHxg7wO1KJCcvr6NfcvYb9i/eX6DartavGaS+clm+bgoYpFpHwlJGRyfPnvE3XpDhmXZrCiA+v8TtSiQrLQu+X3IYpfvXVV4+NXiki4SczM4unO46h24o4Zl60nsf+e4PfkUpcWBb6gnrewVLQMMUiEl4yM7N4ousYLlgcx8wua3l0+jC/IwVFgcfozWy0me0ws2V5rB9qZku8r3lm1jbHut5mttrMkszsgdy+P1zlNkyxiISXx7uP44IfmjO7fRKPfhM6k3mXtMK8qrFAftcXrQMucM61AR4HRgKYWXngVeAS4ExgiJmdeUJpQ8zxwxSLSPh47JJ3SJgdy9xWa3hwfuQWeSjEoRvn3Bwzi81n/bwcT+cDDb3HHYEk51wygJm9D/QDVhQ3rJ9iY2NZtux/f9Ts3x84GVyxYkW++OILv2KJSDE8PmA8F0xrxPwWSfwlcRhRUeX9jhRUJf0r7CZgqve4AbApx7oUb1muzOwWM0s0s8SdO3eWcCwRkYCnrv2YbhPrk9hsLXcmXkulymF5qrJISqzQm9mFBAr9X7MX5dIsz1lOnHMjnXPxzrn4OnXqlFQsEZFjnr3tMzr+Xy0WN07mxu+HUK16Jb8jlYoS+VVmZm2AUcAlzrld3uIUoFGOZg2BLSeyn/wm35aAUJwxTCQUvHjPVNq9WYVV9Tdx9beDiI6p4nekUnPCPXozawxMAK51zv2cY9UC4DQza2pmFYGrgc+Ku5/KlSuza9cuFbJ8OOfYtWsXlStX9juKSEj5z/CvaPlSeZLrbOXSOX2p1/BkvyOVqgJ79GY2HkgAYswsBRgORAE4594AHgVqA695ve2j3iGYo2Z2J/AFUB4Y7ZxbXtygDRs2JCUlBR2/z1/lypVp2LBhwQ1Fyoix//qGxk9lsLnWLi6Y0YvY5tF+Ryp1YTM5uIhIUX3w2g9UvWcne6rs46xp59GuY+R2gjQ5uIiUOZ++s5jK925nX6XDnPZJp4gu8gVRoReRiDNjwgoy/7CerApZ1H23DZ0vbOp3JF9F7q1gIlImzZ2exJ4bVmLOqPxmHBddfrrfkXynQi8iESNx7kY2XvkjJ/1SiawXG3DZ1W38jhQSVOhFJCIsX7SV5f2+peahaux9OpqBN+d6XrJMUqEXkbCXvDqV+Rd/Rb30aDY/UoVr7+7qd6SQokIvImFt84Y0pidMofGuuqy51/j9wwl+Rwo5KvQiErZStx9gwnkTOW17A5befoQ7n+7ld6SQpEIvImEpfe9hxnX+gFYpjVlwQxr3vHyZ35FClgq9iISdQwczeOOc/6P9+mbMu3InD4we4HekkKZCLyJhJSMjkxfi36ZjUnNmX7aZhz+4yu9IIU+FXkTCRmZmFk/Hj6HLyjhmXrSe4Z8P9TtSWFChF5GwkJmZxRNdxtBtSRwzuyYzfMb1fkcKGyr0IhIWHr9oHBcsaM7sDkk8OucGypXTJESFpUIvIiFvxMXvkDAnlm/OWsOD3w2jfHmVrqLQuyUiIe3vV4wnYXoj5p+exP0/DCMqqrzfkcKOCr2IhKx/XPMx50+qz4JmSdy54FoqVdbI6sWhQi8iIemZWybR6b1aLG6czE0//I5q1Sv5HSlsqdCLSMh58c9TaP9WNVbV38TV8wYRXbuK35HCmgq9iISUNx75ipYvV2BtnS1c/k0/6jU42e9IYU8HvEQkZIz95zfEPp1BSq1dXPhVbxo3q+V3pIigQi8iIeH9V7+nzqP72Fk9nfjJCbRodYrfkSKGCr2I+G7S2EVUuW8H6ScdosXEzrTp2MDvSBFFhV5EfDX94xW4OzZwpEIW9d9rR6cLYv2OFHF0MlZEfPPNtDWkDVsFGNXeakFCnxZ+R4pIKvQi4ovEbzaw6epFVDoahXu5IZdc2drvSBFLhV5ESt3yH7eyov931DxUjfR/xjDgxg5+R4poKvQiUqqSVu5k/sVfc0p6TbYMr8I1f+zid6SIp0IvIqUmZf1evrxoKo321GHtX8px84MJfkcqE1ToRaRU7Ni6j0nnTeK07Q1Ydscv3PGPXn5HKjNU6EUk6PbuPsR7XT7izM2NWTAsnXtevMzvSGVKgYXezEab2Q4zW5bH+jPM7DszO2Jm9x23br2ZLTWzxWaWWFKhRSR8HDzwC292fI92G5ox7+qdPPDWFX5HKnMK06MfC/TOZ/1u4I/As3msv9A51845F1/EbCIS5o4cPsqL8eM4Z21zZvfdzMPjr/I7UplUYKF3zs0hUMzzWr/DObcAyCjJYCIS3jIzs3im01i6rIpjZo8NDP90qN+RyqxgH6N3wHQzW2hmt+TX0MxuMbNEM0vcuXNnkGOJSDDt33eEf3QYTbclccw8N5nHZlzvd6QyLdhj3ZzrnNtiZqcAM8xslfcXwm8450YCIwHi4+NdkHOJSJCsWrKdGZdOoVtKHDO7JvPo7Bv8jlTmBbVH75zb4v27A5gIdAzm/kTEX5+OW8yi8+dw+tYGfHftLh779kbKl9fFfX4L2k/AzKqaWfXsx0AvINcrd0Qk/D3/p8lUvHkrUUcrsP3Zavxt3EC/I4mnwEM3ZjYeSABizCwFGA5EATjn3jCzekAicDKQZWZ3A2cCMcBEM8vez3vOuWnBeBEi4p/MzCwe7/kOCTOb8HPdjbT76AI6dmvidyzJocBC75wbUsD6bUDDXFalA22LmUtEwsC2lHTGXfQhCWvimNcyiVtmDiGmblW/Y8lxNPGIiBTL3BlJrP7d93TY1ZRZvTfxyH91PD5U6aciIkU2+uk57Oy/nLr7arDivgxGTL1WRT6EqUcvIkXy5JUf0OmT2mytcYiY/7TmrsGaMCTUqdCLSKHs33eE5y94h26L4ljUZC39pvejWYsYv2NJIajQi0iBVv60jS8vm0a3lDhmdVrL32ZdT6XKKh/hQj8pEcnXpLcXceiPSZx+4FTmXZfKiLdv8juSFJHOnohInp7/42Qq/X4bUZkV2PHv6jz49iC/I0kxqEcvIr+ReTSLx3sFboJaXXcjHT5OIP68xn7HkmJSoReRX9mWks64Cz8kISmOb89cw60zf0fMKboJKpyp0IvIMTlvgprZZxOPfnaTro+PAPoJiggAo5+aQ2q/FdTdV5MV9x/lscm6CSpSqEcvIjxx5Qd0/qQ2W2scpM7Is7hrUCu/I0kJUqEXKcP27zvC8+e/Q7fFcfwYu5Yrpven6Wm1/Y4lJUyFXqSMCtwENTUwE1TntTw4UzdBRSr9VEXKoEljF3HoT0mcfqAB86/fzWNjdRNUJNOZFpEy5t93/pdKt2wjKqs8O5+vzgNjB/gdSYJMPXqRMiIjI5Mne71DwqxYVtXbSPzHCcSfq5ugygIVepEyYNumdMZd9BEJSc35ttUa/jBrKNExVfyOJaVEhV4kws2dnsTPQ7+nw65YZl2awiOf6iaoskY/bZEINvofs0ntv4I6+2qw6oFMRvz3GhX5Mkg9epEI9cSgD+g8MXATVN0323DHwDP9jiQ+UaEXiTD7072ZoLyboAbMuILYuGi/Y4mPVOhFIsjyRVuZefkXdNscx8wua3nwa90EJSr0IhFj4pgfOXJ3MqcdrM/8G3bx2BjdBCUBOisjEgH+fcfnVL51BxWyyrH7+Ro8MGag35EkhKhHLxLGMjIyebLnOyTMjmVVvQ10nHAR7bs08juWhBgVepEwtWVjGv/X/WPdBCUFUqEXCUPfTFtD0rULaL8rltmXpvCwboKSfOiTIRJmRj0xi90DVhKz72RWPZDJcN0EJQVQj14kjDw+cDxdJtVhS80D1HuzLXcM0E1QUrACuwFmNtrMdpjZsjzWn2Fm35nZETO777h1vc1stZklmdkDJRVapKxJTzvM4+1G0W1CfX5qsoHzv7+YXiryUkiF+XtvLNA7n/W7gT8Cz+ZcaGblgVeBS4AzgSFmpk+mSBEtX7SVd858n24/xTGrSzJ3rrhed7pKkRRY6J1zcwgU87zW73DOLQAyjlvVEUhyziU7534B3gf6nUhYkbJmwuiFLE34lrgd9Zk/bDcj5t2oO12lyIJ5BqcBsCnH8xRvmYgUwnO3f85Jt+2kfFY5dr9QgwdGayYoKZ5gdg0sl2Uuz8ZmtwC3ADRurFlvpOzKyMjkyR7vkDAnMBNUxwkX6iYoOSHB7NGnADk/nQ2BLXk1ds6NdM7FO+fi69SpE8RYIqFry8Y0nm85hoQ5sXzbKokrlw1SkZcTFswe/QLgNDNrCmwGrgZ+F8T9iYS12VN+Jvn6BXTY1dS7CepGXR8vJaLAQm9m44EEIMbMUoDhQBSAc+4NM6sHJAInA1lmdjdwpnMu3czuBL4AygOjnXPLg/MyRMLbqMdnUeepdGq7k1n9YBbDn7jG70gSQQos9M65IQWs30bgsExu66YAU4oXTSTyOed4YuAHdJ1Uh5RaBzj1rXbc3r+l37Ekwug6LRGfpKcd5sXz36XbkjgWNl3LwBlXENtc18dLyVOhF/HB8h+3MuvyL+i2JY6ZXdfy4FeaCUqCR58skVI24a2FZPx5HXGH6jP/xj089pZmgpLg0il9kVL07B8+46Q/7MQwdr9YkwfeusLvSFIGqEcvUgpy3gS1ov4GOk+4iPaddX28lA4VepEg27whjfEXfUxCcnPmtl7DnbOvoWb0SX7HkjJEhV4kiGZP+Zl11yVy9u5YZvfdzMMTNBOUlD594kSCZOTfZ7J34GqiD1Rn9UOO4Z8OVZEXX6hHL1LCMjOz+MegD+j66Smk1NpHw9Htub3fGX7HkjJMhV6kBKWnHeYl7yaoxGZruXLGABo3q+V3LCnjVOhFSsjSxC3M6TeD87bEMfPc5MBNUJX0X0z8p0+hSAn4ZFQimfdsoPmhunx/0x4eG3Wj35FEjtGZIZET9Oxtn1H1D6kApL0czV9H6SYoCS3q0YsUU0ZGJk9eNI6EuU1ZWX8DXSZ1p13HXAdyFfGVCr1IMaSs38v47p+QkNycb85aw12zdBOUhC4VepEimvXf1ay/YSHtdROUhAl9OkWKYOSIr0kb/DPRB6qz5hF0E5SEBfXoRQph7+5DvDDkY7rNOJWUWvtpOOZsbuurm6AkPKjQi+Qj45dMnrvlU2InZpGQ3ogfmidx9fSBuglKwooKvUguMjOzeO1vM6g6JpXOqQ1IqpPCiocyuO/vN1GunPkdT6RIVOhFjvPui9+x91+rOWtzLFtqRPH9TXu45/UhREWV9zuaSLGo0It4pn20nCUPzqdjUnOqVKnJ7P5b+PPYAZxco7Lf0UROiAq9lHmJ32xg6p0z6Lq0Ka2j6jPzgnX8fuwVDIit6Xc0kRKhQi9lVvLPqYwb9hmdfziVzllNmNt+HQPe6EOf+FP9jiZSolTopczZnXqQl677hPYza5JwuBnfnbGGTv86n+GX9fQ7mkhQqNBLmXHk8FH+/ftJNP/USNjXiEVNkmnwt3r87dbf+x1NJKhU6CXiZWZm8cpfvqDGuD102XUqP9fdxO57avLnERpKWMoGFXqJaG8/O5f9zyfRdkssm2ums+DWNO5+5XdEVdClklJ2qNBLRJr8/lJWPvQD8cnN2V21JnMGbOWesQOoVr2S39FESp0KvUSU72evZ8adX9JleVNaRtVn5kXruXXsFQxoVMPvaCK+UaGXiJC0Yifv3vQ5XRIb0NE1Zm6H9Qx6sw+XtqvvdzQR3xU4vqqZjTazHWa2LI/1ZmYvmVmSmS0xs/Y51q03s6VmttjMEksyuAhA6vYDjOg9jpXtv+f8+bEsPG0jlT5vyvAFN9FKRV4EKFyPfizwCjAuj/WXAKd5X52A171/s13onEs9gYwiv3HoYAb/vmkiLf5bgYT9jfkxdi1NHm3Hg8N0qaTI8Qos9M65OWYWm0+TfsA455wD5ptZTTOr75zbWkIZRY7JzMzipXumEv1uOufurs+qehvZ+9do7nn4Jr+jiYSskjhG3wDYlON5irdsK+CA6WbmgP8450bmtREzuwW4BaBx48YlEEsizdh/fsOhF5M5e2sTNtXax8I79nH3i9dohieRApREoc9tcG7n/Xuuc26LmZ0CzDCzVc65ObltxPslMBIgPj7e5dZGyqZP31lM0vCFdFjXnNSqJ/PNldu5b8xATqoS5Xc0kbBQEoU+BWiU43lDYAuAcy773x1mNhHoCORa6EWON+/LZGbd/TWdlzfj9Er1mNlzA3e8PYBB9av7HU0krJTE37yfAdd5V990BtKcc1vNrE0XbosAAA9bSURBVKqZVQcws6pALyDXK3dEclq1dDsjOr/F/t5JxK9qzJxO64hbEM9j06/nFBV5kSIrsEdvZuOBBCDGzFKA4UAUgHPuDWAK0AdIAg4Cw7xvrQtMNLPs/bznnJtWwvklguzYuo/XrptA/DcxnH+kKd+1SubCFy5iRI9efkcTCWsWuFgmtMTHx7vERF12X1YcOpjBszdMoOWUisQcqEVis7XEDW9P/+vO9juaSNgws4XOufjc1unOWPFNZmYWL/5xMnXGH6DbnnqsqL+B9Edqc99fdamkSElSoRdfjHpiFpmvbKT99sZsiN7Poj8d4I/PXatLJUWCQIVeStXEMT+y7u+Lab++GTurVWPukB3cO2qQLpUUCSIVeikV33yxhrl/nk2nlc04rdIpzOq9iTvHDmBw3ap+RxOJeCr0ElQrf9rGhzdPpsuiRrS3Rszuso5r3+rL5S3r+B1NpMxQoZeg2JaSzhvXT+ScuXU4LyOW71qvo+fLPXjsgov9jiZS5qjQS4nav+8Izw+bSKuplUk42IQFzZM484lOPHx1d7+jiZRZKvRSIjIzs/j37Z9z6oeH6ba3HstPXc++x2L4y303+x1NpMxToZcTNvKxr+H1zZyzvRHrah/gp3sPcec/r9OlkiIhQoVeiu2jNxPZ/MQS2m1sxo7qVZh3TSr3vnkllSrrYyUSSvQ/Uops1uSf+e7eOXRZHUfFynWY3SeFu94ewJUxVfyOJiK5UKGXQluauIUJt02h66ImnF2uITPPTeb60X3p1yLG72gikg8VeilQyvq9jBo2iY7zTuG8jKZ82yaZPq9czGPn9fY7mogUggq95Ck97TAvXD+B1jOqkHAwlh/ikmjzVBceHaRLJUXCiQq9/EZGRibP3/45DT/6hfPTTmVZg/Uc/Edd7v+TLpUUCUcq9PIrrz/yJVH/2UbHnQ1JjtnM0vuPcPs/dKmkSDhToRcA3n/1e7b/cwVtNzVl68mVmH/9bu5982qiosr7HU1ETpAKfRn31aerWHD/XDr/HMdJJ9Vm9uWb+dPYAdSMPsnvaCJSQlToy5iU9XuZ8cFSNs7bSMXVhzhnTSxtyzdkVrd13Di2P/2a1fI7ooiUMBX6CJa0cidff7iULd9vplLyERruqEaDPafQlHI0pQG7q+5lXtv19H29NyM66VJJkUilQh8hlvywmdmfLGPXwm1U2ZBBox01qJ9ehxaUowWN2F49lY2n7GVtu8NEt6/HeQNbkdClEQP8Di4iQadCH2YyM7P4YfYG5n+2grRFO6i+KYsmO2sRsz+as6gENCGl5nY21ktjdedD1O14KgmDWpPQtp7f0UXEJyr0ISwzM4vZU9awcMoqDizZRc0UiN0ZQ81DJ3M2Vcm0xmyK3s7PjXezvNkhGnZtSPfBZ5EQF+13dBEJISr0IeLI4aN8OWklS6ev4cjyPURvKUfTnadQ7UhVzqEGGeWqsKH2dpbF7cDFHSC2WxMuvvosutev7nd0EQlxKvQ+2L/vCNM/WsbKr5I5uiqNU7ZGEZtal6oZlelMNEcqVGVdzHZ+bLmVcqefTIuEpvS6sjU9dcmjiBSDCn2QpW4/wBcfLCF5zgZYs5+62ysRm1qX6MyKnEsdDlasyrqYHfzQdjMVW9bkzJ7N6TWwFRdXifI7uohECBX6ErRlYxpfjF/CxnkbqbD2EKduP4nGu+rRwJWnAfVIr7yPdXVSmXfOJqq0rk273i3o2fd03X0qIkGlQl9MhblGfX2dXWw4YyPV28bQ8fIzuLRHc40ZIyKlToW+EJYmbmH2x8tITdwauEZ9Zw3qpxV8jbqISChQoc8hMzOLxLkb+W7SCvYu2kG1TZk02VGLOvujaU1FoAmba25nY900VnfSNeoiEh7KbKHPzMzim2lJJP53FfuXpFJjM8TurE2tgzVoRxUyrREptXawptFuVjQ7SIMuDelxVRtdoy4iYafAQm9mo4HLgB3Ouda5rDfgRaAPcBC4wTn3o7eut7euPDDKOfd0CWYvtPyuUY/nZI6Wq8KG2ttY3mwnWXEHie3WmF5Xtqd7w5P9iCsiUqIK06MfC7wCjMtj/SXAad5XJ+B1oJOZlQdeBXoCKcACM/vMObfiREPnp6Br1H8pX411MdtY1HIr1qI6pyU05eKrzqKHrlEXkQhVYKF3zs0xs9h8mvQDxjnnHDDfzGqaWX0gFkhyziUDmNn7XtugFPr0tMO8f/qHxKbWO3aN+qGoaiTX2c4PbVKIalmTM3s0p8eAM+lVvVIwIoiIhKSSOEbfANiU43mKtyy35Z3y2oiZ3QLcAtC4ceMihzi5RmVSo39hS+wmTmodTdveLejR7wxdoy4iZV5JFHrLZZnLZ3munHMjgZEA8fHxebbLz4MrNHm1iMjxSqLQpwA5LxpvCGwBKuaxXERESlFJ3Kb5GXCdBXQG0pxzW4EFwGlm1tTMKgJXe21FRKQUFebyyvFAAhBjZinAcCAKwDn3BjCFwKWVSQQurxzmrTtqZncCXxC4vHK0c255EF6DiIjkozBX3QwpYL0D7shj3RQCvwhERMQnGmFLRCTCqdCLiEQ4FXoRkQinQi8iEuEscC41tJjZTmBDMb89BkgtwTglRbmKRrmKRrmKJhJzNXHO1cltRUgW+hNhZonOuXi/cxxPuYpGuYpGuYqmrOXSoRsRkQinQi8iEuEisdCP9DtAHpSraJSraJSraMpUrog7Ri8iIr8WiT16ERHJQYVeRCTChWWhN7PeZrbazJLM7IFc1puZveStX2Jm7UMk1xlm9p2ZHTGz+0ojUxGyDfXeqyVmNs/M2oZIrn5epsVmlmhm54VCrhztzjGzTDMbFAq5zCzBzNK892uxmT0aCrlyZFtsZsvNbHYo5DKzv+R4r5Z5P8voEMhVw8w+N7OfvPdr2Ant0DkXVl8EhjxeCzQjMLnJT8CZx7XpA0wlMMtVZ+D7EMl1CnAO8CRwX4i9Z12BWt7jS0LoPavG/84ltQFWhUKuHO2+JjBC66BQyEVgSPH/ltZnqwi5ahKYL7qx9/yUUMh1XPvLga9DIRfwIPBP73EdYDdQsbj7DMcefUe8Scedc78A2ZOO53RswnLn3Hwge8JyX3M553Y45xYAGUHOUpxs85xze7yn8wnMCBYKufY779MOVCWf6ShLM5fnLuATYEcpZCpKrtJWmFy/AyY45zZC4P9CiOTKaQgwPkRyOaC6mRmBzs5u4GhxdxiOhT6vyciL2saPXH4parabCPxFFGyFymVmV5jZKmAycGMo5DKzBsAVwBulkKfQuTxdvD/5p5pZqxDJ1QKoZWazzGyhmV0XIrkAMLMqQG8Cv7hDIdcrQEsC068uBf7knMsq7g5LYs7Y0laYSceLNDF5CfFjn4VV6GxmdiGBQl8ax8ILlcs5NxGYaGbnA48DPUIg1wvAX51zmYFOV6koTK4fCYx5st/M+gCTgNNCIFcFoAPQHTgJ+M7M5jvnfvY5V7bLgW+dc7uDmCdbYXJdDCwGLgKaAzPM7BvnXHpxdhiOPfq8JiMvahs/cvmlUNnMrA0wCujnnNsVKrmyOefmAM3NLCYEcsUD75vZemAQ8JqZ9fc7l3Mu3Tm333s8BYgKkfcrBZjmnDvgnEsF5gDBPuFflM/X1ZTOYRsoXK5hBA51OedcErAOOKPYewz2iYcgnMioACQDTfnfiYxWx7W5lF+fjP0hFHLlaDuC0j0ZW5j3rDGBeX+7hliuOP53MrY9sDn7eSj8LL32Yymdk7GFeb/q5Xi/OgIbQ+H9InAY4iuvbRVgGdDa71xeuxoEjoFXDfbPsAjv1+vACO9xXe9zH1PcfYbdoRuXx6TjZnabtz7PCcv9zmVm9YBE4GQgy8zuJnC2vVh/jpVkNuBRoDaBninAURfk0f0KmWsgcJ2ZZQCHgKuc9+n3OVepK2SuQcAfzOwogffr6lB4v5xzK81sGrAEyAJGOeeW+Z3La3oFMN05dyCYeYqY63FgrJktJdBh/asL/CVULBoCQUQkwoXjMXoRESkCFXoRkQinQi8iEuFU6EVEIpwKvYhIhFOhl5BlZjXN7PYcz081s4+DtK/+pTXSY2F4QwXkeXmrmT1rZheVZiYJXyr0EspqAscKvXNui3MuWMMB3w+8FqRtB8PLQJ7DJ4vkpEIvoexpAkMeLDazf5lZrJktAzCzG8xskjdm9zozu9PM7jGzRWY2P3tMcTNrbmbTvIG0vjGz39xGbmYtgCPZN6SY2WBvbPKfzGyOt6y8l2GBBcbHvzXH999vZku99k97y9p5OZaY2UQzq+Utn2Vm/zSzH8zsZzPr5i0/ycze99p/QGA8mOz9jvXyLDWzPwM45zYAtb2b8ETyFXZ3xkqZ8gCB2+TbAZhZ7HHrWwNnA5UJ3AX9V+fc2Wb2PHAdgYHHRgK3OefWmFknAr324w95nEtgMLBsjwIXO+c2m1lNb9lNQJpz7hwzqwR8a2bTCYw/0h/o5Jw7mGPSinHAXc652Wb2d2A4cLe3roJzrqM36NhwAoO0/QE46Jxr4405lJ2nHdDAOdfaew+y8+C1OZfSGXFRwpgKvYSzmc65fcA+M0sDPveWLwXamFk1AhOqfJRjhMlKuWynPrAzx/NvCdx+/iEwwVvWy9tm9qGjGgRGhewBjHHOHQRwzu02sxpATedc9ixKbwMf5dh+9jYXArHe4/OBl7xtLDGzJd7yZKCZmb1MYJjm6Tm2swM4Nbc3RiQnFXoJZ0dyPM7K8TyLwGe7HLA3+y+CfBwiULgBcM7d5vX+LwUWm1k7AuON3OWc+yLnN5pZb4o+HHV2zkx+/X8wtyGa91hgWseLgTuAK/nfmPyVvewi+dIxegll+4Dqxf1mb7C4dWY2GI7NJZzb0LgrCYySideuuXPue+fco0AqgSFlvyAwWFiU16aFmVUl0MO+0QITV2Bm0c65NGBP9vF34FqgoDlS5wBDvW20JjBtIt4Qw+Wcc58AjxAYwTNbCwKjQIrkSz16CVnOuV1m9q13AnYq8GoxNjMUeN3MHgaiCEzb9tNxbeYAz5mZeSM9/svMTiPQi//Ka7+EwGGWHy1wHGgn0N85N83r8Sea2S8ERk59ELgeeMP7BZBMwSOovg6M8Q7ZLAZ+8JY38JZnd8r+BuD9wokjMBqqSL40eqUIYGYvAp875770O0thmNkVQHvn3CN+Z5HQp0M3IgH/IDAhRrioADzndwgJD+rRi4hEOPXoRUQinAq9iEiEU6EXEYlwKvQiIhFOhV5EJML9P+XzJtJnIflZAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "from matplotlib import pyplot as plt\n", "\n", "from numlabs.lab4.lab4_functions import initinter41,rk4ODEinter41,rkckODEinter41\n", "initialVals={'yinitial': 1,'t_beg':0.,'t_end':1.,'dt':0.2,'c1':-1.,'c2':1.,'c3':1.}\n", "coeff = initinter41(initialVals)\n", "\n", "timeVec=np.arange(coeff.t_beg,coeff.t_end,coeff.dt)\n", "nsteps=len(timeVec)\n", "ye=[]\n", "ym=[]\n", "yrk=[]\n", "yrkck=[]\n", "y1=coeff.yinitial\n", "y2=coeff.yinitial\n", "yrk.append(coeff.yinitial)\n", "yrkck.append(coeff.yinitial)\n", "for i in np.arange(1,nsteps):\n", " ynew=rk4ODEinter41(coeff,y1,timeVec[i-1])\n", " yrk.append(ynew)\n", " y1=ynew \n", " ynew=rkckODEinter41(coeff,y2,timeVec[i-1])\n", " yrkck.append(ynew)\n", " y2=ynew \n", "analytic=timeVec + np.exp(-timeVec)\n", "theFig,theAx=plt.subplots(1,1)\n", "l1=theAx.plot(timeVec,analytic,'b-',label='analytic')\n", "theAx.set_xlabel('time (seconds)')\n", "l2=theAx.plot(timeVec,yrkck,'g-',label='rkck')\n", "l3=theAx.plot(timeVec,yrk,'m-',label='rk')\n", "theAx.legend(loc='best')\n", "theAx.set_title('interactive 4.3');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Python: moving from a notebook to a library\n", "\n", "### Managing problem configurations\n", "\n", "So far we've hardcoded our initialVars file into a cell. We need a strategy for saving\n", "this information into a file that we can keep track of using git, and modify for\n", "various runs. In python the fundamental data type is the dictionary. It's very\n", "flexible, but that comes at a cost -- there are other data structures that are better\n", "suited to storing this type of information.\n", "\n", "#### Mutable vs. immutable data types\n", "\n", "Python dictionaries and lists are **mutable**, which means they can be modified after they\n", "are created. Python tuples, on the other hand, are **immutable** -- there is no way of changing\n", "them without creating a copy. Why does this matter? One reason is efficiency and safety, an\n", "immutable object is easier to reason about. Another reason is that immutable objects are **hashable**,\n", "that is, they can be turned into a unique string that can be guaranteed to represent that exact\n", "instance of the datatype. Hashable data structures can be used as dictionary keys, mutable\n", "data structures can't. Here's an illustration -- this cell works:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{(0, 1, 2, 3): 5}\n" ] } ], "source": [ "test_dict=dict()\n", "the_key = (0,1,2,3)\n", "test_dict[the_key]=5\n", "print(test_dict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "this cell fails:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "Traceback (most recent call last):\n", " File \"\", line 5, in \n", " test_dict[the_key]=5\n", "TypeError: unhashable type: 'list'\n" ] } ], "source": [ "import traceback, sys\n", "test_dict=dict()\n", "the_key = [0,1,2,3]\n", "try:\n", " test_dict[the_key]=5\n", "except TypeError as e:\n", " tb = sys.exc_info()\n", " traceback.print_exception(*tb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Named tuples\n", "\n", "One particular tuple flavor that bridges the gap between tuples and dictionaries\n", "is the [namedtuple](https://docs.python.org/3/library/collections.html#collections.namedtuple).\n", "It has the ability to look up values by attribute instead of numerical index (unlike\n", "a tuple), but it's immutable and so can be used as a dictionary key. The cell\n", "below show how to convert from a dictionary to a namedtuple for our case:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "values are -1.0 and 1\n" ] } ], "source": [ "from collections import namedtuple\n", "initialDict={'yinitial': 1,'t_beg':0.,'t_end':1.,\n", " 'dt':0.2,'c1':-1.,'c2':1.,'c3':1.} \n", "inittup=namedtuple('inittup','dt c1 c2 c3 t_beg t_end yinitial')\n", "initialCond=inittup(**initialDict)\n", "print(f\"values are {initialCond.c1} and {initialCond.yinitial}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comment on the cell above:\n", "\n", "1) `inittup=namedtuple('inittup','dt c1 c2 c3 t_beg t_end yinitial')`\n", " creats a new data type with a type name (inittup) and properties\n", " (the attributes we wlll need like dt, c1 etc.)\n", " \n", "2) `initialCond=inittup(**initialDict)`\n", " uses \"keyword expansion\" via the \"doublesplat\" operator `**` to expand\n", " the initialDict into a set of key=value pairs for the inittup constructor\n", " which makes an instance of our new datatype called initialCond\n", " \n", "3) we access these readonly members of the instance using attributes like this:\n", " `newc1 = initialCond.c1`\n", "\n", " \n", "Note the other big benefit for namedtuples -- \"initialCond.c1\" is self-documenting,\n", "you don't have to explain that the tuple value initialCond[3] holds c1,\n", "and you never have to worry about changes to the order of the tuple changing the \n", "results of your code." ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0 }, "source": [ "### Saving named tuples to a file\n", "\n", "One drawback to namedtuples is that there's no one annointed way to **serialize** them\n", "i.e. we are in charge of trying to figure out how to write our namedtuple out\n", "to a file for future use. Contrast this with lists, strings, and scalar numbers and\n", "dictionaries, which all have a builtin **json** representation in text form.\n", "\n", "So here's how to turn our named tuple back into a dictionary:\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "OrderedDict([('dt', 0.2), ('c1', -1.0), ('c2', 1.0), ('c3', 1.0), ('t_beg', 0.0), ('t_end', 1.0), ('yinitial', 1)])\n" ] } ], "source": [ "#\n", "# make the named tuple a dictionary\n", "#\n", "initialDict = initialCond._asdict()\n", "print(initialDict)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why does `_asdict` start with an underscore? It's to keep the fundamental\n", "methods and attributes of the namedtuple class separate from the attributes\n", "we added when we created the new `inittup` class. For more information, see\n", "the [collections docs](https://docs.python.org/3.7/library/collections.html#module-collections)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "outputDict = dict(initialconds = initialDict)\n", "import json\n", "outputDict['history'] = 'written Jan. 28, 2020'\n", "outputDict['plot_title'] = 'simple damped oscillator run 1'\n", "with open('run1.json', 'w') as jsonout:\n", " json.dump(outputDict,jsonout,indent=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After running this cell, you should see the following [json output](https://en.wikipedia.org/wiki/JSON) in the file `run1.json`:\n", "\n", "```\n", "{\n", " \"initialconds\": {\n", " \"dt\": 0.2,\n", " \"c1\": -1.0,\n", " \"c2\": 1.0,\n", " \"c3\": 1.0,\n", " \"t_beg\": 0.0,\n", " \"t_end\": 1.0,\n", " \"yinitial\": 1\n", " },\n", " \"history\": \"written Jan. 28, 2020\",\n", " \"plot_title\": \"simple damped oscillator run 1\"\n", "}\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Reading a json file back into python\n", "\n", "To recover your conditions read the file back in as a dictionary:\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "values are -1.0 and 1\n" ] } ], "source": [ "with open(\"run1.json\",'r') as jsonin:\n", " inputDict = json.load(jsonin)\n", "initial_conds = inittup(**inputDict['initialconds'])\n", "print(f\"values are {initial_conds.c1} and {initial_conds.yinitial}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Passing a derivative function to an integrator\n", "\n", "In python, functions are first class objects, which means you can pass them around like any\n", "other datatype, no need to get function handles as in matlab or Fortran. The integrators\n", "in [test.py](https://github.com/phaustin/numeric/blob/master/numlabs/lab4/example/test.py)\n", "have been written to accept a derivative function of the form:\n", "\n", "```python\n", " def derivs4(coeff, y):\n", "```\n", "\n", "i.e. as long as the derivative can be written in terms of coefficients\n", "and the previous value of y, the integrator will move the ode ahead one\n", "timestep. If we wanted coefficients that were a function of time, we would\n", "need to also include those functions the coeff namedtuple, and add keep track of the\n", "timestep through the integration.\n", "\n", "Here's an example using foward euler to integrate the harmonic oscillator\n", "\n", "Note that you can also run this from the terminal by doing:\n", "\n", "```\n", "cd numlabs/lab4/example\n", "python do_example.py\n", "```" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfAAAAHwCAYAAABZrD3mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOzde3zcVZ0//tdJMrlfe6UdKChIkMLC2AIV3FXwEpTLzhcUdMH1tou6rrvszw1QQFoUhbXrij/3uyreXVi5M4CKVUFdQYu0DNAWWkAuLZ/Se9LcJ5P0fP/4zEkmk/lMZjKfz+ecM309Hw8eD5JMk08mybw/533e5/0WUkoQERGRXap0XwARERGVjgGciIjIQgzgREREFmIAJyIishADOBERkYUYwImIiCzEAE4VQwhxtRDiuwF97t8KIf7Op8/1ihDiXX58riAJIaQQ4piQvta3hBCfz/z/O4QQr2V9zIrniyhsNbovgMgvUsov674Gmh0p5af8+DxCCAngTVLKF/34fGEQQnwRQBzAmwHcIKVcrfeKyBZcgRMRARBCzLigKeYxs/AigCsA/CyAz00VjAGcrCOEuFII4Qgh+oUQW4UQ78y8f7UQ4tbM/x+VSQF/TAixXQjRI4T4lBDiFCHEM0KIXiHEf2Z9zo8KIR4TQnxDCHFACLFFfV6Pa/i4EOK5zOddK4Q4ssBjPyyEeFUIsU8IcU3Ox04VQvwxcz2vCyH+UwhRm/VxKYT4ByHEC5nv94tCiKMz/6ZPCHGnerxKPWe2EvZmUs+XZH2uOiHEvwshtgkhdmXS1g1ZH+/OXMMOIcTHZ/gZLBZCPCCE2C+EeFEI8fc539P6zPXtEkL8R9bH3iaE+EPm+90uhPho5v0/FELcUOhrzvR8CSH+N/Owp4UQA0KIizPv//vMNe7PXPPinOf3M0KIFwC8kOfrqd+jTwghtgF4JDfFn3ncRJo/83t4pxDix5mf2WYhxHKv70lK+SMp5UMA+mf6/omyMYCTVYQQnQD+EcApUsoWAF0AXinwT04D8CYAFwO4GcA1AN4FYCmAi4QQb8957EsA5gFYBeBeIcScPNcQB3A1gAsAzAfwewA/8bje4wF8E8CHASwGMBfA4VkPGQfwL5mv+VYA7wTwDzmf5mwAywCsgLtSuwXAJQCOAHACgA9lPfawzOeKAvgIgFsyzxkA/BuAYwGcDOCYzGOuy1zn2QD+FcC7M8/XTHvOPwHwWuZ7ej+AL2fd8HwdwNellK0AjgZwZ+ZrLAHwEIBvwH3eTgbw1AxfJ5fn8yWl/KvMY06SUjZLKe8QQpwF4EYAFwFYBOBVALfnfM443J/98QW+7tvhpri7irzO8zNfpx3AAwD+s/DDiUrHAE62GQdQB+B4IURESvmKlPLPBR7/RSnliJTylwAGAfxESrlbSunADbyxrMfuBnCzlDItpbwDwFYA5+T5nJ8EcKOU8jkp5RiALwM42WMV/n4AP5VS/q+UMgXg8wAOqg9KKTdIKddJKceklK8A+DbcYJHt36SUfVLKzQA2AfillPIlKeUBuAExlvP4z0spU1LK38FNy14khBAA/h7Av0gp90sp+zPX/cHMv7kIwA+klJuklIMAVud9NgEIIY4A8DYAV2ae26cAfBfuTQoApAEcI4SYJ6UckFKuy7z/EgC/llL+JPMc78v826IV+XxluwTA96WUT2ae/5UA3iqEOCrrMTdmnpPhAp9ntZRycIbHZHtUSvlzKeU4gP8GcFKR/46oaAzgZJVMcdLlcAPMbiHE7dkp0Tx2Zf3/cJ63m7PeduTU6T6vwl1h5joSwNczadxeAPsBCLgr2lyLAWzPuv5BAPvU20KIY4UQPxVC7BRC9MENqvPK+B56Ml8j93uYD6ARwIas6/5F5v3TrjPz77wsBqBuArIfr77/T8Bd6W8RQjwhhDg38/4jABS62ZpRkc9X7rVOfC9SygG4z3/2z2p77j/Ko5jHZNuZ9f9DAOpFMPvndAhjACfrSCn/R0r5NriBVMJNDfshmlmpKksA7MjzuO0APimlbM/6r0FK+Yc8j30dbuACAAghGuGm0ZVvAtgCt3K6FW5qXmD2OoQQTXm+h71wg/3SrGtuk1Kq4D/lOjP/zssOAHOEEC05j3cAQEr5gpTyQwAWwP3Z3J25pu1wU+rlKPX52gH39wQAkLmOuepaM4oZyZj9mEG4N0Pqc1Zj8kaIKDQM4GQVIUSnEOIsIUQdgBG4QWncp0+/AMA/CSEiQogPwN3z/Hmex30LwEohxNLMNbVlHp/P3QDOzRRv1QL4Aqb+3bUA6AMwIIQ4DsCnffg+rhdC1Aoh/hLAuQDuklIeBPAdAF8TQizIXHdUCKH2dO8E8FEhxPGZm4xVXp9cSrkdwB8A3CiEqBdC/AXcVfdtmc97qRBifuZr9mb+2Xjm4+8SQlwkhKgRQswVQpxc4vc20/O1C8Abs97+HwAfE0KcnPmd+TKAxzPp99l6Hu6K+hwhRATAtXC3dWYl8/tWD/f3oibznFaXcX10iGAAJ9vUAbgJ7opyJ9yge7VPn/txuAVcewF8CcD7pZT7ch8kpbwP7sry9kwadxOA9+b7hJl968/ADSSvA+iBW/yl/CuAv4FbgfwdAHeU+T3szHyNHXAD5qeklFsyH7sS7pGldZnr/jWAzsx1PgS3yO+RzGMemeHrfAjAUZmvcx+AVVLKX2U+djaAzUKIAbgFbR/M7JVvA/A+AJ+Du+3wFErfG57p+VoN4EeZbYKLpJQPw607uAfu8380Jvf9ZyVTe/APcPf9Hbgr8tcK/qPCvgP3RvRDcIsshzFZT0DkSUzd8iM6NGWOM/1dJjVvJSHEOwDcKqU8fKbHEpH9uAInIiKyEAM4ERGRhZhCJyIishBX4ERERBZiACciIrKQVZ2B5s2bJ4866ijdl0FERBSKDRs27JVS5m0UZFUAP+qoo7B+/Xrdl0FERBQKIYRnW2Om0ImIiCzEAE5ERGQhBnAiIiILMYATERFZiAGciIjIQgzgREREFmIAJyIishADOBERkYUYwImIiCzEAE5ERGQhBnAiIiILMYATERFZiAGciIjIQgzgREREFmIAJyIispBV88CJiIh0SiQdrH5gM3qH0wCAjsYIVp23FPFYNPRrYQAnIiIqQiLpoPuup5E+KCfe1zOURvfdTwNA6EGcKXQiIqIiXP/g5inBW0mPS6x+YHPo18MATkRENINE0kHPUNrz473DaSSSTohXxABOREQ0ozVrt/ryGD8xgBMREc1gR++wL4/xEwM4ERHRDNobIzM+pq1h5sf4iQGciIiogETSwcDI2IyPGxwdC3UfnAGciIiogDVrt+atPs+VHpeh7oMzgBMRERVQyt52mPvgDOBEREQFLG5vCOSx5WIAJyIiKqC7qxO1NVPDZaRKIFItpryvIVKN7q7O0K6LAZyIiGgG2bG6ozGCNR84CWvefxJa6qon3l8fCTekMoATERF5SCQdrLx3I4bTByfeN5L1/7l90VfeuzG0SnQOMyEiIvKwZu1WDKfHp7xvOD0+UW2eHcyzPxbGYBMGcCIiIg9eVeWFqs3DqkRnCp2IiMiDV1X54vaGgh8LAwM4ERGRh+6uTtTlVKCravPurk40RKrzfiwMTKETERF5iMei2LqzH9/83Z8h4K6uu7s6p+xxX//gZvQMpbGgpQ5Xv+/Noex/AwzgREREBR05txEA8Psrz8ThHY1TPhaPRbGorR4X37IO//6Bk/BXx84P7bqYQiciIvKQSDr40s+eAwBc9O0/5j0iFu1w97wdjhMlIiLST50B70+5k8h29I7kPef9p5f2AQBW3rsRZ9z0SGjnwBnAiYiI8pjpDDjgBvlrEpsn3nZ6h0Nr5sIATkRElEcxZ8CLCfJBYQAnIiLKo5hz3rNp9OIXBnAiIqI8urs6Ue9xBlzR2cyFAZyIiCiPeCyKz5x1zMTb0fYG3HjBiVPOeets5sJz4ERERB6OX9QKAEh85gycfET7tI+rYK6jmQsDOBERkYddfSkAwMLWOs/HxGNRzG+pwyXffRxf/2AMbz16bijXxhQ6ERGRh119IxACmNfsHcABYEGL+/Hd/SNhXBYABnAiIqK8EkkH3/n9S5ASeMea3xY8272gpR4AsKc/FdblMYATERHlUl3YhkbdM94zNWh5ZMsuAMANP3sutG5sWgO4EKJdCHG3EGKLEOI5IcRbdV4PERERUFqDlkTSwdX3bZp4O6xubLpX4F8H8Asp5XEATgLwnObrISIiKqlBi65ubNoCuBCiFcBfAfgeAEgpR6WUvbquh4iISCmlQYuubmw6V+BvBLAHwA+EEEkhxHeFEE0ar4eIiAhApgtbpHAXNkVXNzadAbwGwFsAfFNKGQMwCOCq3AcJIS4TQqwXQqzfs2dP2NdIRESHoHgsin9+55sm3s7XhU3R1Y1NZyOX1wC8JqV8PPP23cgTwKWUtwC4BQCWL18uw7s8IiI6lJ10uNt57Sd/v6JgcxYV1Nes3YodvcNY3N6A7q7OwLuxaQvgUsqdQojtQohOKeVWAO8E8Kyu6yEiIsq2d3AUADCvuXbGx8Zj0VDap2bT3Ur1swBuE0LUAngJwMc0Xw8REREAYN+A25Rl7gxd2HTRGsCllE8BWK7zGoiIiPLZO5BCdZVAe0NE96XkpfscOBERkZH2DYxiTlMtqqqE7kvJiwGciIgoRyLpIJF0sKc/FVpr1FIxgBMREWVRfdBHxg4CCK81aqkYwImIiLLoao1aKgZwIiKiLLpao5aKAZyIiCiLrtaopWIAJyIiytLd1Yn6muL6oOvEAE5ERJQlHovis0X2QddJdyc2IiIi48SWFNcHXSeuwImIiHL0DKYBAB1NZnZhAxjAiYiIptk/5A4ymdM48yATXRjAiYiIcvRmJpG1M4ATERHZY//QKJrralBbY26YNPfKiIiINOkZHDV6/xtgACciIpoikXTw0Kad2L5/2NhBJgADOBER0QQ1yCRl+CATgOfAiYioBImkg9UPbEbvcOaYVWMEq85balyTk9kqNMjEtO+RAZyIiIpybWIjbl23bcr7eobS6L77aQAwLsDNhi2DTACm0ImIqAiJpIPbcoK3kh6XWP3A5pCvKBi2DDIBGMCJiKgIa9ZuhSzw8d7htJH7xKXq7upEfcT8QSYAAzgRERXBKSKFvGbt1hCuJFjxWBSfe/dksDZ1kAnAPXAiIppBIulAAAVX4EBxQd4Gp7xhDgDgex9Zjne+eaHmq/HGFTgRERU0U/pcEUBFpNEPZCrs2xrYyIWIiCxWbAW2RGWk0RnAiYioIpRSgW3icatSMYATEVFFOPO4+dPeJzwea+Jxq1L1ZQJ4KwM4ERHZKpF0cM+GqfvaAsDpR89BQ6R62vvzBXvbHBhOo66mCvU5359pGMCJiMhTvtaiEsAr+4Zx4bLotPffs8GxvpDtwFDa+PQ5wABOREQFFGot+pste6a9X/UNt9mBYQZwIiKyXKHWojb1DS9WIungka278cLuAaNHiQIM4EREVEB3Vyci1VNL1lRrUZv6hhdDjRIdtWCUKMAATkREM8gOFB2NkYnWot1dndMK2UztG16MQqNETcQATkREeakVaWp8sg/bSPrgxP/HY1HceMGJaKmbDOK5g0BsYtuWgL3PNBERBarYFWn64GSA7xlKG512LsS2LQEGcCIiyquYFematVunrMoBs9POhdg0ShRgACciIg/FrEhtSzsXYtMoUYABnIiIPHR3daLWowJdsS3tPJNTM6NEv/u3y/HYVWcZG7wBBnAiIvIQj0XxnqWHAXDbpOZbkVZaJfrEIJNG8xu51Oi+ACIiMlMi6eCRLbsBuCvq7q7OaStS9fYXfvos9g+OYn5zHa45581Gr1wLUQG8nZ3YiIjIRuoI2dCoW4VeqKlJPBbF3//lGwAAewZSWLN2q5VV6IA9o0QBBnAiIsqjlKYmiaSDrz/8wsTbpncwK+SAJaNEAQZwIiLKo5Tq8ko6StZnyShRgAGciIjyKKW6vJKOktkyiQxgACciojzcI2TFNTWppKNkDOBERGS1eCyKC5a5leReR8iUSjpKZlMA5zEyIiLKa2FLPYQAtn7xvait8V7vqaC+6oFNODA8hkVt9bjy7OOsO0qWSDp44pX9SI9LnHHTI3mPzZmEAZyIiPJyeoexoKWuYPBW4rEo6mqq8OnbnsR3P7IcSxe3hXCF/lHH5tKZyWuqkh6AsUGcKXQiIpomkXTwwNM7sKsvhTNueqSoI2GLMnver/eOBH15vrNtFjjAAE5ERDnUanR0zD0aVuy57sVt9QCA1w/YV31uYyU9AzgREU0x29Xooy/sBQB8/v7NRa/aTWFjJT0DOBERTTGb1Wgi6eCaxKaJt23rxtbd1Ym6GntmgQMGBHAhRLUQIimE+KnuayEiotmtRm3cQ84Wj0Vx+bveNPG26bPAATOq0P8ZwHMAWnVfCBERuavRK+5+BqPjk+1RZ1qN2riHnOuUo9xZ4D/++Kn4q2Pna76amWldgQshDgdwDoDv6rwOIiKaFI9F8dexxQBmbuKi2LiHnKt/ZAwA0FJvwtp2Zrqv8mYAVwBo0XwdRESUZVFrPaoE8PwN70VN9cxrve6uTqy8d+OUNLrpe8i5+kbcSWQt9XZ0YtO2AhdCnAtgt5RywwyPu0wIsV4IsX7Pnj0hXR0R0aFtZ98I5jXXFRW8AXfVfuMFJ6K90Q1+C1vrjN9DzqVW4K2WrMB1ptDPAHC+EOIVALcDOEsIcWvug6SUt0gpl0spl8+fb/6eBBGR7RJJB/c/tQO7+4tv4gK4Qfz//2AMAPCND73FquANZKfQuQIvSEq5Ukp5uJTyKAAfBPCIlPJSXddDRESTTVxSJTZxURa01gEAdvfb142tfySN6iqB+oj2A1pFsSNPQERkiETSweoHNqN32N0v7WiMYNV5S61bbXopdBysmO9xQYvbjW13XyqQ6wtS/8gYWuprIITQfSlFMeI2Q0r5Wynlubqvg4iokGsTG3H5HU9NBG8A6BlK4/I7nsK1iY0ar8w/5R4H+93W3QCAL/z0Weu6sfWPpK2pQAcMCeBERKZLJB3ctm6b58dvW7fNqmDlpZzjYImkg6vvs7cb20BqDC11dux/AwzgRERFuf7BzZAFPi4Ba7qOFVJOS1Hbu7H1ZVLotmAAJyKaQSLpoGcoPePjHIu6jnmJx6L48IojARTfxEWxuRtbIungyVd78PjL+61J/dtzq0FEpEmxK0gBNxDYXtB29IJmAMBjV51VUie1xe0NeW9iTO/Gpirvxw66ORaV+gdg9M+SK3AiohkUu7KulDT63n63gnxuc21J/667qxMNkeop77OhG5utqX+uwImICkgkHQig4P53tkpIo+8ZSKGtIYK6muqZH5xFrVZXPbAJB4bHsKitHleefZzRq1jA3tQ/V+BERAWsWbu16OANTKbRbbZ3IIV5Ja6+lXgsitXnLwUA3PZ3pxkfvAF7B7EwgBMRFVBoFZav3YftafRE0sGvn92NP+8ZnHUx19wmtxvbvsFRvy8vEN1dnaifZeW9TgzgREQFeK3Cou0Nnitz01OvXlQxl5oDPttz3HOa3NX7vgE7Ang8FsXl73rTxNulVN7rxABORFTAmcdNH6KkVmdRS1OvXvwq5lLFb/stWYEDwClvmAsA+MHHTsFjV51lfPAGGMCJiDwlkg7u2TB19SkAXLgsingsmrfqWiB/0LeBX8Vckytwe/qhD6TsGiUKMIATEXnKtyKVAH6zZQ8AN/V64bLotI/fs8GxspDNr2KuhzbuhADw1V89b01TlP4Rt1GPLaNEAQZwIiJPxaxIVTDPZsMZ4nzKaaOqqH10VR9gSz/0yVngXIETEVmvmBWprWeI84nHorj0tCUASm+jqtjaFMXGFbg9txpERCE787j5uDVnAlnuitTW9qFejl7QAgD4w8qzsKit9O/B1hua/pExVAmgqba05jU6cQVORJTHTAVsSqUVsu3NFJ6ps9ylsrUpSv/IGJrraiBEvtP9ZmIAJyLKY6YCNqXSCtn29KfQ3hhBbc3swoOt/dD7RtJWpc8BBnAiorxKSQVXUiGb20Z1dqtvwL2hufGCE9HW4O7QLmqrt6IpSr9ls8ABBnAiorxKSQXbuu+bK5F08PBzu/Hi7oGyjn/FY1Fcf/4JAOzph94/kmYAJyKqBN1dnaipmrof6pUKtnXfN5tfbVSViWYulnRjc1fgTKETEVUEkdXtvKMx4pkKtnXfN5vfx79s6oeeSDrYsrMfj2zZbU3jGYDHyIiIplGr0fTByfeNZL+RQwX11Q9uRu9QGgta6nD1+95sRepY8XsbQO2jm94PXf2sxw+6N2sq8wDA+J8fV+BERDlmsxqNx6L46OlHAQB296ewZu1Wa1ZygP/bAB1Nbjp6/6DZ/dBtbTwDMIATEU0zm9VoIungW7/788TbtrQQVfxoo5pN9UP/91+a3Q/d5gJEBnAiohyzWY2uWbt1WprdlpUc4GYQPnbGUQBm30ZVsakfus0FiAzgREQ5urs6EakurgJdsXklpxx3WCsA4OHPvb2smdg2paX9zjyEiQGciChHPBbFijfMgUDxq1GbV3KKKjhTFeSzZdPNTDwWxWffeczE2+VkHsLGAE5ElCORdPCnV3og4Qbg7q7OGV/QK+Eo2f7BUVRXCbSWeR7atpuZU46cA8BtOlNO5iFsDOBERFnU/m1qrLSGJqqFaEejG/wWttZZs5JT9g+NoqMxgqqq8gZ62HYzM5ByZ4E319l1stquqyUirRJJB6sf2Ize4fTE+zoaI1h13lKrAlUhhfZvZ/oe47Eooh0N+MC3/oivvP8kvP1YuyaS7R8YRUdjeelzYPL89HX3b0LfyBgWt9XjirOPM/Z3ZCKAs5UqEVWiRNJB911PTwneANAzlMbldzyFaxMbNV2Zv8rdv312xwEAwEe+/yejj0/ls39otOz9byUei+LK9x4HAEh85gxjgzdg7wqcAZyIinL9g5uRPig9P37rum0VEcTL2b9NJB3c9NCWibdNPj6VK5F08OSrPXj85f2+3Xi0N7g3A7k3faYZZAAnokqVSDroGZr5Rfi2ddusCFaFzOYImeKm3+07C672/cdy2omW+7Nsz9QD9BjeTnVgZAxCAI211TM/2CAM4EQ0o+sf3FzU4yRgfLCaSTwWxTsye9elNjSx6fhUtqDObbc1uAHc9BV4f2oMzbU1EKK84r2w2ZUvIKLQFbv6VhzDg1UxGutqcMScBvz+irNK+neL2xvyfv+mHp9Sgrrx6Mjspx8o4fdHh8HUmHUFbABX4EQ0g1JXYQKwOo2eSDp4aONObN8/XPJesG3Hp5Sgzm23Z1bgPUOGp9BTY2iybP8bYAAnohmUuqK2OY2u9oJHx0s7A66os+Dtlp0F7+7qRG0A7UQba6sRqRbGp9AHUuPWFbABDOBEVEAi6cBrV1CtrvKxNY3ux15wPBbFf37oLQCAr38wZnzwBtxrvvS0JQDKH2SS7f6ndmD8oMQ3f/tno4/UDYykGcCJqLKsWbsV+Q6OCQCrz1+KqEeK1dY0ul97wZtfd8+Cf/CWdUYHrmzHLmwBADx21Vm+tBNV2Qx18tDUI3WJpINnXjuAR1/ca83PSmEAJyJPXoFLwl21dXd15l2h25pG92MvOJF08LVfPT/xtqmBK9c+nwaZKDZMJAvq+FxYGMCJyJPay82lVt7xWDTvCh0w/+hUPuWcAVdsnQveMziKxtpq1Ef8OQttw5E6G24yCmEAJ6K8EkkHAyNj094fqRZTAppXGr2twB65qeKx6ET/8tnuBdsQuPLZP+hPH3TFholktv6sFAZwIsprzdqteVunNtXWTAlo3V2diOSZXjU4OmZNKjJbY20NlsxpxMs3nTOrvWAbAlc+fvZBB+w4Umfrz0phACeivLxWIQdyjgTFY9G8TTDS49KaVKSSSDp4aNPr2LZ/aNYFTTYErnx6Bv0N4OpIXWvmd2NxW71xR+q6uzpRH8DxubAwgBNRXl773/lWJ70enbZsSUUCkwVN6fHyCppyA9ciAwNXrkTSwaYdffjd83t8rcTOnkh2n4ETyeKxKD73nqnbQab/rLLZd/CNiAJX7P63YmsL0WzlzAHPpR5/+R1P4b8/cRqOWdDs23X6Td24jOdUYgPwJZBNTCQbSmNha33Zn89vK944FwDwnb9djncfv1Dz1ZSGK3AimqbY/W8lX9pYADjzuPlBXaLv/C5oWtTmBqvXD5idhQi6EltlcnoNbadq6yxwgAGciPIodv9biceiuHDZ1MAuAdyzwbGmkM3vgib1717vHZn1NYUh6ErsiZGihg40YQAnoooym2D2my17pr3PpjO1fpwBz/b4S/sAAFfc84zRHb6CrsRuzxxNOzBs5gp8UAVwTiMjokrQ3dWJ6qrSgpntZ2rjsSi6lrp7oOX2A08kHXz+/skZ6iZ3+ApqkImieuZ7FTrq1m/xCty+KyaicMjJPfCOxghWnbe0YDCrhEK2OU11aGuI4OlV7ynr8/hZEBe0eCyKp7f34Ad/eBUC7s+ru6vTt+tUE8mMTaGPMICXTAhxBIAfAzgMwEEAt0gpv67reoiKkUg6WP3A5injEYsJbjZxq5KfwXhWDVtua9B8urs6sfLejVMCl01nagHg9QMjE8Vn5bAtG/Gmha0AgD+sPAuL2vy94RJCoL2x1ugUenWVQH3EvoS0ziseA/A5KeWbAawA8BkhxPEar4eooETSQfddT0+bbdwzlMbldzyFaxMbNV2Zv9zVY+m9vNX558asF0KbXhQTSQe/3bobW3b2l71nbVuHr95McFVHvvyUSDroGRzFT/603chagIHUGJpqqyGE1+Bcc2n765JSvi6lfDLz//0AngNQGUsYqkjXP7g579Eq5dZ12yoiiJe7esxeufcMpY3d+83mVxMXxbZubL1DadTVVKGh1p9BJooN074GUmNoqbevbz9gSBGbEOIoADEAj+u9EqL8rk1sLGoP77Z124x6cZqNclaPa9ZuRWrMvklcfp+FVtkIVcC1sLXO6A5fvUP+DjJRbJj2NTAyZuX+N2BAABdCNAO4B8DlUsq+PB+/TAixXgixfs+e6cdUiIKWSDq4bQ+J+Q0AACAASURBVN22oh4rAax+YPOMjzNZvuYrxa4ebdv7VYK47ngsiq9edBIA4NsfXm5s8AbcTIlX69xy2PD7MJAaQ1Odv5mHsGgN4EKICNzgfZuU8t58j5FS3iKlXC6lXD5/vj1dnahyXP/gZs+Z1/n0DqetXYUnkg7u2TD12gWAC5dFiwpAtu39KkFd94IWtyBud5/ZzVx6h0YDCeA2/D4MpMbQzBR6aYRbMfA9AM9JKf9D13UQFZJIOrM6/mLrKjxfylMif5OWfGzb+1X8buKiLGitAwDs7k+V9XmC1juUDiSFbsPvw0BqDM1cgZfsDAAfBnCWEOKpzH/v03g9RNPMdq/O1lV4uSlPWyvR47Eo3rv0MADlN3HJ9tgL7o3PtYlNRlZgA+5N6p/3DOChTTt9v0b1+9DWYOZktkTSwUt7BvDzjf5/72HQtnMvpXwU7t8KkbHyNSZRLl2xBD975nXPFbqJTTtm4lczlnyV6IA/062CMqe5Ds11Ndh0fZcvny+RdHBNYno3NsCc50Gd+VeHK4K4xngsivGDEp+762ncftkKHDm3yZfPWy5VIR/k9x4082+NiTRJJB3PO8z2hghuiJ+IVect9fz3hYK/qbq7OlFTYgvVXLZWou88MILDfGjiothQgT3bM/+las1U4/cNTx9Rq4sNP5+ZMIATeVizdmve4jUBYPX5buCOx6Lo8Cj+EYB1KTkA01qolprytKHyOJ/X+/zpwqbY8DyEdY2tmUEhXtPsdLDh5zMTBnAiD14raImpKbZV5y3Nu1KXmP0eug4TTTdKbKGay4bK41yJpIONr/Xi9y/s9W0v1IbnIaxrbMvc5PaNmBPAbfj5zIQBnCiPQunzaM4feDwW9TxmZtPdvF8pRRsqj7Mlkg6uyrMPXG4Qt+F56O7qRG11cJPIlNZ6lUI3J4B3d3WiLsApbGFgACfKo1D6PN8feG5QV9oa7Dlf6ldKcbLy2P3eD2s1q/I415q1W6dlGvzYC7XheYjHorj41CMA+Ft9n0vtgZuUQo/HovjMmUdPvB3U9x4kO/vHEQXMK2jlps+V7q5OdN/19LRe6YOjY0gkHSteFPwcBxqPRbF1Zz+++bs/Y2ffyEQwNPF5CHIvNB6Loj5ShU/d+iS+99HlWLq4rezP6bdjFzQDAB6/5p0TjWf81lRbjeoqYVQKHQCWHTkHAHDHZStw2hvnar6a0nEFTpSHV1cqr5V2PBZFc/30++H0uLRmH9yPCnQlkXTwg8dennjbxCEWStB7ofNVNzZDm7moY5BBTCJThBBora8xqgodcJu4AMj7t2sDBnAqWSLp4OTrf4mjrvoZjrrqZ4h94ZdGvjDPViLpYGBk+gtNpFoUDGa9HufBbdkHj8eiOCHaiuoqUXY6dc3arRix5ChZUF3YlAUtbje2PX1mBvDeoTSa62pQWxNcOEgkHfSNjOG/171qVMMU9XfOYSZU8RJJB2/+/EO4/I6npszErsR52PnGhjbV1hQMZl4rNlv2wRNJB5ucPowflFjc3oDurs5Zp7xtOqITj0VxZqc7ZyGIfeB1L+0DAFxxzzNGBS+ld2g00N9Rdbph3MCRooOjDOB0CEgkHXTf9fS0pg/ZKn0e9kwFON1dnYhUTa9dV/vgJlMdufya22zbEZ32xlosaKnDyzedg8euOsu34J1IOrju/und2Ez6fegZGkVHU3AB3OSGKf2ZFXgTAzhVsusf3Jx3VZqrEoK41/73TMHH5n1wvzty2XCEKtvu/hQWtvpfwGVy8FJ6h4MZZKKYnI0ZTI0hUi2mHSezhZ1XTaG6NrGxpIlct63bZtQKoxSz3f9WbN0H9/tF1oYjVNl29aUm9qr9ZHLwUnqH0oGm0E3OxriTyGrgDse0DwM4FZRIOrht3baS/o2E3eM0Z7P/rZj8YlVIENcdj0XxocwZ412Zo2Qm3tglkg627uzDw1t2+75HbfrvQyLp4NV9g/jpM68Htj9vcjZmYGTM2vQ5wABOBSSSDj5359OeXcYKqbRxmsU2oMj3YiUAnHnc/HIvLVDdXZ2o9ukImeIeJXsFgHtTZ+L+b1Bd2BSTg5fXJDK/fz4qG6P6oS82aKSoWoHbigGc8pqoHJWzCd8uG1fh5a6Y4rEoLlw29YVJArhng2NU4MoVj0Vx7IJm1PhwhEyxYSpZUF3YFJODV1iTyAD3eeg++zgAQOIfzzDi+wcYwKlC5Su+ySYA3Hzxybh0xRLPx9i2Ck8kHQympu9/l7pi+s2WPdPeZ1rgypVIOnh+9wDGfDhCptiw/xvGNcZjUVxzzpsBAHd+6q3GBK+wfz7qJsakZi6DqTFrm7gADODkodAsawHgkhVLEI9FcUP8RM9xmoA907hUxqE3J1V+KIzTVGlkv8/pmr7/C4R3jfOa3QK5fQOjvn7ecoT98zGxH3o/V+BUaQpN4qoWAl+7+GTcED9x4n2rzlvq+bkK3QiYxCvj0Fhk8Vo2GwJXtqDSyCbv/ypBd2FT5mYC+N4Bc7qxhTWJTFGV7qb0Q08kHbyyN9gCvqAxgNM01z+42XMS11cvOmlaQIvHop6rcAFY8Yfh56rZtkK2oDIGufu/iwza/1WC7sKmzG1yz1mbtAKPx6L4wCnu9xnk966YNFJUZdyCLuALmr25Aw0SSQerH9g8kWatEsBB6f7i+7FnaIJE0vE88+01iQtwV+H/csdT0wK/hLvCM/258XsS1/pX9+PWrON3qpBt+ZFzjHsu/Pzec8VjUVRXCXz2J0n88GOnovOwlrI/p986Guswv6UOT1zzrsC+hkqh7x00ZwUOAEfPc38eT37+3ehoCq6ZCwC0Nqg9cP0BvFCDHdP+PgvhCrxI1yY2TusBnn33dvkdT2Hpdb+w7g4uV6G0qdckLsB9ofaqV7chjX7mcfOnbRuUk060qZAtiCNk2V7eOwAA6Lr5f41MVe7qH8HCVv+buGRbu3knBICv/GKrUc9B73AaQkzuTwdpYgWep1FS2GyrU/HCAF6EaxMbp6ymvAyOjls/1KNQsJ3pBd0rwJueRk8kHdyzwZlyAyIAXLgsesgM8zh2QTMi1f4dIVMSSQf/9ds/T7xtYqpyd18qsDnYwGS6Vv1+mfQc9A6NorU+Mu0GLgj1kWrU1VQZsQK3rU7FCwP4DIoN3tlsbSVaqHitvSEy4wt6d1dn3n+v0uimypdOk8i/ii6WTS8QiaSDF3YPID3u3xEyJehz1uVKJB1s2dmHRwLowqaY3A+9Zyhd8BSJ31obIkZUoXd3dU7rf25agWUxGMALmE0bUcDeVqJr1m71LF5bfb53pbliaxo9iNWyLYVsfk8hy2VyJiLoLmyKyc9B79Ao2gMcZJItkXTQMziK25/Yrn0bIR6L4tPvOHri7aAL+ILCAF6AVzV2MWxrYgJ4v6AUKl7LZWMaPag+4DZ0ZAu6G5fJmYiwsgMmPwe9Q2nP6Xt+UtsIQd0ozsbyI+cAAO785Ft9HSEbJgZwD6VO4MrHtlW41x9yoeK1XDam0f0uYFNsKGQLenVo8lnwsFbGJj8HPUOjgY4SVUzcRhjIdF1kI5cKU0zq/NIVS3DzxSejIeL9FPYOp60paCt3jKZSKI1uQsowVxAFbIrJqVMl6NWhyWfBw1oZT45WdZ+Dwwx5DhJJB07PMO5LOoGntE38W2AAr1Azpc4vXbEEN8TdP8DnvvjegkUgthS0lTtGM5vXij3ImcOzFUQBm2Jy6lTp7upETYBHyAAVwP4CAPD9j56iPXApYXVhAzDRdhgAfvzxU7U/B2r/P6zKeBP/FtTcg6a66hkeaS4G8ByFGpkAbjV2dhtRoHArUZNTx9m8isxmUzHa3dWJSJ5jKYOjY8bdzAS5MjA5darEY1GcsLgV1T5OIctncbt7TMvpMSf7EI9FcVYIXdiUuc1uqtqEdqphnw4w8W9hYgVu8TATe688INc/6L1v7VWNHY9Fcf2Dmz0Dv8kV2MDk8bF8WYfZdiLL93ykx6VxnY6C7kIGACvvfWaiUKy+wJaLLmNS4m3HzMOPPn5qYF9jo9MLAPi7H683qnNhR1Md5jXXYf21wXVhU+YbNNAk7JS2+llfd/8m9I2MYVFbPa48+zitvwMDqTFEqgXqargCrwgzrb7VBK58Vp231PMMtckV2EDh42OzvUPu9XgeTdr/BYIrYMuWvTPRM5TWXn2b7b4nX8Nmpw+/e35PYPugiaSDG3++ZeJtEyqQlV19wXdhU0waaKIjpR2PRXH1+9yxqvd8+nTtN3CDlk8iAxjApyiUPsqXOs8Wj0VxicdsbNPPhftxfCyX1wuBSfvgQRawKWvWbkVqzMxGJomkg5X3Bd8hLOijauXY1ZfCgpZwAvjvtuwGAFz/4LPaz0GHPYlMMWmk6MDIGJoYwCtHoVR3MY1MCgV4k8+F+3F8LJcN++BBFrApJlbfKmHtg5r6HKgubL/ZGlz2IftrXZ3YNPG27ixEdp+CMPb/FZMmkg1wBV45ym0jqhQKeiasOHL5dXwsVzwWzVscovbBTRBGYDGx+lYJK7Ca+ByE1YVNMfEc9NHzmwEAT69+T2iNTNpMWoEzgFeOQjOwi1l9K4WCnonFbH4eH8tl+j54GIHFxOpbJazAauJzEHYVtolZiN6hNKqrBFpCDGITI0UNmEg2mGIKvSLMdgZ2PvFY1PNcuInFbF4vIH7cIZu+Dx5GAZtq4tGYVX1uSiV6d1cnqkXw56DVc6C2aha21mlvZBJ2QDUxC9EzNIq2hgiECH4SmcIUur8KvpIIIdqEEBcLIf4/IcS/ZP6/PayLC8tsZ2B78apIN/FMuNf+tx8vLCbvg4dRwJZt3MBK9HgsiiPnNqK2uirwfdB4LIr/+zdvAQDcfHFMewVy2AHVxCxE73A4fdCztWS21ZhC94dnABdC/C2AJwG8A0AjgCYAZwLYkPlYxShnBnY+tkzlCmr/WzF5HzyMArbsr2VaJXoi0z7zpb2DqKkW+NrFJwe+D/rsjj4AwIe+s86IKuywurABZraUPTCURnvI2bCa6io019Wgb0RvAE8kHezuS+GO9fono5Wj0O3HNQCWSSl7s98phOgA8DiAHwd5YWEp1MSklOK1XFGPBiEqja57BQIEu/+tmLoPHmYK1bT9TzUZSt3ADI2OY+W9bs/+oH4vE0kHX/3V5A2LKhoL8msWEo9F8dDG17H22V0QgO9z0L2+5uj4QVxx9zO485NvxRFzGgP7WsXoGRrFwtb60L9ua30N+ob17YGrEbq5xycBPb+L5SiUQveKawczH6sI5c7A9mLDVK4g978VE/f+Cn39IK7LtOdAR0V02EVjxZjTXIc5TbV4+aZzQqvCnmdQO9VeDStwwD0LrjOFbnJfglIVCuBfAvCkEOKbQoirM/99C25a/UvhXF7wvFLa5TQxAexIowe5/63k2/sTcAvIdApzT9K0/U8dGQHTshAAsKd/JLQmLsrcJjPaqSaSDnb0DuPeECaR5X7dl/YM4tfP7dKWujbxd3G2PAO4lPJHAE4B8DqANIAUgN8CWC6l/GEYFxe0Qme/y2liMtPnMKEaPej9byW7YYQiAdyzwdH+HNTVTP76dzRGAi3gMqkSXUdGwLQsBADs7k9hQcgpZDXQZN+gvhV42JPIsr/uyns3YnT8YKhfN5eJv4uzVfBVREq5H0AtgE8COA1Ab+a/ihBED/BsJqfRw9j/VvIVhulMWakXkt6sNF5uejcIplSi68gImJaFAIDdIbZRVeZN9EPXtwLXtZ1hSjOb7q7OKTfvgP7fxdmacRkgpbwWwJsAfA/ARwG8IIT4shDi6ICvLXBBpc8Vk9PoYex/z/S1dKWsdO0Bm1KJrjIC6kUsjDaa6mu2ZRp5HKa5Cvu+Da9hZ98I7t7wWqip3F9s2gkB9/fhUEshm/I6EI9F8Ym3vQFAuG1kg1BUHk9KKQHszPw3BqADwN1CiK8EeG2BCjp9PtPn0p1GD2P/e6bPqauhC/eA3RexI+Y04uylh4VWwBWPRbHm/ScBAG758DJtL5iJpIOViY0Tb4edQg47dZ1LVwrZpNT1SUe47Uwe/OzbQvv9D8KMAVwI8U9CiA0AvgLgMQAnSik/DWAZgAsDvr7ABJ0+V0xMo4e1/62Y1tDlUN8DTiQdnH7Tw3hx9wAefXFvqD+DaIf7/To9+jJQTCHrmURm0jbKYMp9/avYRi5Z5gG4QErZJaW8S0qZBgAp5UEA5wZ6dQEKOn2umJhGD3P/GzCvocuhvAesVoE7ekcAuN2owlwFPr3dLaH59G1PMoUc0tfNpWsSmUktdQdUAM/zumSTYvbAr5NSvurxsef8v6TghZU+n+lz6kqjh7n/rZjW0KWuZvI3IMgKdGXixatB74uXzlVgIungiz99duJtppDD+br56JhEBrh/Bzdd8BcAgO995BRtqeuBQ2gFXnHCSp8rpqXRw9z/nulzh/3iNVmBPrmFEEYFOuC+eH3k9CMBALv6UlizdmvowUvnKtCUBhrdXZ2oqQqvjWr21zUhCwO4XdjCnkSmqIlk/Ronkg2mxlBdJaZVo9vG7qufpbDS50qhNHrYK9Cw978VUxq66F6Bfvt/X5p4W8cKVOeNlEkp5Le+cQ4E9KSQVSW+zn7oqgtbmJPIlImJZBr7oQ+MuINMdHz/fjrkAnjY6fOZPnfYldhh738rpjR00b0C1d1OVOcq0JQsDAA010fwxvlNobZRBdy/gy/89QkAgB9//FRtKeTe4TTaQp5EppgwUnQgNW59+hzQHMCFEGcLIbYKIV4UQlwVxtcMO32umFKJrWP/WzGhocuhvgJVq0A1iSvMM7AmpZB396ewoCX8QR7AZDOXfYP6mrn0Do1q6YMOTKbQ+zSm0AdSaQbwcgghqgH8XwDvBXA8gA8JIY4P+ut6vVgGlT5XTKnE1rH/rZgQwM48bv60DMyhtgKNx6Joa6jFB085IvTV540XnIgOA6qQd/ePYEFruF3YlIl2qhq7sfUOpdHRWKvla6vAqXMFPpgaR1Nd9cwPNJzOFfipAF6UUr4kpRwFcDuAvw76i3q9WAaZPle8KrHDOk6ma/9b0R3AEkkH92xwpmRgBIALl0UPmRVoIungrTc+jL0DKfx84+uhF9HFY1H81yXLAAD/cdHJWoL3fU++hu37h3H/Uzu0HGWb0+QGzv0a+6H3DulLoauZ4DqL2AZSY2jiCrwsUQDbs95+LfO+QJm4BxjWcTJd+9+K7kK2fAVsEvlT+0FQK9Cm2snnIMyhJqoC//UD7hnwvpFwz4Ari9vd1LWOI4SJpIOV94XfhS3bnEY1UlR3Cl3PChzIzATXWcSWGkOL5WfAAb0B3Otk1dQHCXGZEGK9EGL9nj3lv9CqF9Foe0PofXB1HyfTuf8N6C9kMyGFDwDjWTdRYQ41MaUT2BMv7wcAdN/9TOgrYBMKCWuqq9DRGNE2kezu9dsxODqO7z/2srZmOi31Ec0p9DE01TKAl+M1AEdkvX04gB25D5JS3iKlXC6lXD5/vj8rtXgsiseuOktLBarO42Q6978VnYVsulP4QCaAaBpqYsINTCLp4PP3b554O+wVsCnPQf/IGG5dty30AJpIOrg2sWnibR0ZiETSwct7B/HLZ/XMBE8kHew8MIK7Qh5kEwSdAfwJAG8SQrxBCFEL4IMAHtB4PaHQdZxM9/63ovMF1IQ9aJ3fvyk3MDqzALqfA7WNMZbJwoQdQHXeQAL6Z4K7Xz/8WehB0RbApZRjAP4RwFoAzwG4U0q5ufC/sp+u42S6978V3S+gYbdQzaXz+z/Ub2AAfV3YFN03MLqff93fvyndAP2i9Ry4lPLnUspjpZRHSym/pPNawqLrOJnu/W9FVyGbzhaq2XQGUVX/oQKYjjnIum/g4rEozjhmbuhd2BTdAVT386/7+9f99f12yHViM4GO42Qm7H8D+grZdN/5KyqIqqNE85vDOwudSDr4ytotGDso0VxXg+6uztCPcZmQBWhtqMWSuY2h18AA+gNod1fnRBMfJcznX/f3r/vr+40BXIOwj5OZsv+t6ChkM+nOOx6L4pN/9UYAwJ6BcIaa6B4jqkzcwGRuKOe3hN/MZXffCBa06GniovsGJh6L4n0nLAKgJwOh+/vv7uqcNsBEVzdAPzCAaxD2cTJT9r8VHcHUpDvvRNLB1379/MTbYRTSmJKBANwg8v2PnQoAuPH/hN+JbY/GNqq53egWaLiBOXxOA6qrBF668X2hZyB0j9WNx6L4xF++AYCeGxi/MYBrEPZxMlP2vxWvoBlkJb7OFqq5dJxFNikDAWQ1czkQ/tff3Z/CfE0rcMD9+//mpW43uq9dHH43uh6Nk8gA9/v/8gUnAgB+pGGgS+yIDgDAA//4ttBvYPzGAK5JmMfJTNn/VsKuxNfdQjXXoZ6BAIDHnt8LALju/s2hnsW984ltGEiN4Yd/eEXrGeC5mRoIHQNNDmhso6pMTiQLv53qQMpduOQrJrYNA7gmYQUx0/a/gfAr8XW3UM2lI5jq3nvMlkg6uFpDMxHdTWSyzVUTyQbC78bWOzyqbZCJMjGRTEMWcCDlvhZwmAnNWlhBzLT9b8WrEj/MLQRd6WMdwVTnGNFcuvbj16zdipTGJibZ2hsiqBJ6JpL1ZlLoOrWoFbiGfuhqQcNxolSWMI6TeX0uXfvfSpirUNPSxyqYttSFO9QkHouitT6CD50a7hjRXLpuqEy6kauqEpjTVKslha5zEpnSWq9vBT6YGkOVwLSbaBsxgGsU9HGyRNLJW+1e6GuHJcyGLiYVsGVLhzjURI0R3Tc4ip89E/4Y0Wy6bqhMupFLJB30DqXxkz+F3w+9d0h/Cl2twHWMFFWjRHUV8fmJAVyjoI+TrVm7NW+1u8h8bZ3CauhiWgGbEmYluiljRBVd+/G626gqOvuhq0lk33tU3yQyAKitqUJDpFpPCj01VhHpc4ABXKtCx8n8SKN7pQZl5mvrFkZDF9MK2JQw07kmnQEHJrcQVCX2vObaUPbj47EoTj6iDdVC7xlgXT8PEyaRZWttqNFShT7IAE5+8TpO5kca3ev4mNfXDFsYQcykfc9sYaZzTXwO4rEofnLZCgDAdectDS2INtTW4ITD27W0UVV0/Tx0TyLL1VIf0bYCb2IAJz8ElUY38fhYrjCCmEn7ntnCTCOb+hwsass0cwnxRmJ3X0pbG1VF18/DtBu51voabQG8pQLOgAMM4NoF1ZXN1ONj2cIoZDO1gE2lkec1u2nkuU3BpZFNOgOe7eHndkMAuOmhLaHtx+7u19cHXdH18zDpRi6RdLB5Rx8ee3Ff6Hvxg6kxNNUygJNPgujKZurxsWxBF7KZWsCmxGNR/NNZxwBwO3IFNdRE3Syo+i0T+j+rQi71swl6PzaRdHD6jQ+jZyiNB57eobUKf9pEupD6oeueRKaon706kx/2XvzACFPo5CO/u7KZfHwsV5CFbKYWsCmJpIMbH9oy8XZQL2SJpIN/+8UWHJRu2lLHGNFcYRZyTUxiy1Th92uuwgfcIP7tD7v90L/6gZNC+XnonkSm6C6qZAqdfFWoK9vqBzbn+ReFmXx8LFeQ+3Km7fnlcl/Igi0qMu0ImXIoV+Erk/3Qw2unGu1oQI2mSWSKzr/L+558DX0j+nvh+4UB3BBeXdl6h9Ml/5J5pc9NOT6WLcjJZCbt+eUTxguZqcHrUK/CB7L7oYfXja13OI32Rn2TyAB9f5eJpIOV922ceFv3MTo/MIAbotAvbykvtoXS56YcH8sW5FAXUwvYlDBeyEwNXqzCd7czItUi1HaqB4bSgY7tLYauIj4dY3yDxgBuiEK/vKU0dbEpfQ4EN9TF9AI2IJwXMlODF6vwASEy/dBDnEjWMzSKds1tVNXPviPTp2JBSEV8pt7MloMB3BDxWHTiFzpXKU1dbEqfK0FMJjO9gA0IZ6iJqcELcL//+/7hDADAFWcHV1innmf1PJhQhQ+4f9P7BkZx5/rXQtuP7R1Ke77OhCkei+LmD8YAAN+89C2h/CxMvZktBwO4QVadt7Sspi62pc+VIPbBbbrbDnKoiYlHyLId1lYPIYAdvSOBfh23jWo73rKkXeskNkVXP/QDw2m0NehdgSuTE8nCaafa3dWJ2pqpIc+Um9nZYgA3SLm90a9/cLNV6XMliH1wW+62g96XM/EIWbafPfM6BICvP/xC4KvQXX0jOCzT/U03XcWFbgpd/wocCH8meDwWxd+uOBKA3mN0fmIAN8xse6Mnkg56PFLRJqfPgWD2wU0vYFOCzBSYeoRMUdenEhBBrkKllNjZN4LDWs24gdORIbp7/XYMGTCJTGltyKzAQxwpeuLhbQCAX3/u7UZkYsrFAG6YQr3RC50JLxToTE6fK37ug9tQwKYEmSkw9QiZEtb1JZIOTr/pEQyNjuPuDdu1By4g/AyRaZPIAKBVrcBD7A6p5o9zGhkFolAavdCZ8EIpdtNWnfn4uQ9uQwGbEmSRmel1AGFcn6lZiLCLC02bRAYA9ZFq1NZUhTrQZDDFAE4BK7RizrcKvzaxMc8jXe0NEeNWnfn4uQ9ueuDKFmQluul1AGFcn6lZiLD7oZv6N9FaHwl1JvhgagxCAI211TM/2AIM4AYqdBfeO5yeErATSQe3rduW97ECwOrzl/p9eYHwcx/cq0jHlMCVTxCV6CYfIQPCuT5TAxfg/s5/9yPLAQBfufAvAr3RNvVmLuyRov2ZSWQ6O9H5iQHcQIXOhAPAreu2TQRxr8pzwPzitVx+7IPbMAc9V5CV6HU1ky9UHY0Ro6pu1Sp0fma855wAmrmYGrgU1Q99b8DNXEyZRJarpSEysS8dhsHUWMWkzwEGcGOtOq/wyvnWddtw9MqfeVaeA3YUr2XzYx/chjnouYJYJaq9396s9GTuTYIJ4rEoHvrnvwQAfPasY3z/GZmehVD90PcH3E41HoviPccvBGDWEarW+ppQi9gGUmNoqquM9DnAAG6smVbhADDuVWIgZwAAIABJREFUtfSG+We/8/FjH9wr6Jk0Bz1XEKtEU/d+8/n9825x4fUPPuv78Sa1yq/PNPAwJXApTbXVqKupCqUf+qK2BjTWVuPlm84x4ghVIulg/Sv78dT23tCOtQ2kxrkCp3B4dWYrxiUrlmj/Ay2VH/vgNu5/B7FKNHnvN1si6eDq+4I93hSPRfGmhS14+7HzjQhc2YQQmNtUG3gKHQD2D42iQ3MfdEVliNQ43bCOtQ2mxvK+xtiKAdxg8VgUl6xYUvK/a2+I4Ib4iQFcUfC89sGL6URn4/43EEwluul7v0pYmQK3iYsZXdiyJZIO9gykcO+TTuCr0J7BUXQ0mdGFTVeGaGDELWKrFAzghrshfiIuLSGI21R5no9XgClmoIuN+9/Z/KxEN33vVwk6U+A2cXkYe/pT+Pmm17Wf/86mVqHp8XD6oe8fShuzAteVIRrgCpzCVkoQtzF1nq1QJ7qZ7s69Vukm738rQVSim1yBrgSZKVABUg1K6TekiYsS9iq0d2h04ty5broyRAOsQicdZgriTbXVuPnik61NnSuzHehSaBKbaWnjfPxckdhSgQ4EmykwvZAv7FXo/kFz9sB1ZIjue/I1HBhO48d/fNWIXvB+YAC3yA3xE3HzxSejPetYVUdjBDdffDI2f+Fs41ZXszWbgS5r1m61chKb4ueKxPTAlU3t/6vz0POa/TsLbnohX5ir0Hs2bEf/yBh++IdXjAheuT/3+c3BdqJLJB2svG+yAZYJveD9wABumXgsiqdWvQev3HQOXrnpHCSve0/FBG5lNml0r9W5Lc1s/FyRmB64csVjUdzxyRUAgGvPOd63n5fphXxhrUITSQfXBFzpPxvxWBTf/vAyAMC/X3RSoH+nQY/t1YUBnIxTahq9UPrclmY2akXSmlVgM9tKdNMDVz6L2txr23HAv5sM0wv5wlqFmjjIRGnNZBP7A26nattNbbEYwMlIpaTRvdrJ2pI+z5YeL78S3fTAlc+vnt0FIYCv/GKrbyneiSYuETObuADuNX7/o6cAAL4c0LWZHLwmR4oG207VxpvaYjCAk5GKnYueSDqe7WRtSZ8rfu5d21CBrqiiO5m5d/EzxRuPRXHMgmYjm7goc5vdFfj+wWCauZgcvFob3IxT0ANNurs6UVszNdyZflNbDAZwMlKxc9Gvf3D6eFXFlvS54sdKyaYKdCXoorsdvSOIdpj7uzC3ye2HvncgmHaqpg4yUddRUyUC74cej0Xx4dPcUzwm9YIvV+UciKOKE21v8CxOW/3AZqx/dX/BYS4mvECVYrHH91vKSqlQMDT1xSqoFG8i6eArv9iC/YOj+OnTO3DqUXOMfA4aaqvRWFuNfQEF8Hgsil9u3omfb9oJAff3qbur04jnQgiB1oZIKCNFj1/cBgD4zb++A0fNawr864WBK3Ay1kxz0W/1mIMOuO1kTXiBKkW+vWsB4Mzj5hf9OUze7/QSRIp3oonLAbeJS59hTVxyzW2uDSyFDgCHtTWgyaBBJtla6mtCGSk6kBqb+HqVggGcjFXMRDYvNraTjceiuHBZdMrevwRwzwan6MBj8n6nlyCK7mw6C59IOth5YASJp3YEdka7d2gUHYZ0YcvVWh8JZaSoqnRnK1WikMw0Fz0fG1ffym+27Jm2919K4DnzuPnTiv9M2e/0oqrFF7S4e8F+FN3ZkokIqx/6foPaqOZqbahBXwgr8P7UGGprqlBXw3ngRKEodRVu+zCXcgJPIungng3OlBsAAeDCZVHjb2jisSh+86/vAAD83V++sezrtSUTEVamoGdwFO2GtFHNFd4KfAwtFdQHHWAAJwuUMhfd9mEuXgGmrWHmm5h8wUDCXdXb4FfP7kKVcL+PclPJtpyFDytTsH9oFHNmuR0VpETSwe+e34MXdg8E3uJ1YGSsova/AQZwskCxc9EvXbHE+mEu3V2diFRNv10ZHB2b8cXNlrRxPiqVfNCns+AqLa/O/pp6bCisTEHPYNq4PXD1Mx8adW86g27xWmmjRAEGcLLETNPYKiF4A27gyfcikx6XM6ZVbUkb5xNEKjkei2JRWz3OO2mxcZXXShiZgrvXb8dAagw/eMyMQSZK2IWG/SNptNSZl4UoBwM4WaPQNLZKCN5Kr8fZ9plW0jYWsCl+Zw8SSQen3/QwXt03hN9s2W1M0MoV5DQ2wH0erk2YN8gECD9j1D9SeStwLd+NEGINgPMAjAL4M4CPSSl7dVwL2SUeM78gq1xeDV0K7YPbXMAG+NPERlGpWbW6G0i5Z8ABM1vrqnav537jUdwQPxFnn3CYb5+70CAT3c+Fnz/zYvRzD9w3vwJwgpTyLwA8D2ClpusgMs5s9sFtL2DzM5Vs0xlwZV6ze4Run8/NXEyuiwi70HAgxSp0X0gpfymlVAf/1gE4XMd1EJloNvvgJr9QF0Olkhe2uoGsvWH2Z8FtfC46mtzsit/tVE2ui1A/83nNwWwfZLvvyddwYDiNH/3xVaPqAMplwh74xwE8pPsiiExS6j54u8cRIRNeqIsVj0VxZddxANxWuWvWbp3VC63JQctLXU01WuprsH/Q3wBu8iATwP2Z//cnTgMAfPGvTwgkeCeSDlbet3HibZPqAMoVWAAXQvxaCLEpz39/nfWYawCMAbitwOe5TAixXgixfs8eO9KBROUq5Tx4IulgIE8nq0i1MOaFuhiJpINrfCi4suUMeK55zXXYO+BvCj0ei+Idx7q99E2dwtWa+Z0OaqDJmrVbp03kM31LpViBbQhIKd9V6ONCiI8AOBfAO6WUXpMjIaW8BcAtALB8+XLPxxFVku6uTnTf9TTSB6f+yqt98OwX4DVrt057HAA01dYY9UI9E78mqanHXnH3MxgdP4ioQdO3vCSSDpyeYby8dxDJbY/4er3zWuowv6UOT1xT8CVZm9bMdlHfcDDtVG3cUimWlhS6EOJsAFcCOF9KOaTjGohMVso+uNfI1QMhtKf0k58vtPFYFHOaanHBW6LGngFXVNX86Li7SvQ7xbt3YHTimJqJmmprUCWCW4HbuKVSLF174P8JoAXAr4QQTwkhvqXpOoiM5bUPnh2wE0nHs82sbS9Qfr3QJpIOTr/xYezsG8Gvnt1l/F5n0FXz+wZSE1XuJqqqEmipjwQ2UrS7qxO11VNDnQ1bKsXQVYV+jJTyCCnlyZn/PqXjOohM5hW4BDARlK5/cPO06WXqMba9QPmxd507B7zf8DngQPAp3n2Do5jbbO4KHHBndAc10CQei+JvTjsCgLl1ALNVWYfiiCpId1cn/uWOp6YFaAlg9QObAQA9Hqt0CTOblhSirvfGh57Drr4U2hsiWH3+0pK+D7/20cMUdEOTvf0pzG0ydwUOZCaSBZRCB4DjF7UBAH5/5Zk4vKMxsK8TNhOOkRFRHvFYNO/qGnCPWV2TdTQmV9Sy9LkSj0VxVRlHyWwsWAqyav7OJ7ZjcHQc33/sZaPPP7c21ARWxAa4s8ABsBc6EYWnUCAeHB33/Jht6XMlkXRwdRlHyWwsWFINTeZn9qnnNvnT0CSRdPD5+83sg54r6BV4f+ZzV1ovdAZwIoPNJhC3N0SMTRfPpNyCLlvPgMdjUdz3mdMBwLcjZGvWbkXKow+6SRJJB4++uBdbdvYHliUYGBlDY201qvO0KLYZAziRweKxKDo8uqx5WX3+0oCuJnjlpsDVala9UNtUsDS/xV2B7+n3p5mLDdsJYc0Er8RBJgADOJHxVp1XfEC2efUNlJ8CTyQdfGXtFowflGipqzG+gUu2uppqtDVEsNunAG7DdkJYg2cGUmNorrBBJgADOJHxil2FC9i9+gbKS4FPHCHrzRwhS5l/hCxbIulgMDWG/17nz8CN7q5O1FSZ2wcdCC9L0DeSRkt9ZRWwAQzgRFZYdd5Sz4YtyiUrlliz2vSiUuBtDZOrpfpIcS9TNo4RVdTNx1imJa4fqeR4LIq3HTMXAuaefw4rSzCQYgqdiDSJx6K4ZMUSz49fumIJboifGOIVBWt0bPIAXc9QuqhgZsOer5egbj46mupw+JwGvHzTOUa2lA2r6JB74ESk1Q3xE3HzxSejPWsiWUdjBDdffHJFBe/ZBjMb9ny9BHXzsXfA7CYuYc0EHxipzD3wyvuOiCpYPBY1bhXlt9kGs+6uTqy8d+OU4G/anq+XoLqx7RsYxeL2+rI+R9DisSiOWdCMc7/xKL70f05E19LDfP8a/dwDJyIKXjkr6bqayUqBjsaIcXu+XoJKJZu+AlfUnPsg+qHfu+E1DI6O43uPmt2NbjYYwInIKLMJZqoIrDerHedI+qDn400TRCr54EGJ/RYMMgHcTmwA0OfzRDK3s99ky2GTu9HNBgM4ERllNpXoNlegK/FYFLdftgIA8Plzjy87c3D7E9swdlDiv377Z+NXnqrFab/P7VTXrN067UbOtt+LQhjAichIpVSi21yBnm1+i7tfXW43tkTSwfUPPjvxtukrz+oqgeY6/weaVMrvhRcGcCIyTqkrapsr0LO11tegtqaq7G5stvRBz9ZaX+P7QJNK+b3wwgBORMYpdeVk6xCTXEIILGipK3sFbuPKs7Uh4nsRW3dXJ2qrp4Y5G38vvDCAE5FxZrNysrUCPVsi6WBX3wjuSzpl7VvbuPIMYqRoPBbF+5e5vwOmdqMrBwM4ERkn34paADjzuPnTHmt7Bbqivo/0ePntVG3og54tkXTwzGu9WPfSft8L7t60sAUAkLzu3UZ2oysHAzgRGScei+LCZdEp/d8lgHs2ONNe3CuhAh3w9/uIx6I4MdqG6iph/MpT3biMZPbs/S64O5BJy7MTGxFRSH6zZQ9kzvtUQMsORDbu9+bj9/cRqanCsiM7cOcn31rOZQWu0I2LHzccfcNuG9Wa6spbr1bed0REFaHYgNbuMWrV5P3efPzet97dN4KFrWa3UQWCvwHrG0mjtQIHmQAM4ERkKK/A1ZY1zCWRdDCQp3tXpFoYu9/rxc9KeikldvWlsLDF/DaqQRfc9Q2n0dpQeX3QAQZwIjJUd1cnIlXTp6APjo5N7I+uWbsV6YO5iXagqbbGyP3eQlQHuoWtbtBtn2UlfSLp4PSbHsFwehx3rt9ubPMWJegjgH0jDOBERKGKx6ITLTazpcflRGFXvglewGThkm3isSh+130mAOATZ7xhVsF75b0b8fqBEQBub3GTO7ABkzcuc5rcnu3zW+p8Lbg7MDw20Wu90jCAE5GxeofyB2KndxiJpIPp63OXbfvf2eoj1ZjTVIsdmSBcClsr8uOxKG758DIAwFc/cJKv2RM3hc49cCKiUHkFYgHg+gc3T6tSVx+zbf87WyLpoH8kjZ/8aVvJZ6JtrshXxYi9PmdP3CI2rsCJiELV3dWZd5Ut4Q44yUcC1u1/K+U2c7GxA5vS1uCm0A8Mjfr2Oe/b8Br6R8bwwz+8YvxEttlgACciY8Vj0byr7EKiFgQrL+WmwG3uCa9OF/hVv5BIOlhZwbPAAQZwIjJcKQHZ9vR5uSlwVRCm5qeb3IEtV21NFRprqz3rHkpV6bPAAQZwIjJcKQHZ5vQ54E8KPB6L4s2LWnHGMXOt6/3d3hDxbQ/c5nqAYjGAE5HR4rEoOjy6reWyOX0O+JcC392XwsIW87uw5WprrPVtBW5zPUCxGMCJyHirzls642NsT58D/jRzOXhQYnf/CBa22RfA2xsiODDsTxGbOwvcnolss8EATkTGK2YVfsmKJVali73EY1H87xVuM5ePnV5aM5dE0sHp//YI0uMS//P4q1YVbCWSDpLbe/DEKz2+VIzHY1Fc8JbDAVTmLHCA08iIyBKrzluK7ruezts69dIVS3BD/EQNVxWMuppqzGuuxesHit+vVUfQVBX7gWG3Cxtgfl3AxEjR9NSRokB5137MgmYAwNOr31ORZ8G5AiciK8RjUaz5wEloz+pr3dEYwc0Xn1xRwRtwA9qB4TRuf2J70atRW7uwAcFde99wGkIAzbWVuVatzO+KiCpSPBY1fjVZLq9mLkDh1ajNVddBXXvfyBha6mpQlWcoTiXgCpyIyCCzXY3aXHUd1LX3DafRVuQJBhsxgBMRGWS2q1Gbu7AFde0Hhiu3DzrAAE5EZJTZrkbVEbRI5uiUTVXX6trnNbv90Oc11/py7ZU8yARgACciMko5q9F4LIo5TbX4wLLDrevCFo9F8eOPnwYAuCF+gi/X3jc8VrGjRAEGcCIio6jV6KJMI5bW+pqiVqOJpIPTb3wYu/pS+MXmnVadAVcmRor61I2NK3AiIgpVPBbFH1e+E+2NEZx70uKigvfKezdix4ERAED/yJiVk7f8ngneN5yemHJWiRjAiYgMlEg6GEyN4X8e3zbjWXCbz4Bna4hUo7a6ypcV+D0btmNwdBzfffTlipwFDjCAExEZx+ssuFcQsvkMeDYhBNoay++Hnkg6uCaxaeLtSpwFDjCAExEZp9QVtc1nwHO1N0TKXoEfCrPAAQZwIiLjlLqitvkMeK72xvIDeKVkJGbCAE5EZJhSV9Sqcr2myr4z4NkSSQcbXzuAP760r6x960rKSBTCAE5EZJhSV9SJpIM1a7dg7KBEU537OBuD98p7N2JkbOpEstkE8e6uzomGNoqtGYlCGMCJiAyjVtTZk9fqI/lfrlXgc3rdI2SDqXErC7b8rKSPx6I458RFACp3FjjAaWRERMZKjU0WYvUMpfNOJSsU+GwKWH7vWx8xpxFVAnjxS+/jNLIgCCH+VQghhRDzdF4HEZFpil2RVkrBlt/71j1Do2hriFRs8AY0BnAhxBEA3g1gm65rICIyVbGBuVIKtvyupO8ZSqOjqdaPSzOWzhX41wBcAUBqvAYiIiMVG5gr5QhZ7kSyuU3lTSTrHRpFRyMDuO+EEOcDcKSUT+v4+kREpislMNfVTL6UdzRGrC3YiseiuP2yFQCA6847vqzvYf9gGh2NldsHHQiwiE0I8WsAh+X50DUArgbwniI/z2UALgOAJUuW+HZ9REQmU8Fr9QOb0Ds8BmB6JbqqQM/eK8/tQGabOU11AIB9A+W1U+0dGsXSxa1+XJKxAluBSynfJaU8Ifc/AC8BeAOAp4UQrwA4HMCTQoh8wR5SyluklMullMvnz58f1OUSERkpNTa5y/j/2rv3IDvr+o7j789md3OFhJBwOwGCXEKHkLI0AmNGaxAM9cZKi1MHWyxOVQqOykyQONx0GBsN2IvT0kGLYKvcBFZQS8CCoKhAkg1JMAQpt2QTSELuCbmQfPvHeU48JLvZsyfZfZ7nPJ/XTCbnPPvss9/f/mb3s7/n8vtV7kSvPCLWKIuYVBs1tIUmlW9C2x9rt2xv+BH4gJ9Cj4iFEXFYRIyPiPHAMuD0iHh9oGsxM8uy3gK6Ue5Ar9bUJA4Z1sqbm+sP8HvmLGXrjl1891eNuxIZeCIXM7PM6i2gG+UO9D2NHt7KmjpPoXd0dnFNAVYigwwEeDISX512HWZmWdNTEI9MZmibevJY9nzKOY93oO9p9PBW1tQ5Ap81e8nu6Vgr8n5ZoSepB7iZmXVv+rQJtHQzEcnm7W9zdcdC7p3b9Y7ncAX85Z+VcnkHerXRw1t5c/O2uj63ES8r9MQBbmaWUe1tJUYM2fthoR07gzueWrrX9fEAHnt+1QBV1z86Ort4/IVV/N+qzXVdv27UywrdcYCbmWVYT2tj74zu58DK80iz8ljclu3lP0zquX5dlJXIwAFuZpZpfR055nmkeSAei2tvK3HunxwONPZKZODVyMzMMm36tAl8+a75Nc05rWT/vDpQ16+PGDmUEYObWfS1aQeirMzyCNzMLMPa20o1LxgRkOuR5oG6fr12y3ZGNfgkLuAANzPLvFKNAVbrfll1oBZmWVuAhUzAAW5mlnm1BlieT5/DH1ckO+LgIUD5efd6rl8XYSlRcICbmWVee1up13m9Rw1tyfXp84r2thJPXDkVgEumHFdXm9YVYB50cICbmeXCdR89Za9Z1yoEXP+xUwaynH7V2tzEQUOaWVPnZC5rN/sUupmZZUR7W4mLzjpmrxAXcNFZxzTE6LvaocPrW9DkvrlL2bD1bW77zSsNvZAJ+DEyM7PcuKH9VCYfO5pZs5ewfN1bHDVqKNOnTWi48IbydKp9XVK0o7OLr3azkAnk++78njjAzcxypL0t/3Od96ajs4vfL9/A1rd3MWXmozX/kTJr9hK27uh+IZNG/J75FLqZmWVGZTrVyopifZlOtUgLmYAD3MzMMmR/plMt0kIm4AA3M7MM2Z9RdJEWMgEHuJmZZcj+jKLb20pMnTAWaPyFTMA3sZmZWYZMnzaBGfctfMdp9L6Mog8dMZgxI1qZc/W5/VViZjjAzcwsMyqj5W/8fDErN27jkGEtXPfRU2oeRa/etJ0xIwb3Z4mZ4VPoZmaWKe1tJR654s8BuGzqCX06Bb560zYHuJmZWVoOHtJMa3MTKzf2bTrVcoA3/jSq4AA3M7MMksRhBw1mVV8DfKNPoZuZmaVq7EGDWblxa8373/X0a7y1Yyff+/XLDT8POjjAzcwsgzo6u1i8YgNPvvhmTWHc0dnFtQ88t/t9X2ZwyysHuJmZZcru6VR31D6d6qzZS9j2dvfzoDcqB7iZmWVKPdOpFm0edHCAm5lZxtQTxkWbBx0c4GZmljH1hPH0aRNobirOPOjgADczs4yZPm0CQ1sGvWNbb2Hc3lZi0riRDGpSIeZBB0+lamZmGVMJ3W8+9Dwr1m/l4CHNfP38ib2GcWtzE21Hj+LHl75nIMpMnUfgZmaWOe1tJX474wOMGNzMBaePq2kk/caGbRw+csgAVJcNDnAzM8usI0YO4fX1vU/mcv+8ZbyyejM/W7CiEJO4gAPczMwyqqOzi6VrtvDQc6/vM5Q7OruYcf9CInlfhElcwAFuZmYZVJnMpTI5y75CedbsJbsnfalo9ElcwAFuZmYZ1JfJXIo4iQs4wM3MLIP6EspFnMQFHOBmZpZBfQnlIk7iAg5wMzPLoL5M5tLeVmJSqViTuIAncjEzswyqhO+3Zj/P8nVbGTG4mRva9zGZi+DM40bzo78/awCrTJdH4GZmlkntbSV+c9UHOGb0MKaefFi34d3R2cWUmY8y77V1zF+6ruEfHavmADczs8zq6OzijQ1befDZ5Xs9C1551KwrubFty/adhXj+u8IBbmZmmdTbs+D1rBveSBzgZmaWSb0FdFGf/65wgJuZWSb1FtBFff67wgFuZmaZ1FtA17NueCNxgJuZWSb1FtDtbSX+8YJTaR1UjrKiPP9d4QA3M7NMqgT0qKEtu7cNafljbHV0djFr9hK279zFsNZysBclvMETuZiZWcZV7kIHWLtlBzPuW8icV9dw79yu3Te5VR4hAwoT4h6Bm5lZZvV0J/odTy0t9CNk4AA3M7MM6+lO9J0Rfdq/ETnAzcwss3q6E13dbi3OI2SQYoBL+oKkJZKek/SttOowM7Psmj5tAi1Ne8d1d+PvlkEqzCNkkFKAS5oKnA9MiohTgBvTqMPMzLKtva3EiCG13W89vLW5MDewQXoj8EuBmRGxDSAiVqZUh5mZZdy6LTtq2m/9W7Xt1yjSCvCTgPdKekrS45LenVIdZmaWcbVe1y7S9W/oxwCX9AtJi7r5dz7l588PAc4CpgN3S+r2ngRJn5U0R9KcVatW9Ve5ZmaWUdOnTejxprUKJfsVSb8FeEScExETu/n3E2AZcF+UPQ3sAsb0cJxbImJyREweO3Zsf5VrZmYZ1d5W6vamtWpBcSZwqUjrFHoHcDaApJOAVmB1SrWYmVnGlXo5Pd7bxxtRWgF+K/AuSYuAO4GLI3p4Kt/MzApvX6fRi3j6HFKaCz0itgOfSuNrm5lZ/rS3lZjz6hp++LvX3nE6XcBFZx1TuNPn4MVMzMwsJ25oP5XJx45m1uwlLF/3FkeNGlq4FciqOcDNzCw32ttKhQ3sPXkudDMzsxxygJuZmeWQA9zMzCyHHOBmZmY55AA3MzPLIQe4mZlZDjnAzczMcsgBbmZmlkMOcDMzsxxygJuZmeWQA9zMzCyHHOBmZmY55AA3MzPLIQe4mZlZDjnAzczMckgRkXYNNZO0Cnj1AB5yDLD6AB4vTW5LNrkt2dQobWmUdoDb0pNjI2Jsdx/IVYAfaJLmRMTktOs4ENyWbHJbsqlR2tIo7QC3pR4+hW5mZpZDDnAzM7McKnqA35J2AQeQ25JNbks2NUpbGqUd4Lb0WaGvgZuZmeVV0UfgZmZmuVTYAJd0nqQlkl6UdFXa9ewPSa9IWihpvqQ5adfTF5JulbRS0qKqbaMlPSLpD8n/h6RZY616aMv1krqSvpkv6UNp1lgLSUdLekzSYknPSfpisj13/bKPtuSxX4ZIelrSs0lbvpZsz2O/9NSW3PULgKRBkjol/TR5PyB9UshT6JIGAS8A5wLLgGeAT0bE71MtrE6SXgEmR0TunqGU9D5gE/CDiJiYbPsWsCYiZiZ/XB0SEV9Js85a9NCW64FNEXFjmrX1haQjgSMjYp6kg4C5QDvwaXLWL/toyyfIX78IGB4RmyS1AL8GvghcQP76pae2nEfO+gVA0hXAZODgiPjIQP0OK+oI/AzgxYh4KSK2A3cC56dcUyFFxBPAmj02nw/cnry+nfIv3MzroS25ExErImJe8nojsBgokcN+2UdbcifKNiVvW5J/QT77pae25I6kccCHge9VbR6QPilqgJeApVXvl5HTH+pEAA9Lmivps2kXcwAcHhEroPwLGDgs5Xr21+WSFiSn2DN/erOapPFAG/AUOe+XPdoCOeyX5FTtfGAl8EhE5LZfemgL5K9f/hm4EthVtW1A+qSoAa5utuXyr7/ElIg4HfgL4LLkVK5lw83A8cBpwArgpnTLqZ2kEcC9wJciYkPa9eyPbtqSy36JiJ0RcRowDjhD0sS0a6pXD23JVb9I+giwMiLmpvH1ixrgy4Cjq96PA5anVMt+i4h0rMtxAAAE/0lEQVTlyf8rgfspXyLIszeSa5eVa5grU66nbhHxRvKLahfwXXLSN8l1yXuBH0bEfcnmXPZLd23Ja79URMQ64JeUrxnnsl8qqtuSw36ZAnwsuQ/pTuBsSf/NAPVJUQP8GeBEScdJagX+Gngg5ZrqIml4cnMOkoYDHwQW7fuzMu8B4OLk9cXAT1KsZb9UfogTHycHfZPcYPSfwOKI+HbVh3LXLz21Jaf9MlbSqOT1UOAc4Hny2S/dtiVv/RIRMyJiXESMp5wjj0bEpxigPmnuj4NmXUS8LelyYDYwCLg1Ip5Luax6HQ7cX/49RTPwo4h4KN2SaifpDuD9wBhJy4DrgJnA3ZI+A7wGXJhehbXroS3vl3Qa5Us0rwCfS63A2k0B/gZYmFyjBPgq+eyXntryyRz2y5HA7clTNE3A3RHxU0m/JX/90lNb/iuH/dKdAflZKeRjZGZmZnlX1FPoZmZmueYANzMzyyEHuJmZWQ45wM3MzHLIAW5mZpZDDnAzM7MccoCbZYykUZL+oer9UZJ+3E9fq13Stf1x7HpI+qWkyfv4+I2Szh7ImsyyygFulj2jgN0BHhHLI+Kv+ulrXQn8ez8duz98B7gq7SLMssABbpY9M4HjJc2XNEvSeEmLACR9WlKHpAclvSzpcklXSOqU9DtJo5P9jpf0ULJC3a8knbznF5F0ErCtso68pAslLZL0rKQnkm2DkhqeSVaI+lzV518paWGy/8xk22lJHQsk3V9ZTSoZWX9T0tOSXpD03mT7UEl3JvvfBQyt+rq3JfUslPRlgIh4FThU0hH99c03y4tCTqVqlnFXAROTlZoqy2BWm0h5WcwhwIvAVyKiTdI/AX9LeXnDW4DPR8QfJJ1JeZS956nnKcC8qvfXAtMioqsyTzXwGWB9RLxb0mDgSUkPAydTXuP4zIjYUvnDAfgB8IWIeFzS1ylPJ/ul5GPNEXGGpA8l288BLgW2RMQkSZOq6jkNKEXExOR7UKmHZJ8plBcoMSssB7hZ/jwWERuBjZLWAw8m2xcCk1ReOvM9wD3JHPkAg7s5zpHAqqr3TwK3SbobqKxA9sHkmJVT+COBEymH7/cjYgtARKyRNBIYFRGPJ/veDtxTdfzKMecC45PX7wP+NTnGAkkLku0vAe+S9B3gZ8DDVcdZCRzV3TfGrEgc4Gb5s63q9a6q97so/0w3AesqI/h9eItyIAMQEZ9PRusfBuYni0qI8oh6dvUnSjqP8oIT9dS9k3f+7tnrOBGxVtKfAtOAy4BPAJckHx6S1G5WaL4GbpY9G4GD6v3kiNgAvCzpQigvqZmE4Z4WAydU3kg6PiKeiohrgdXA0ZRX7LtU5TW1kXSSysvWPgxcImlYsn10RKwH1laub1NeBexx9u0J4KLkGBOBScnrMUBTRNwLXAOcXvU5J5HxZSbNBoJH4GYZExFvSnoyuXHtf4B/q+MwFwE3S7oaaAHuBJ7dY58ngJskKcrLEs6SdCLlUff/JvsvoHy6e57K5+NXAe0R8VAyQp8jaTvwc8rLdF4M/EcS7C8Bf9dLnTcD309Onc8Hnk62l5LtlUHGDIDkD4kTgDl9/YaYNRovJ2pWYJL+BXgwIn6Rdi21kPRx4PSIuCbtWszS5lPoZsX2DWBY2kX0QTNwU9pFmGWBR+BmZmY55BG4mZlZDjnAzczMcsgBbmZmlkMOcDMzsxxygJuZmeXQ/wNAdsH7ENM6GQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import json\n", "from numlabs.lab4.example.do_example import get_init,euler4\n", "#\n", "# specify the derivs function\n", "#\n", "def derivs(coeff, y):\n", " f=np.empty_like(y) #create a 2 element vector to hold the derivitive\n", " f[0]=y[1]\n", " f[1]= -1.*coeff.c1*y[1] - coeff.c2*y[0]\n", " return f\n", "#\n", "# first make sure we have an input file in this directory\n", "#\n", "\n", "coeff=get_init()\n", "\n", "#\n", "# integrate and save the result in savedata\n", "#\n", "time=np.arange(coeff.t_beg,coeff.t_end,coeff.dt)\n", "y=coeff.yinitial\n", "nsteps=len(time) \n", "savedata=np.empty([nsteps],np.float64)\n", "for i in range(nsteps):\n", " y=euler4(coeff,y,derivs)\n", " savedata[i]=y[0]\n", "\n", "theFig,theAx=plt.subplots(1,1,figsize=(8,8))\n", "theAx.plot(time,savedata,'o-')\n", "theAx.set_title(coeff.plot_title)\n", "theAx.set_xlabel('time (seconds)')\n", "theAx.set_ylabel('y0');\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "##### problem coding A\n", "\n", "\n", "\n", "As set up above, do_example.py\n", "solves the damped, harmonic oscillator with the (unstable) forward Euler method.\n", "\n", "1. Write a new routine that solves the harmonic oscilator using [Heun’s method](#eq:heuns)\n", " along the lines of the routines in [lab4_functions.py](https://github.com/phaustin/numeric/blob/master/numlabs/lab4/lab4_functions.py)\n", "\n", " Hand in a fresh notebook with the code and a plot.\n", "\n", "\n", "##### problem coding B\n", "\n", "1. Now solve the following test equation by both the midpoint and\n", " Heun’s method and compare. \n", " \n", " $$f(y,t) = t - y + 1.0$$ \n", " \n", " Choose two sets\n", " of initial conditions and determine if \n", " there is any difference between the two methods when applied to\n", " either problem. Should there be? Explain by analyzing the steps\n", " that each method is taking.\n", " \n", "2. Add your answer as new cells to the problem A notebook\n", "\n", "\n", "##### problem coding C\n", "\n", "1. Solve the Newtonian cooling equation of lab 1 by any of the above\n", " methods. \n", "\n", "2. Add cells that do this and also generate some plots, showing your along with the parameter values and\n", " initial conditions." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mathematical Notes \n", "\n", "\n", "\n", "\n", "\n", "\n", "### Note on the Derivation of the Second-Order Runge-Kutta Methods\n", "\n", "A general s-stage Runge-Kutta method can be written as,\n", "\n", "\n", "$$\n", " \\begin{array}{l}\n", " k_i = h f(y_n+ {\\displaystyle \\sum_{j=1}^{s} } b_{ij}k_j, t_n+a_ih), \n", " \\;\\;\\; i=1,..., s\\\\\n", " y_{n+1} = y_n + {\\displaystyle \\sum_{j=1}^{s}} c_jk_j \n", "\\end{array}\n", "$$ \n", " \n", " where\n", "\n", "${\\displaystyle \\sum_{j=1}^{s} } b_{ij} = a_i$." ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0 }, "source": [ "In particular, an *explicit* 2-stage Runge-Kutta method can\n", "be written as, \n", "\n", "\n", "$$\n", " \\begin{array}{l}\n", " k_1 = h f(y_n,t_n)\\\\\n", " k_2 = h f(y_n+ak_1, t_n+ah)\\\\\n", " y_{n+1} = y_n + c_1k_1 +c_2k_2\n", " \\end{array}\n", "$$\n", "\n", "where \n", " \n", "$b_{21} = a_2 \\equiv a$. \n", " \n", " So we want\n", "to know what values of $a$, $c_1$ and $c_2$ leads to a second-order\n", "method, i.e. a method with an error proportional to $h^3$.\n", "\n", "To find out, we compare the method against a second-order Taylor expansion,\n", "\n", "\n", "\n", "$$\n", " y(t_n+h) = y(t_n) + hy^\\prime(t_n) + \\frac{h^2}{2}y^{\\prime \\prime}(t_n)\n", " + O(h^3)\n", "$$\n", "\n", "So for the $y_{n+1}$ to be second-order accurate, it must match the\n", "Taylor method. In other words, $c_1k_1 +c_2k_2$ must\n", "match $hy^\\prime(t_n) + \\frac{h^2}{2}y^{\\prime \\prime}$. To do this, we\n", "need to express $k_1$ and $k_2$ in terms of derivatives of $y$ at time\n", "$t_n$.\n", "\n", "First note, $k_1 = hf(y_n, t_n) = hy^\\prime(t_n)$.\n", "\n", "Next, we can expand $k_2$ about $(y_n.t_n)$, \n", "\n", "\n", "\n", "$$\n", "k_2 = hf(y_n+ak_1, t_n+ah) = h(f + haf_t + haf_yy^\\prime + O(h^2))\n", "$$\n", "\n", "\n", "\n", "However, we can write $y^{\\prime \\prime}$ as, \n", "\n", "$$\n", " y^{\\prime \\prime} = \\frac{df}{dt} = f_t + f_yy^\\prime\n", "$$ \n", "This allows us\n", "to rewrite $k_2$ in terms of $y^{\\prime \\prime}$,\n", "\n", "$$k_2 = h(y^\\prime + hay^{\\prime \\prime}+ O(h^2))$$\n", "\n", "Substituting these expressions for $k_i$ back into the Runge-Kutta\n", "formula gives us,\n", "$$y_{n+1} = y_n + c_1hy^\\prime +c_2h(y^\\prime + hay^{\\prime \\prime})$$\n", "or \n", "$$y_{n+1} = y_n + h(c_1 +c_2)y^\\prime + h^2(c_2a)y^{\\prime \\prime}$$\n", "\n", "If we compare this against the second-order Taylor method,\n", "we see that we need, \n", "\n", "\n", "$$\n", " \\begin{array}{l}\n", " c_1 + c_2 = 1\\\\\n", " a c_2 = \\frac{1}{2}\n", " \\end{array}\n", "$$\n", " \n", "for the Runge-Kutta method to be\n", "second-order.\n", "\n", "
\n", "If we choose $a = 1/2$, this implies $c_2 = 1$ and $c_1=0$. This gives\n", "us the midpoint method.\n", "\n", "However, note that other choices are possible. In fact, we have a\n", "*one-parameter family* of second-order methods. For example\n", "if we choose, $a=1$ and $c_1=c_2=\\frac{1}{2}$, we get the\n", "*modified Euler method*,\n", "\n", "\n", "\n", "\n", "$$\n", " \\begin{array}{l}\n", " k_1 = h f(y_n,t_n)\\\\\n", " k_2 = h f(y_n+k_1, t_n+h)\\\\\n", " y_{n+1} = y_n + \\frac{1}{2}(k_1 +k_2)\n", " \\end{array}\n", "$$\n", " \n", "while the choice\n", "$a=\\frac{2}{3}$, $c_1=\\frac{1}{4}$ and $c_2=\\frac{3}{4}$, gives us\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
Heun's Method
\n", "\n", "$$\n", " \\begin{array}{l}\n", " k_1 = h f(y_n,t_n)\\\\\n", " k_2 = h f(y_n+\\frac{2}{3}k_1, t_n+\\frac{2}{3}h)\\\\\n", " y_{n+1} = y_n + \\frac{1}{4}k_1 + \\frac{3}{4}k_2\n", " \\end{array}\n", "$$" ] }, { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 3 }, "source": [ "## Glossary \n", "\n", "\n", "- **driver** A routine that calls the other routines to solve the\n", " problem.\n", "\n", "- **embedded Runge-Kutta methods**: Two Runge-Kutta\n", " methods that share the same stages. The difference between the solutions\n", " give an estimate of the local truncation error.\n", "\n", "- **explicit** In an explicit numerical scheme, the calculation of the solution at a given\n", " step or stage does not depend on the value of the solution at that step\n", " or on a later step or stage.\n", " \n", "- **fourth-order Runge-Kutta method** A popular fourth-order, four-stage, explicit Runge-Kutta\n", " method.\n", "\n", "- **implicit**: In an implicit numerical scheme, the\n", " calculation of the solution at a given step or stage does depend on the\n", " value of the solution at that step or on a later step or stage. Such\n", " methods are usually more expensive than implicit schemes but are better\n", " for handling stiff ODEs.\n", "\n", "- **midpoint method** : A two-stage,\n", " second-order Runge-Kutta method.\n", "\n", "- **stages**: The approximations\n", " to the derivative made in a Runge-Kutta method between the start and end\n", " of a step.\n", "\n", "- **tableau** The tableau for a Runge-Kutta method\n", " organizes the coefficients for the method in tabular form.\n", "\n" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "all", "encoding": "# -*- coding: utf-8 -*-", "formats": "ipynb,py:percent", "notebook_metadata_filter": "all,-language_info,-toc,-latex_envs" }, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.6" }, "latex_envs": { "LaTeX_envs_menu_present": true, "autoclose": false, "autocomplete": true, "bibliofile": "biblio.bib", "cite_by": "apalike", "current_citInitial": 1, "eqLabelWithNumbers": true, "eqNumInitial": 1, "hotkeys": { "equation": "meta-9" }, "labels_anchors": false, "latex_user_defs": false, "report_style_numbering": false, "user_envs_cfg": false }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": true, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "165px" }, "toc_section_display": "block", "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }