Permalink
Cannot retrieve contributors at this time
Name already in use
A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
CompMech01-Getting-started/notebooks/03_Numerical_error.ipynb
Go to fileThis commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
1139 lines (1139 sloc)
95 KB
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| { | |
| "cells": [ | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "###### Content modified under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2020 R.C. Cooper" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "# Numerical Error\n", | |
| "## Freefall Model Computational solution\n", | |
| "\n", | |
| "<img src=\"../images/freefall.png\" style=\"width: 200px;\"/> \n", | |
| "\n", | |
| "Here is our first computational mechanics model. \n", | |
| "\n", | |
| "An object falling is subject to the force of \n", | |
| "\n", | |
| "- gravity ($F_g$=mg) and \n", | |
| "- drag ($F_d=cv^2$)\n", | |
| "\n", | |
| "Acceleration of the object:\n", | |
| "\n", | |
| "$\\sum F=ma=F_g-F_d=mg - cv^2 = m\\frac{dv}{dt}$\n", | |
| "\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "### Define constants and analytical solution (meters-kilogram-sec)\n", | |
| "\n", | |
| "We define parameters for our problem as the acceleration due to gravity, g, drag coefficient, c, and mass of the object, m. Once we have defined these parameters, we have a single variable whose derivative $\\frac{dv}{dt}$ is equal to a function of itself $v$ i.e. $\\frac{dv}{dt} = f(v,~parameters)$. \n", | |
| "\n", | |
| "**parameters:**\n", | |
| "\n", | |
| "g=9.81 m/s$^2$, c=0.25 kg/m, m=60 kg\n", | |
| "\n", | |
| "**function:**\n", | |
| "\n", | |
| "$\\frac{dv}{dt} = g-\\frac{c}{m}v^2$\n", | |
| "\n", | |
| "We can solve for the analytical solution in this case. First, consider the speed of the falling object when acceleration is $\\frac{dv}{dt}=0$, this is called the terminal velocity, $v_{terminal}$. \n", | |
| "\n", | |
| "$v_{terminal}=\\sqrt{\\frac{mg}{c}}$\n", | |
| "\n", | |
| "Now, substitute this terminal velocity into the equation and integrate to get the analytical solution v(t):\n", | |
| "\n", | |
| "$v(t)=v_{terminal}\\tanh{\\left(\\frac{gt}{v_{terminal}}\\right)}$. \n", | |
| "\n", | |
| "## Exercise:\n", | |
| "\n", | |
| "Calculate the terminal velocity for the given parameters, g=9.81 m/s$^2$, c=0.25 kg/m, m=60 kg.\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 306, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "c=0.25 \n", | |
| "m=60\n", | |
| "g=9.81 \n", | |
| "\n", | |
| "\n", | |
| "def v_analytical(t,m,g,c):\n", | |
| " '''Analytical solution for the velocity of an object released from rest subject to \n", | |
| " the force of gravity and the force of drag with drag coefficient, c\n", | |
| " \n", | |
| " Arguments \n", | |
| " ---------\n", | |
| " t: time, the independent variable\n", | |
| " m: mass of the object\n", | |
| " g: acceleration due to gravity\n", | |
| " c: drag coefficient\n", | |
| " \n", | |
| " Returns\n", | |
| " -------\n", | |
| " v: the speed of the object at time t'''\n", | |
| " \n", | |
| " v_terminal=np.sqrt(m*g/c)\n", | |
| " v= v_terminal*np.tanh(g*t/v_terminal)\n", | |
| " return v" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Inside the curly brackets—the placeholders for the values we want to print—the `f` is for `float` and the `.4` is for four digits after the decimal dot. The colon here marks the beginning of the format specification (as there are options that can be passed before). There are so many tricks to Python's string formatter that you'll usually look up just what you need.\n", | |
| "Another useful resource for string formatting is the [Python String Format Cookbook](https://mkaz.blog/code/python-string-format-cookbook/). Check it out!" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "If we print these values using the string formatter, with a total length of `5` and only printing 2 decimal digits, we can display our solution in a human-readable way.\n", | |
| "\n", | |
| "```python\n", | |
| "{:5.2f}\n", | |
| "```\n", | |
| "where \n", | |
| "\n", | |
| "- `:5` prints something with whitespace that is 5 spaces total\n", | |
| "- `.2` prints 2 significant figures after the decimal\n", | |
| "- `f` tells `format` that the input is a floating point number to print\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 347, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "at time 0.00 s, speed is 0.00 m/s\n", | |
| "at time 2.00 s, speed is 18.62 m/s\n", | |
| "at time 4.00 s, speed is 32.46 m/s\n", | |
| "at time 6.00 s, speed is 40.64 m/s\n", | |
| "at time 8.00 s, speed is 44.85 m/s\n", | |
| "at time 10.00 s, speed is 46.85 m/s\n", | |
| "at time 12.00 s, speed is 47.77 m/s\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "for t in range(0,14,2):\n", | |
| " print('at time {:5.2f} s, speed is {:5.2f} m/s'.format(t,v_analytical(t,m,g,c)))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Analytical vs Computational Solution\n", | |
| "\n", | |
| "The analytical solution above gives us an exact function for $v(t)$. We can input any time, `t`, and calculate the speed, `v`.\n", | |
| "\n", | |
| "In many engineering problems, we cannot find or may not need an exact mathematical formula for our design process. It is always helpful to compare a computational solution to an analytical solution, because it will tell us if our computational solution is correct. Next, we will develop the **Euler approximation** to solve the same problem." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "# Define numerical method\n", | |
| "## Finite difference approximation\n", | |
| "\n", | |
| "Computational models do not solve for functions e.g. v(t), but rather functions at given points in time (or space). In the given freefall example, we can approximate the derivative of speed, $\\frac{dv}{dt}$, as a finite difference, $\\frac{\\Delta v}{\\Delta t}$ as such,\n", | |
| "\n", | |
| "\n", | |
| "$\\frac{v(t_{i+1})-v(t_{i})}{t_{i+1}-t_{i}}=g-\\frac{c}{m}v(t_{i})^2$.\n", | |
| "\n", | |
| "\n", | |
| "Then, we solve for $v(t_{i+1})$, which is the velocity at the next time step\n", | |
| "\n", | |
| "$v(t_{i+1})=v(t_{i})+\\left(g-\\frac{c}{m}v(t_{i})^2\\right)(t_{i+1}-t_{i})$\n", | |
| "\n", | |
| "or\n", | |
| "\n", | |
| "$v(t_{i+1})=v(t_{i})+\\frac{dv_{i}}{dt}(t_{i+1}-t_{i})$\n", | |
| "\n", | |
| "Now, we have function that describes velocity at the next timestep in terms of a current time step. This finite difference approximation is the basis for a number of computational solutions to ordinary and partial differential equations. \n", | |
| "\n", | |
| "Therefore, when we solve a computational problem we have to choose which points in time we want to know the velocity. To start, we define time from 0 to 12 seconds\n", | |
| "\n", | |
| "t=[0,2,4,6,8,10,12]" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 352, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import numpy as np\n", | |
| "#t=np.array([0,2,4,6,8,10,12])\n", | |
| "# or \n", | |
| "t=np.linspace(0,12,7)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "Now, we create a `for`-loop to solve for `v_numerical` at times 2, 4, 6, 8, 10, and 12 sec. We don't need to solve for `v_numerical` at time 0 seconds because this is the initial velocity of the object. In this example, the initial velocity is v(0)=0 m/s." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 353, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "array([ 0. , 19.62 , 36.03213 , 44.8328434 , 47.702978 ,\n", | |
| " 48.35986042, 48.49089292])" | |
| ] | |
| }, | |
| "execution_count": 353, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| } | |
| ], | |
| "source": [ | |
| "v_numerical=np.zeros(len(t));\n", | |
| "for i in range(1,len(t)):\n", | |
| " v_numerical[i]=v_numerical[i-1]+((g-c/m*v_numerical[i-1]**2))*2;\n", | |
| "\n", | |
| "v_numerical" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "Let's print the time, velocity (analytical) and velocity (numerical) to compare the results in a table. We'll use the `print` and `format` commands to look at the results. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 354, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "time (s)|vel analytical (m/s)|vel numerical (m/s)\n", | |
| "-----------------------------------------------\n", | |
| " 0.0 | 0.00 | 0.00\n", | |
| "\n", | |
| " 2.0 | 18.62 | 19.62\n", | |
| "\n", | |
| " 4.0 | 32.46 | 36.03\n", | |
| "\n", | |
| " 6.0 | 40.64 | 44.83\n", | |
| "\n", | |
| " 8.0 | 44.85 | 47.70\n", | |
| "\n", | |
| " 10.0 | 46.85 | 48.36\n", | |
| "\n", | |
| " 12.0 | 47.77 | 48.49\n", | |
| "\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "print('time (s)|vel analytical (m/s)|vel numerical (m/s)')\n", | |
| "print('-----------------------------------------------')\n", | |
| "for i in range(0,len(t)):\n", | |
| " print('{:7.1f} | {:18.2f} | {:15.2f}\\n'.format(t[i],v_analytical(t[i],m,g,c),v_numerical[i]));" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "## Compare solutions (plotting)\n", | |
| "\n", | |
| "We can compare solutions in a figure in a number of ways:\n", | |
| "\n", | |
| "1. plot the values, e.g. $v_{analytical}$ and $v_{numerical}$\n", | |
| "\n", | |
| "2. plot the difference between the values (the absolute error) e.g. $v_{numerical}-v_{analytical}$\n", | |
| "\n", | |
| "3. plot the ratio of the values e.g. $\\frac{v_{numerical}}{v_{analytical}}$ (useful in finding bugs, unit analysis, etc.)\n", | |
| "\n", | |
| "4. plot the ratio of the error to the best estimate (the relative error) e.g. $\\frac{v_{numerical}-v_{analytical}}{v_{analytical}}$\n", | |
| "\n", | |
| "Let's start with method (1) to compare our analytical and computational solutions.\n", | |
| "\n", | |
| "Import `pyplot` and update the default plotting parameters." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 358, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "import matplotlib.pyplot as plt\n", | |
| "\n", | |
| "plt.rcParams.update({'font.size': 22})\n", | |
| "plt.rcParams['lines.linewidth'] = 3" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 360, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0, 0.5, 'velocity (m/s)')" | |
| ] | |
| }, | |
| "execution_count": 360, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEaCAYAAAAsQ0GGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3hUxfrA8e9sGiSBJJTQQgcpkQDSQQXEQpduVwQEBcu9V0WRCyp6FeWqFxRQVATrT4xILyJIB2lq6L0TehJSSLLJzu+Ps9lNTzbZ3bT38zz77M6c2XMmbd/MnClKa40QQgjhbKairoAQQojSSQKMEEIIl5AAI4QQwiUkwAghhHAJCTBCCCFcwrOoK1BcVKlSRderV6+oqyGEECXK7t27r2qtq2Z3TAKMVb169di1a1dRV0MIIUoUpdTpnI5JF5kQQgiXkAAjhBDCJSTACCGEcAkJMEIIIVxCAowQQgiXkAAjhBBlVcQC+OhWeCPQeI5Y4NTTyzBlIYRwRMQCWDsFYs5BQAj0mAxhw1x3Pa3BkgraAjo102vrs7ZY89O/tmTMz5x37DfYMh1Sk4zrxJyFpc8br5309ShZrt/Qtm1bLfNghCiDLKmQHAdJsZBkfU6OzT594S84tdn4wE6jTBDUAHwrZQoCOt3rXD740+dnF0hw82d0QG345758F1dK7dZat83umLRghBCu5+z/+i0We1CwBYfYPNJxkHQj0/E4MMcX7mvTFrh+DK4X7jTFRsw5p51KAowQwrUiFhhdL+abRjrmLCx5Dq6fhDodMn7YZwgAObUm4oxHWaJMoDyMZ5NHutfZ5JtM9vImj3SvjXytPEjVJkyXIjBZkrNeKyDEadWWACOEcA1LKlw+CCvH24NLmpREWP9O0dQrO97+4FPB/uzjDz4VM6UrwJYZkBid9f1+wfDAN5k+zPMbENKCgClTQEiXr1S21b6ZnEpUQjJRCclEJ5itr81Ex1ufrceiEszExBnHY26a0Rr6mzYz1esLfFW6IONV3mhdOokEGCGEcyTHw/ndcOYPOLMNzu00WiSu4u2fNQB4V8iUtgaK3NLe/sYHf34E1M7YGgPjQ/m+/0CdjgX+UlItmpibZqLikolOSCIq3pw1aKRLpz0npVgKfM0lltvBDOM9F1BTXSOhfHX8e09x6oAFCTBCiIK5EQlntxsB5ex2iIzIePM7Lx5eULtjNgGggv2RU9rbz/gP393SPnxzuJ+kteamOZWoBDNR8ekDghEkohKSickUNKISzNxINFoV7lKhnCdBvt6c9u3DRN+BBPl6cX+rWnRvGuzU60iAEULkLa27K31AiT6T9/v8go0P4Yt7wWK253uVh34zXDu818nMqRZOX4vnmOkOjrX8hfPRiUaA2J5M9LqNtqCRXIhWhaO8PUwE+noR5Ottew7y8yLQ15sgXy8Cy1vz/axpX28Cy3vh6eGeKZASYIQQWSXHw7ldcPYPOLM9/91dVZsZN+5rdzSeg+ob9w/cPXekEG4mp3L8ShzHLqd7XInj1NV4Uiyua2ZULOdJkJ+3LQikBYTMQSN9MPH19kDlcH+mOJAAI4QoWHeXZ3mo1cYeUGq3g/JB2ZcNG1bsAkpMgpljV2JtQeSo9fl89M1CdVfl2apICxq+9rwAN7Yq3EkCjBBlTUG7u/yrQe0Oxs3s2h2hegvw9HZ9fQtBa82V2CRb8EgfTK7GJTl8vpoB5WgY7E+jYH/qVfajkp+3PZD4GS2P4t6qcCcJMEKUdgXt7gpuni6gdICgejkOly1qFovmXNRNW4vk6CWjW+vY5ThiE1McOpdJQd3KfjSs6k/jav40qmoElIbB/vj7yEemI+S7JURpc+OCEUjSAsrFvc7t7ipCySkWTl2Lz9IaOXElzuEhu96eJhpU8aORtUXSOLiC0TKp4ouPZxGMUCuFJMAIUZI5q7urRpgxbLiYSEhO4fjleI5ejs1wo/30tQRSHbzR7u/jaQsijYKNFknjav6EBPniYSqeLbLSQgKMEMVZ5tFXXV+BwDoOdncpCG5WLLu7ouKTbV1Zaa2R49Yb7Y6q4u+dqVvLaJFUq+gj90SKiAQYIYqrbNfwejbv93mWh5C29oAS0g7KB7q2rvlw9FIsW45dtd1wP34ljqtx2ayFlYdageXTdWvZWyaBvsV7wEFZJAFGiOLq10lZ1/DKTvrurjodoXrx6e46fiWO5RGRLI+I5PCl2Hy/z8OkqFfZN0PXVuPgCjSo6oevt3xslRTykxKiuNEads2FuIs5l2nzZLHr7kpz+lo8yyIiWRYRycHI3LvvynmZaFAl42itRsH+1K3sh7dn6ZsXUtZIgBGiOIm7Yixlf2RlzmUCakO//7mvTvlw9noCy/caLZW952OyLVPOy0T3JsG0rhNoG7FVK7A8JrnRXmpJgBGiuDjyKyweC/FX0mUqMuxo6OTl1AvjQvRNlkdEsmxvJH+fzWYJe4yhwN2bVKVPWE16NA3GT+aRlCny0xaiqJlvGvdbdn6eMb/D08Zs+fVTi80aXpduJBpBJeICe87kEFQ8TNx5S1X6htWgR7NgKpQrHveDhPtJgBGiKEVGwM+j4Ophe55/NRgwCxrdbaRbP1o0dbO6HJvIqn0XWfZ3JDtPX892nS5Pk+KOxlXoG1aTu5tXI6C8BBUhAUaIomGxwLaPYe1bGZexb9IH+s8AvypFVzfgWlwSK/ddZFnEBf44mX1Q8TApujSqQt8WNbg3tJoMExZZSIARwt1izsOip+HkRnuely/0fBdue6LIRoRFxSezev9FlkVEsu3EtWxnzJsUdGpYmb5hNbkvtDqV/CSoiJxJgBHCnfb/Akv/kXFf95qtYdAXUKWR26sTk2Bm9YGLLI+IZMuxq9nud6IUdKhfiT5hNel1a3Wq+Pu4vZ6iZJIAI4Q7JN6Ala/A39/b85QJbv8XdHvVrRMjbySa+e3AJZZFRLLp6BXMqdmv7dWuXhB9rUEluGI5t9VPlB4SYIRwtTN/wMKnIPq0PS+gDgz6DOp2dksV4pJSWHvQCCobDl8hOTX7lYdvqxNIn7Ca9GlRg+oBElRE4UiAEcJVUlNg4/uwcRrodB/oYQ9A72lQLsCll09ITmHtwcssj4jk98OXc1zOvmVIgNFSaVGdkCBfl9ZJlC0SYIRwhWvHYeFoOL/LnucTAH0/hBZDXHbZRHMqvx+6zLKISNYeukSiOfugcmutivRpUZO+YTWoXUmCinANhwKMUqoy0B1oDVQDAoEo4DKwB1ivtb7m7EoKUWJoDX9+a9xvMcfb8+veDgM/hcDaTr9kojmVjUeusCwikt8OXiIhOfvNxZrVqEjfsBr0blGD+lX8nF4PITLLM8AopTyBocBYoBPG2hXZjaPUgFZKbQVmAeFaa8f2Ks143XeACdbky1rr/+ZQ7mHgGSAM8AAOAV8Bs7XWjm1xJ0RhJFw3ltc/uNSeZ/KEu/4NnZ8Hk/N2SUxOsbDp6BWWR0Sy5sAlYpOy/1O7pZo/fVrUpE9YDRoF+zvt+kLkR64BRin1GPAOUBMjqFwCtgEHgOvADaAiUBloDnQEbge6AO8rpV7TWn/raKWUUu2A8RhBK8dJAUqpmRiBLxFYC5iBHsAnQA+l1FCt89orVggnOP47LHoGYiPteZUbw+DPjWHITmBOtbDl2FWWRUTy6/6L3Mhhr/kGVf3oG2Z0f91SrYJTri1EQeQYYJRSfwBtMYLKB8B8rfX+vE6olLoVGA48DMxXSj2rte6Y3woppXyAedbr7gAG5FBuMEZwuQjcqbU+as2vBvwODASeBabn99pCOMycaOw4uX1mxvy2I+De/4B34e5vaK3ZcuwayyIusGr/RaITzNmWq1fZl75hRkulafUKsoOjKBZya8HUBp4H5mits/+tzobWeh/wklJqAjAGeM3BOk3BaA31BwbnUi6t++yVtOBivf4lpdQzwHrgVaXUx9JVJlzi0gFjHbHL6f7v8q0C938CTXoV+vSnr8Xzys8RbD9xPdvjtSuVt92oD61ZUYKKKHZyCzANtdaOb4xtZQ1Knyilvszve5RSHYAXge+11kutrZTsyoUAbYBk4Kdsrr1BKXUeqIXRbbe1AF+CENmzWGDHZ7DmdUhNsuc3vhfunwn+wYU6fapFM3/rKaatPsxNc8Ye3lqB5ekTVoM+LWoQFhIgQUUUazkGmMIEl4KcRylVDpiPcW/nhTyKp3Vq78/l/DsxAkxrJMAIZ4m9CIvGwvG19jzPcnDv29BuVKHXETt2OY5Xfo5g9+koW56HSfFw+zoMvK0WrWsHSlARJUZxmgfzH6AJ8KDW+moeZetbn0/nUuZMprJCFM7BZcZukzfTdVlVb2GsIxbctFCnTkm18Pmmk3z02xGS002IbFq9AtOGtKRFiGsnZQrhCvkOMNY5MA2BE+kDgFKqFvAe0BI4BUzWWv/pSCWUUp2BfwCLtNY/5uMtaeMt43MpE2d9znEYjVJqNDAaoE6dOvm4rCiTkuNh1QTYMz9dpoLOzxlDkD0Lt/jj4YuxvBz+NxHn7FsNe5oUz97ViLHdGsne9KLEcqQFMwH4J0aX01WwjfjaDNTBGE4cCtyulArTWp/Nz0mVUuUx5q3cwBgVlq+3WZ+zX6Uvn7TWc4A5AG3bti3UuUQpdX43/PwUXD9uz6tYy5g0Wf/OQp3anGph9vrjfLzuaIYFJ1vUCuD9IWE0q1GxUOcXoqg5EmC6Y7ReItLlPQjUBdZhzJfpjzHy7FnglXye9x3gFmCE1joyr8JWsdbn3GaOpR2LzaWMENmzpMLmD43tii3p5puEDoS+H0H5oEKdft/5GF4Oj+Bg5A1bnreHiX/c05jRdzTA00NaLaLkcyTA1AL+ypTXB6MV8ZTW+iSwTinVF+hJ/gPMQMACPKGUeiLTsbSO7Wes5z2mtR6F0RUHRnDLSdqaHKdyKSNEVlGn4ZcxcGabPc+7grFAZcsHC3UjPykllY/XHmP2huMZNvRqXSeQaUPCaBQsEyNF6eFIgAnC2jWWTifgsDW4pPkTYza9I0xA11yON7A+AtNdAyBUKVU+h5Fk7TKVFSJvEQtg+YuQZG9ZULsDDJoDQfUKdeo/z0QxPjyCo5fjbHk+niZevq8JT3apj4dJRoeJ0sWRAHMTsG0UrpSqg9GqyTzPJRnI9z6qWut6OR1TSs0DniDTWmRa67NKqT3AbRjrpH2d6X1dgRCMWf7bECIvN6ONwLIv3J6nPIzNwG7/F3gUfMBlojmVD9cc4YtNJ0i/YWT7+pV4b3CYLDwpSi1H/moOYNzAr2IdRfYIRvfYxkzlamMs8+Jq72JMsnxPKbVVa30MQCkVjLHYJsBUmcUv8nRqM/zyNMSkG5cSVB8GfwEhbQt16p2nrjM+PIKTV+0DHn29PXi1V1Me7VAXk7RaRCnmSID5GuODe5e19dAH4wb64rQC1smStwEbnFnJ7Gitw5VSszFWUt6rlPoN+2KXFYFFGIteCpG9lGRY/w5s/h8ZBiS2fgx6TgWfgq8+nJCcwvurDjN/2yl0ulPf3qgK7w5qIXuwiDLBkQAzB2PZlccxhiXHAiO11uk6q+kP+OKGAAOgtR6rlNoMjMO4h5O2XP9cZLl+kZsrR2DhKIj8255XPgj6zYDm/Qt16q3HrvLKwgjOXrffGqzg48nEPs14oF1tmYkvygyltWPTP6z3XoKBQ1rruEzHWmGM7NqutXZHN5nTtG3bVu/atSvvgqJk0xp2zYXVEyEl3diQBt1hwGyoWKPAp45NNPPuykN8/8eZDPndm1TlnUEtqBFQvsDnFqK4Ukrt1lpn25ec23L9jwLLtNbR6fO11mewL8NCpmN/kXUosxDFQ9wVY6mXIyvteR7ecPeb0OFpMBV87sn6w5d5beFeLsQk2vICynvxer/mDGxdS1otokzKrYvsa8CslPod+AVYrLW+6J5qCeFkR36FxWMh/oo9L7g5DPocqt9a4NPGJJh5a/kBwnefy5B/b/NqvD3gVoIrlivwuYUo6XILMC8Cg4C7gXuBmUqpbRjB5pdMc1+EKJ7MN+HXSbDz84z5HZ6Bu98Ar4IHgDUHLjHxl71cjrUv2V/Jz5s3+4fSN6yGtFpEmZfbcv0fAR9Zh/0OxAg23TC2Q56mlNoL/IwRbPa5oa5C5E/EAmOXyZhzYPLIuNSLfzXjXksjR+cC212PT+aNJftZ8veFDPl9w2rwZv9QKvsXbvFLIUoLh27yK6UqYowUG4TRqvHFGN95AiPYLNJab3dBPV1ObvKXEhELYOnzRssls6Z9jVFifpULfPrlEZFMXryPa/HJtrwq/j68PeBWet5avcDnFaKkyu0mv8OjyNKdtDzQCyPY9MZYxkUDkRjdaIuA30vKUGEJMKXER7dmnDCZpnwQjD9Z4HXErsQmMXnxPlbuy3gbctBttZjctzmBvvlevEKIUqVAo8jyYl3/ayGwUCnliTHBcRBGC2ccxtL7bwBvFfQaQjgs5lz2+TejCxRctNYs+us8by49QHSC2ZZfvWI53h3Ugu5NC7c9shClmVN2tNRapwCrgdVKqaeB2zHu2+R3+X0hCi/2ohFEsmuVB4Q4fLqLMYlM/GUvaw9dzpD/UPvaTOjdjIrlvApaUyHKBKdvmayNPrdN1ocQ7pGaAuEjIbseWa/y0GNyvk+ltWbBrrO8vewgsUn2AQK1Asvz3uAwbm9cJZd3CyHSFCjAKKVCgJpAjmM8tdaZF8EUwnXWvwOnN9vTflUh/qrRcukxGcKG5es056ISmLBwL5uOZtyZ4vFOdXmlZ1P8fJz+P5kQpZZDfy1KqYeAN4GGeRTVjp5biAI7ugY2fWBPd3sNuuV3vzuDxaL5bscZpq44SHxyqi2/bmVf3h8cRocGBR95JkRZle8goJR6BGN2vwKuAyeBuFzfJISrxZyDhU/Z0w26w50vOXSK09fieeXnCLafuG7LUwpGdqnPi/c2oby3h7NqK0SZ4kgrI+1fwnHAHK11am6FhXC5lGT4aTjcjDLSFWoYS7+Y8hcQUi2aeVtPMW31IRLN9ns3Dav6MW1oS26rE+SCSgtRdjgSYBoDm7XWs11VGSEcsvZNOLfTeK08YMhX4F81X289djmO8eF/s+eMfS1XD5NizJ0NeL5HY8p5SatFiMJyJMBcA867qiJCOOTgUtiWbj+5u1+Hup3yfFtKqoXPN53ko9+OkJxib7U0rV6BaUNa0iIkwBW1FaJMciTArAa6KqWULuj0fyGc4fpJWDTOnr6lF3R6Ls+3Hbp4g/HhEUSci7HleZoUz97ViLHdGuHtWfDl+oUQWTkSYF4HdgL/VUq9Yp1cKYR7mRPhpycgyRokAurAgFm57uViTrUw6/fjfPL7Ucyp9v+Nbq1VkWlDWtKsRkVX11qIMinfAUZrfU4p1QVYCgxUSq0DzgHZrTWmtdayRIxwvl8n2rc5NnnB0HngWynH4vvOx/ByeAQHI+07e3t7mHjh7saMubMBnh7SahHCVRwZpmzC2COmCWACRmRTTGMMY9bIGmTC2faGw84v7On7/gMhbbItmpSSysdrjzF7w3FSLfZWS+s6gUwbEkaj4Aqurq0QZZ4jXWQTgGcAM0Yr5hgyD0a4y5UjsOR5e7r5/dB+dLZFz0UlMGLeTo5csv96+niaePm+JjzZpT4eJtkITAh3cCTAjADigS5a6wgX1UeIrJITjPsu5ngjXakB9P8429WRr8cn8/iXOzhxNd6W175eJd4bEkb9Kn7uqrEQAscCTA2M/V0kuAj3WvEyXD5gvPbwgaHzoVzW4cQJySmMmLfTFly8PUxM7NOMxzrWxSStFiHczpEAcx7IZptAIVzoz2/hr2/t6d7vQ42wLMXMqRbGfbeHv84aEyeVgv892IreLWq4q6ZCiEwcGULzf0A3pZS/qyojRAaX9sPydOuKhT0Itz2RpZjWmtcW7uX3w1dseW/0C5XgIkQRcyTAvA0cAJYppW5xUX2EMCTFwoInIMXaaK7aFPp+mO19lw9+PcJPu+07WY7r3pAnOtdzU0WFEDlxpItsJUZAuh3Yr5Q6Te7zYHo4oX6iLNIalr4A144aaS9f476Ld9ab9PO3nuKT34/Z0kPahPDSvU3cVVMhRC4cCTDd0r32ABpYH9mRpWREwe2aC/t+tqf7/g+Cm2YptmJvJG8s3W9Ld29SlXcHtUBl08oRQrifIwGmu8tqIUSaC3/Bqlft6duegJYPZCm27fg1/vF/f5G2Kl6r2oHMfOQ2vGRmvhDFhiNLxWxwZUWE4Ga0Md8lNdlIV2sBvd7LUuxg5A1Gf72L5FSjd7ZBFT/mDm+Hr7dsoipEcSL/7oniQWtYPA6iThlp7wowbD54lc9Q7FxUAsO/2kFskrHWatUKPswf0Z5Kft5urrAQIi8SYETxsH0WHFpmT9//MVRumKFIVHwyj8/dwaUbSQBU8PFk/pPtqV3J1501FULkU44BRim1VSl1Z2FOrpTqqpTaUphziDLg7A5YM9mebj8GQgdmKHIzOZUR83dy4op9lv5nj7eheU1Zal+I4iq3FkwD4Hel1O9KqUeUUuVzKWujlCqvlHpMKbUeWAfUd0I9RWkVfw1+Gg4W6/ZCNW+DezMuxJ2SauHZ7/fw5xn7LP0PH2hJ54ZV3FxZIYQjcrsr2hhjk7HngDuB2UqpbcA24CDGFso3gIpAZaA50Mn68MVYdfkDZNl+kROLBX4ZAzesO3GXCzT2d/H0sRXRWjPxl32sPXTZlvd63+b0Davp5soKIRyVY4DRWscCLymlPgaeBZ4E7gHuzuEtaZMPrgIfA7O11medWFdR2mz5CI6tsacHfgZBdTMU+WjNEX7cZf81eqZbQ4Z3kUaxECVBnuM6tdangZeVUv8G7sCYcNkKCAYCgGjgMrAH+B3YorU2u6rCopQ4tRnWvW1Pd3kBmvTMUOSb7aeZsc4+S3/QbbUYf5/M0heipHBkHkwS8Jv1IUTBxV2G8BGgrasM1ekEd03KUGTVvkgmL95nS3drUpX3BofJLH0hShAZpizcy5IKP4+EuEtG2rcyDJkLHl62In+cuMbz6WbptwwJYJbM0heixJG/WOFeG96DkxutCQWDPoeK9hv2hy7eYNTXu0hOMVo39WWWvhAllgQY4T7H1sKG9+3pruOhkX3R7fPRNxk+dyexifZZ+l+PaE9lf5/MZxJClAASYIR7xJyHhU9hW2i7/p3Q9RXb4eiEZJ6Yu4OLNxIB8Pfx5Kvh7WSWvhAlmAQY4XqpZuOmfsI1I+1fDQZ/CSYPwJilP3L+Lo5djgPAy0Mx57E23ForoKhqLIRwAgkwwvXWToGz243XymTc1PcPBoxZ+s/98Ce7T0fZin84rBWdG8ksfSFKOgkwwrUOrYCtM+zpuyZBvdsBY5b+pMX7+O3gJdvhyX2b06+lzNIXojSQACNcJ+o0LHranm58L3T5hy35v9+O8sMO+yz9MV0bMOJ2maUvRGmR7wCjlHpaKZV1U3QhspOSZCximRhjpCuGGEvBmIxfue/+OM30tUdtxQe1rsUr92XdFlkIUXI50oKZBZxTSk1XSskngcjdr5Pgwh7jtcnTWMTStxIAq/dfZNIi+yz9O2+pyntDwjCZZJa+EKWJIwFmEeCHsbryfqXUb0qpgUop6WYTGe3/BXZ8Zk/f8xbUbgfAzlPXee6HP7FYRyuHhQQwW2bpC1Eq5fuvWms9CKgHvA1cAu4CwoHTSql/K6WquaSGomS5dhwWP2dPN+0LHZ8B4MilWEbO22mbpV+vsi9zh7fDz0dm6QtRGjn0b6PW+oLWejJQB3gI2AzUAt4EziilflBK3eH8aooSwXwTFjwBybFGOqge3D8TlOJC9E2emLuDG9ZZ+lX8vfl6RAeqyCx9IUqtAvVLaK1TtNY/aq27Ai2Az4BEYBiwXin1t1JqdH53wRSlxMpX4NJe47WHNwydD+UDiUkw88TcHUTGGLP0/bw9mPdke+pUlln6QpRmhe741lrvB6YC32BsOqYwgs5sjO6zMYW9higB/v4/2DPfnu45FWq2ItGcyqivd3I03Sz9zx5rK7P0hSgDChVglFI9lVJLgOPAM0ACMAcYCiwBKgGzlFL/LGxFRTF2+SAsS/cjvnUItB1hm6W/85R9lv5/h7bk9sYyS1+IssDhu6tKqUrACOBpoD5Gi+UkxjDmL7XW0daiPyul2gDrgHHAR06psShekuKM+y7mBCNd5RboNx0NTFq8nzUH7LP0/92nGfe3qlU09RRCuF2+A4xSqj0wFuM+iw9GYPkN+BhYpnXa9lB2WuvdSqkVwBDnVFcUK1obLZerh420Z3njvouPPzN+O8oPO87Yio6+swGj7mhQRBUVQhQFR1ow1tUKiQO+Aj7RWh/Mx/viAQ9HKyZKgD3zYe8Ce7rvh1CtOT/sOMNHvx2xZQ9oVZNXe8rcXCHKGkfuwRwH/gmEaK3H5TO4oLUepbWWWXSlTWQErBhvT7d+FFo9zJoDl5j4y15b9h2Nq/D+kJYyS1+IMijfLRitdWNXVkSUIIkx8NMTkJpkpINDodc0dp26zrPf77HN0m9RK4DZj7bB21P+vxCiLHJksct1SqmX81HuJaXUusJVSxRbWsOS5+D6CSPt7Q/D5nM0ytg0LMk6S7+udZa+v8zSF6LMcuSvvxtwKh/lmgBdC1IZUQLsmAMHFtvT/WcQ6RXCE59vJeamGUibpd+eqhVklr4QZZkr+i58gFQXnFcUtXO7YfVEe7rdKGIa9Gf43J1cSDdL/6vh7albWXZ2EKKsc2qAsa6s3Aa46szzimIg4bqxv4vFaKVQoxWJd73FU1/v4vAlY+0xT5Pi08fa0CJEZukLIfLoIsvmXkrPXO6veAKNgGrAghzKiJLIYoFFz0CMdV6LTwCpQ+bxj/CD7Dh13Vbsv0NbckfjqkVUSSFEcZPXPZhu6V5roLr1kZs/gVcKUSdR3GydAUdW2ZJ6wCwmb4xj1f6LtryJvZsxoLXM0hdC2OUVYLpbnxXGki+rgPdyKJsMnNdan8nheI6UUl7AnUBvoAtQF6gMXAG2YUzqXJ/L+x/GWAstDGNS5yGMyaCztdYWR+sj0jm9FdZOsac7PcsnF5rw3R/2iZSjbq/PU3fKLH0hREa5Bhit9Ya016zjL+oAACAASURBVEqpDcD69HlO1BVYY319EdiNsQJAc2AwMFgp9ZZ1L5oMlFIzMZawSQTWAmagB/AJ0EMpNVRrLYMOCiLuCoSPgLRvX0h7fgwYwQeLDtmK3N+qJq/1blZEFRRCFGeOTLTsnnepArMAPwPTtdab0h9QSj0AfAdMUkr9rrX+Pd2xwRjB5SJwp9b6qDW/GvA7MBB4FpjuwrqXTpZUWDgKYiONdPlKbGo1jQk/24PL7Y2qME1m6QshclAsplhrrddprYdkDi7WYz8C86zJRzMdnmB9fiUtuFjfcwmjywzgVevoNuGIjf+FE+ttyaO3f8hTiyNts/RvrVWRTx+TWfpCiJzl2IJRSqV1R32itb6eLp0fWmv9VuGqlsGf1ueQtAylVAjGkOhk4KdsKrBBKXUeY0vnjsBWJ9andDuxHta/a0teb/M8Q9f6kWg2hijXqeTLV8Pbyyx9IUSucvuEeANj5Nj/AdfTpXPrD0k7rgFnBpi0ddAi0+W1tj7v11rfzOF9OzECTGskwOTPjUj4eRTGjxCSQjozYN+dRCckA1DZz5v5MktfCJEPuQWYKRifMlczpd1KKVUdGG5N/pzuUH3r8+lc3p42oq1+LmUEQMQCWPsmxJyzZVn8gnnyxhjOxBjBxdfbg6+ebEf9KjJLXwiRtxwDjNb6jdzS7qCU8gS+BQKAtVrrpekO+1uf43M5RZz1uUIO5x8NjAaoU6dO4SpbkkUsgKXPg9neENTAL3Rn62UvwJilP/vRNoSFBBZRJYUQJU1xv0P7KcaQ47NkvcGf1lVX4FaV1nqO1rqt1rpt1apleAb62ikZggsY39wOsWtt6feHhNH1ljL8PRJCOKzYBhil1HRgJMYQ5B5a64uZisRan/3JWdqx2FzKiHTdYunVVNcAmNCrKYNuC8m2jBBC5MSR/WDGKaVSlVJ9cynT11pmTGEqpZT6AHgeYyZ/j/RDkNM5ZX2um8upamcqK7Ljk20PIhd0ZUZ0qc9omaUvhCgAR1owA4HLwPJcyqzACAqDClohpdT7wL+Aa8A9WusDORRNG7ocqpQqn0OZdpnKisyuHIakuCzZCdqb32qO4d99mqGUTKQUQjjOkQDTFNintc7xnod13a+9QIHWDlFKTQVeBqIwgsvfuVzrLLAH8AaGZnOurhjzZi5irGcmMrNYYMnzGAspQLL2xKLhnKUKX1X6Jw+NelFm6QshCsyRmXJVgfX5KHcZuMPRiiil3sJYhTkaI7jkp9XxLsYky/eUUlu11ses5woGZlnLTJUFL3OwZx6c3Q5ACh70S36bw7oOzWtU5McxHfHx9Cja+gkhSjRHAkw0kJ+xvCHYhwfni1KqP/Bva/IY8FwO3TKHtNZT0xJa63Cl1GyMZWH2KqV+w77YZUVgEcailyKzG5Gw5nVbcnZKPw7rOpT38uDTR9tQoZxXEVZOCFEaOBJg9mCsTtw4h5vuKKUaA50AR1dcrpTudVvrIzsbgKnpM7TWY5VSm4FxGKsypy3XPxdZrj9nK8dD0g0ATurqfJIyAIBXezWlTmXfoqyZEKKUcCTAfAXcByxWSg3SWh9Kf1Ap1QRYiPEB/5UjldBaz8O+oKXDtNbfA98X9P1lzqHlcHCJLfmaeSRJeNOxQSUe65jboDwhhMg/R5brX6CUegToh9EdtQ2jpQDQBOiMEVyWWz/wRXGUeAOWv2RL/pjSjW2WUHy9PWTpfSGEUzm6HO4Q4L/A08Dt1kcaM8aN9ZedUzXhEuvegtgLAFzVFXkn5WHAmExZu5J0jQkhnMehAKO1NgMvKKXeBu7CPsnxNLBOa33FyfUTznR2J+z43JacYn6cGPzp1KAyj3SQrrHcpKSkcP36dWJiYkhJSSnq6gjhdJ6engQEBFCpUiU8PZ2zFUeBzmINJD86pQbCPVKSjQUtrUu3/Z7akiWWTvh6e/D+kDDpGsuFxWLh7Nmz+Pj4UKdOHby9vWXyqShVtNYkJydz7do1zp49S926dTGZCr+SWLFdi0w42dYZcNlYFCFB+zApZQSgmNC7mXSN5SEqKgpPT09q1KiBj4+PBBdR6iil8PHxoUaNGnh6ehIVFeWU8zocYJRSzZVSnymlDiul4qyPw0qpT5VSoU6plXCua8dhw/u25AcpQzinq9K5YWUeaV+GtynIp7i4OAIDAyWwiFJPKUVgYCDx8bntgpJ/DgUYpdRIjPkwozB2mfS1Phpj7Kuy21pGFBdaw9IXIDUJgAhLfeal9sTP24P3BkvXWH4kJibi6yutPFE2+Pr6cvNmTpsEO8aR1ZQ7AJ9h3Lf5CeiJEVhuwZgf86P12KfWsqI4+Os7OLUJgBRtYoL5KVLxkK4xB1gsFqf0RwtREphMJiwW58xPd+Qm/0sY+1A9pLVekOnYMWCNUmohRqB5ERjmlBqKgou7DKsn2pJfpPZmv67H7Y2q8EgH6RpzhHSPibLCmb/rjvxbdjuwM5vgYqO1/gnYQQEWuxQusGoCJEYDcMZSlekpg/D38WTq4BbygSmEcDlHAkwljJZKXo6RcW0xURSOroF94bbkxJSR3KQcr/VuRkiQdI0JIVzPkQBzHWiUj3INrWVFUUmKg2X/siUXpt7OJksYdzSuwkPta+fyRiGEcB5HAsxWoJ1SKsfdKpVSA4AOwJbCVkwUwvp3IeYMAFHan7fNj1q7xsKka0yUSfXq1UMpxalTp9xyPaVUkf2tnTp1CqUU9erVK5Lrp+dIgPkAYxr4j0qpr5VSPZRSDZRS9a2v5wMLMLZH/MAVlRX5cOFP2D7Llnzb/CjXqcjEPs2oFZjTztJCiPwaPnw4SinmzZtX1FUp9hxZTXmrUuo5YDrwiPWRngJSgOe01rJFcVFITTG2QLZugbM5NZSfLXdwR+MqPNhOusaEcJeDBw8WdRWKBUcXu5ytlNoCvADcCdSyHjqPsRnYDK11hHOrKPJt+yy4aHz7E7UXE1NG4u/jJV1jQrhZ06ZNi7oKxYLDs8e01hFa65Fa68Zaa1/ro7HWepQElyIUdQp+f8eWnJ4ymNO6OpP6SteYcL0//viDl19+mbZt21KtWjW8vb2pWbMmQ4YMYfv27VnKv/HGGyileOONN7h06RJjxowhJCQEHx8f6tevz6uvvkpiYmKW98XGxjJnzhwGDBhAo0aN8PX1xd/fn9atW/Of//wn3zPQtdY0btwYpVS29UszaNAglFLMmjXLdm9j/vz5ADz55JO2ey2Zu8xyuwdjNpuZM2cO3bt3p1KlSrZFVPv27ct3332Xoezp06d599136d69O7Vr18bHx4dKlSrRvXt3vv+++G+75Zw1mUXR0hqW/RNSjD+ug5Y6fJ7am663VGVYW+kaE643ceJE1q9fT2hoKO3bt8fHx4fDhw/z888/s2jRIn744QeGDh2a5X1nz56lTZs2aK3p3LkzN27cYPPmzbz33nscOHCAJUuWZCj/999/M2bMGIKDg2nSpAlt27bl2rVr/PHHH/z73/9myZIlbNiwgXLlyuVaX6UU48aN45///CezZs2iY8eOWcqcP3+epUuXUqFCBR577DGSkpJ44okn2Lx5M8ePH6dLly40amQfWJv+dU6ioqLo06cP27Ztw8fHhy5duhAcHMyFCxfYsmUL+/bt45FH7HcfvvnmGyZNmkTDhg1p2rQpXbp04dy5c2zatIn169fzxx9/MH369DyvW1QkwJQGe3+C4+sAsGjFq+ZRlPcpJxMq3aTeq8uLugoFdmpqH6ec56WXXuK7776jWrVqGfKXLl3K4MGDefrpp+nTp0+WNd3mzp3LqFGjmDlzJt7e3oBx/6J9+/YsXbqULVu20KVLF1v5evXqsXbtWrp165Zh+Z7o6GgeeughVq1axfTp03nllVfyrPOTTz7JpEmTWLBgAR9++CFVqlTJcPyzzz4jJSWFxx9/nAoVKlChQgXmzZvH8OHDOX78OKNGjWL48OEOfZ+GDx/Otm3b6NSpE+Hh4dSsWdN2LDExkd9//z1D+fvuu4+BAwcSGppxHeGjR4/So0cPZsyYwcMPP0yHDsVzda4cu8iUUnML8fjSnV9EmRZ/DVa9akvOS72Pv3UjJvVtTo0A6RoT7tGzZ88swQWgX79+DB06lOvXr2f58ASoXbs2M2bMsAUXgGbNmvHYY48BsHbt2gzlQ0JCuOuuu7KsDRcYGMiMGTMACA8PJz8CAgJsLZO5c+dmOGY2m/n8c2NzvrFjx+brfHn566+/WLJkCf7+/ixevDhDcAEoV64cvXr1ypDXrl27LMEFoHHjxkyaNAnI/9dbFHJrwQwvxHk1IKsqu8Ov/4aEawCc15X5IGUo3ZpUZWjbkCKumChrrl69yrJly9i3bx/R0dG2nT/37dsHwJEjR+jTJ2OL6a677qJ8+az/CKXdJL9w4UKWY1prtmzZwsaNGzl37hw3b95Ea43W2nad/Hr22WeZPXs2n376KS+99JItcC1cuJCLFy/SrVs3mjdvnu/z5WbVqlUA3H///VStWjXf70tMTGT16tXs3LmTK1eukJRkrIweGRkJOPb1ultuAeZJt9VCFMyJ9fC3/UbfJPOTmMpV4N1B0jXmTs7qZirJPvvsM/71r3+RkJCQY5kbN25kyatTJ/tFVytWrAiQ5Ub/pUuXGDRoEFu3bnXoOjlp3rw5d999N7/99hurVq2id+/eAMyaZcwlGzduXL7PlZfTp08Djo0w27ZtG8OGDePcuXM5lnHk63W3HLvItNbzC/Nw5xdRJplvwtJ/2JLLUjuyznKbdI0Jt9u1axfPPPMMZrOZadOmcejQIeLi4rBYLGitmTBhAoCthZGeo9sgjBo1iq1bt9KlSxfWrFnD5cuXSU5ORmtt+8/eUc899xxgDyr79+9n48aN1KxZkwEDBhTonM6QkJDAwIEDOXfuHCNHjmTXrl1ER0eTmpqK1prVq1cD2X9fiwu5yV9SbXgPok4CcEP78qb5cbo3qcrQNtI1JtwrPDwcrTXPP/88L730Upbjx47lZ43cvMXHx7NixQo8PDxYtmwZgYGBTrlO3759qV+/PitXruTUqVPMnDkTgNGjR+Pp6byPyLp16wJw+PDhfJXfuHEjly5dok2bNnzxxRdZjjvr++pKBdpFSSkVoJS6Wyn1kFKqs7MrJfJwcR9smWFLvpPyMInlqvDuIJlQKdzv+nVjbdvatbMOib9y5Qpr1qxxynViYmKwWCxUqFAhS3ABsswhyS+TycTYsWOxWCxMmzaNb7/9Fk9PT0aPHp1t+bQBCWn3mPLrvvvuA2Dx4sVcvXo1z/K5fV+BEjEPxtEtkwOUUnOBy8Bq4FuM7ZPTjo9VSl1QSmUdVC6cw5IKS54DnQrAH5am/Jjajdf7hVI9IPex/0K4Qto9ha+//pq4uDhbfmxsLCNGjCA6Otop16lWrRpBQUFER0dn+XBdtWoVH374YYHPPXLkSHx9fZk1axaxsbEMHDiQGjVqZFu2Vi1jARNHl4Np3bo1/fr1s50/7SZ9msTERFauXGlLp31f161bx6FDh2z5FouFKVOmsGVL8V9T2JEtk/2A9Rijy6KAlRjrj6W3CqgOFF3HZWm343O4sAeAJO3Ja+aRdG9ancG31crjjUK4xpNPPknt2rXZs2cPDRo0YNCgQQwcOJB69eqxa9cuRowY4ZTreHh4MHGisUPrI488QufOnW1zQHr16sW//vWvPM6Qs6CgIB599FFbOreb+/fffz8mk4n//e9/3HfffYwcOdJ2bygv8+bNo127dmzevJkGDRpwzz338PDDD9OtWzdq1KjBM888Yyt722230a9fP27cuEGrVq3o1asXDz74II0bN+att95i/PjxBf563cWRFsxLQEuMVksDrXXfzAW01ieAI8BdzqmeyCD6LKydYkvOTBnAFZ+6MmpMFKmgoCB27drF6NGj8ff3Z/ny5ezatYtBgwaxZ8+eHLt4CuLFF18kPDycjh07sn//fpYtW4aHhwfffvst//nPfwp17nvuuQeA0NBQunbtmmO5Vq1a8eOPP9KuXTu2bt3K3Llz+fLLL/M1XLhSpUps2rSJjz/+mNtuu40dO3awcOFCTp48yR133MHUqVMzlA8PD2fq1Kk0atSI9evXs3btWkJDQ9m8eXOWOTPFkcrvCASl1D4gEGiotU6y5lmAeVrrEenK/QqEaq1L1L/Ubdu21bt27SrqauRMa/j+AThqjBw5YqlFn+R3mTq0DYPlxr5LHTx4kGbNmhV1NYSLDRw4kEWLFjFr1qwMLYmyyJHfeaXUbq112+yOOdKCaQDsTAsuubgKVHbgvCI/DiyyBReACeZR3Nm0JoOka0yIQtu9ezdLliyhcuXKPP7440VdnVLDkTF4ZiA/d5FDgLg8S4n8uxkFK+z9rd+m9OCoTyhrpGtMiEIZNWoUcXFxrFixwnbz3M/Pr6irVWo4EmAOA62VUuW01lnX0QaUUkEY92n2OKNywmrN6xB/GYBLOpD3Uh7izUGhVKsoo8aEKIwvv/wSk8lE3bp1mTx5stPWHRMGRwJMODDV+vhHDmXeAfwxtk4WznBqC+yxL4ww2TycDs3qMbC1dI0JUVjFeRZ8aeBIgPkEeAJ4TinVFlhoza+nlHoGGAp0BfYCspqyM5gTYekLtuTq1LZs9+nCmoHSNSaEKP7yHWC01glKqXuBn4DOQCfroa7WhwJ2AwO01snOrmiZtPlDuHYUgFhdntfNT/DmoFCCpWtMCFECOLTQjtb6PNBZKdUT6I0xsswDOIsx8XKRljanc1w+iN70oW0m6/spD9CieXPub1Uz17cJIURxUaCV3LTWqzBm7QtXsFhg6QsoixmA3ZbGLPPuxeqBt0rXmBCixHBkqRiZaeYuu7+Cs38AYNYeTDCP4o37WxBcQbrGhBAlhyMTLfcppbYrpZ5WSmVdylQ4x40L6N9etyU/Te1HvWZt6d9SusaEECWLIwHmCtAemAlEKqX+TynVU0mfjXOtHI9KigXghKU633gN5W3pGhNClECOBJiaQF/gZ2t6GLAcOKeUmipdaE5wcBkcXGpLvpYyion3t5auMSFEiZTvAKO1tmitV2ithwE1gHHATuvr8RhdaH9IF1oBJd7AssK+G+CPKd0IaNZdusaEECVWgXa01FpHa61na607Ak2B94ALQDuMLrQLzqtiGbF2CqZYYwOiK7ois7we5+0BMqFSCFFyFSjApKe1PqK1noAxJ2YGxoRLn8Ket0w5uwO9077n9hTz47w0oBNVK8i3UYiSrl69eiilOHXqVJFcXylVZP+oFmgeTHpKqVCMXS4fAapZs28W9rxlRkoyqYufwwNjfuq61FakNBtI37Dst2sVQoiSokABRilVCXgYY22y27BvnbwVmAf86IzKlQlbp+Nx1dhvO0H78IHXGObLWmNClBpr167FbDZTq1bZW6A23wFGKeWBsTzMExijybwwAss54BuMnS2PuqKSpdbVY1jWv2/rp/wgZSjPDO5GFX/pGhOitGjYsGFRV6HIOHIP5jywCBgEWDBaKT2BulrriRJcHKQ1KYufw2Qx1gX929KAi82eoG+YjBoTQMQC+OhWeCPQeI4ovjtgpO/j//HHH+nUqRP+/v5UqFCBHj16sHnz5izvOXXqFEop6tWrl6/z5pQ/b9482rZti5+fH9WrV2fkyJFcuXIFgMTERF5//XVuueUWypUrR506dZg4cSJmsznHa65evZr+/ftTrVo1vL29qVGjBg899BB79+7N9WtISUnhv//9Ly1btsTPz4/AQPtA2tzuwWitWbBgAb169SI4OBhvb29q1apFjx49+OSTTzKUvXLlCtOnT6dnz57Ur1+fcuXKERAQQMeOHZk5cyapqak5fl1FRmudrwdGUNkGjAEC8vu+kvJo06aNdqvdX2v9ekWtX6+ozZMD9YNvfqqvxCa6tw4iXw4cOODeC/79o9ZvV7P9fujXKxrpv390bz3yCdCAnjRpkjaZTPrOO+/Uw4YN002bNtWA9vb21lu3bs3wnpMnT2pA161bN8/z5pQ/fvx47e3tre+55x49aNAgXb16dQ3osLAwHRsbqzt37qyDgoL0gAEDdK9evbSvr68G9FNPPZXt9Z5//nkNaE9PT92pUyc9dOhQ3bp1aw3ocuXK6eXLl2f7NdSpU0f3799fe3t767vvvls/+OCDunPnzrZydevW1YA+efJkhvcnJSXp/v37a0B7eHjoLl266IceekjfddddOjg4OMvX/s0332hAh4SE6G7duukHHnhAd+vWTfv4+GhA33///dpiseT7+5gbR37ngV06h89VR+7BNNNaHy5kPBMAcZcxr5qIlzX5RWpvHh3cX7rGSqo3Alx/DfNNWPiU8XCmN2KcdqqZM2eyY8cO2rRpA4DFYuHpp5/m888/Z/LkyaxZs8Zp1wKYP38+f/31F82aGXO8o6Ki6NSpExEREXTq1InAwEBOnjxJQIDx8/nrr79o164dX3zxBRMnTqRu3bq2c3366afMmDGD0NBQwsPDadq0qe3YokWLGDp0KI888ggnTpwgKCgoQz3OnDkDwP79+2nUqFG+6//yyy+zZMkSbrnlFhYvXpzhmqmpqSxfvjxD+TZt2rB9+3Y6dOiQIT8yMpLevXuzePFiFixYwAMPPJDvOriaIxMtJbg4SfLy8XglG3/YZyxVOdRkHH1k1Jgo4d58801bcAEwmUy8/fbbAGzatCnXrqmCmDJlii24AAQFBfH0008DcODAAebMmWMLLgCtWrWid+/eaK3ZsGGDLT81NZUpU6YAsGDBggwf9AADBgxgzJgxREdH8+2332Zbl3fffdeh4HL58mVmz56NyWRi4cKFWa7p4eFB//79M+Q1a9YsS3ABqFGjBu+//z4A4eHh+a6DOxR6mLJw0JFf8T74iy35nucYpgxsk8sbhCgZ+vbtmyUvODiYoKAgoqKiuHbtGtWrV3fa9Xr27JklL+1Dvm7duhmCT5rGjRsDcOGCfS74X3/9RWRkJKGhoTRv3jzba3Xt2pWZM2eybds2nnvuuSzHBw4c6FDd161bh9lspkuXLoSGhub7fSkpKaxbt45t27Zx8eJFEhMT0VoTG2usX3jkyBGH6uFqEmDcKSmOxEUvkLay2MLU2+k7+FEqS9dYyebEbibAuKG/9HmjWyyNV3noNwPChjn3Wk5Up06dbPMrVqxIVFQUiYmJTr1eSEhIljx/f/8cj6U/nr4uJ06cAIwurrymB6QNIEgvODiY8uXL56/SVqdPnwbI0nLJzZEjRxgwYAAHDx7MscyNGzccqoerSYBxo6Q1b1EuwfjP6br2Z0fjF5naQrrGRCZpQWTtFIg5BwEh0GNysQ4uYHSJOYvFYinU9RypS9roq1q1anH33XfnWja7gOBocCmoIUOGcPDgQfr378/48eNp1qwZAQEBeHh4cOTIEZo0aZI2IKvYkADjLud347Vrji053WM44wffXnT1EcVb2LBiH1AKw9vbG4C4uLhsj6f9h+8OtWvXBox7GfPmzXPLNdMGGBw+nL9b24cOHWLv3r0EBwezcOFCPDw8Mhw/duyY0+voDM77l0PkLNVM7E9jMWH8V7Yp9VY6DRxHJT/vIq6YEEWjatWqeHt7c+3atWy7nVasWOG2urRv357KlSvz559/uu2D+q677sLLy4utW7fm2uWV5vr16wDUrFkzS3AB+O6775xeR2eQAOMGNzd+TIVoYzmYRO3FukYT6NlCJlSKssvLy4s77rgDgMmTJ2fo2tm8eTOTJ092a10mTZpEamoqAwYMYMeOHVnKxMfH88MPP+QrGORHcHAwTz/9NBaLhcGDB2e5OZ+amsrSpfa9oRo3bozJZGLfvn1s3LgxQ9mvvvqKH374wSn1cjbpInO16yfx2DjVlvzcYxjPDbm3CCskRPEwZcoUNm3axKeffsqGDRsIDQ3l9OnT7N69m9dee802xNkdXnjhBU6fPs1HH31Ehw4dCAsLo2HDhlgsFs6ePcuhQ4dISEhg5cqV2Y5OK4hp06Zx/PhxVqxYQWhoKJ06dSIkJITLly+zd+9eLl++bAu8VatWZezYsXzyySd0796drl27Ur16dfbu3cu+ffuYMGEC7777rlPq5UzSgnElrbn24zi8dRIAByx1uWXABOkaEwLo3Lkza9eupUePHpw9e9bWLfb111/z1ltvub0+H374IRs2bODBBx8kKiqK5cuXs379ehISEujXrx/fffedrdXlDD4+PixdupRvvvmGO++8k3379hEeHs6hQ4cICwtj5syZGcpPnz6dOXPm0LJlS3bs2MHKlSupVq0aK1euZPTo0U6rlzOp4jbqoKi0bdtW79q1y6nnjN/5LX7LxwFg0YqP6s3mxScfcuo1hOsdPHjQaf+1ClESOPI7r5TarbVum90xacG4Svw19KrXbMn/M/VmxLDBRVghIYRwLwkwLnJhwb/wTzUm4J3TVag24C2CpGtMCFGGSIBxgbgDv1Lz9CJbekXtF+nRsuzuCSGEKJskwDhbcgKJv7xgS65RnRn28KgirJAQQhQNCTBOdvLnyVQxG8vBxGhfvPtOI9BXusaEEGVPqQkwSqmHlVKblFIxSqk4pdQupdQ4pZTbvsYbJ3ZT+/CXtvTKGuPo2uZWd11eCCGKlVIRYJRSM4HvgLbAJmANcAvwCRCulMq6toKzWVKJWvAMntblYHarUHo+9rLLLyvcQ4bzi7LCmb/rJT7AKKUGA2OBi0CY1rqv1nog0Bg4CAwEnnV1PQ4u/i91E42F65K0J0n3fUCgnyzDXxqYTKZ8re4rRGlgsVictjJ2iQ8wwATr8yta66NpmVrrS8Az1uSrruwqi75wnLp/f2hLrwt+nM4dO7nqcsLNypUrR0JCQlFXQwi3SEhIcNoWBCU6wCilQoA2QDLwU+bjWusNwHmgOtDRJZWIWIDv5x3xxdjA6BJBdH7cfWsoCdfz9/cnOjpauslEqae1Jjo6Gj8//OGUXAAACplJREFUP6ecr0QHGKC19Xm/1vpmDmV2ZirrPBELsCwai7dOtmVVMcUTcHK50y8lik5QUBApKSlERkaSlJQkgUaUOlprkpKSiIyMJCUlhaCgIKect6Svplzf+pzb7kRnMpV1mtQ1r+NhMWfI87AkGzsRluLNosoak8lE7dq1uX79OmfOnCElJaWoqySE03l6ehIQEEBwcLDT7sGU9ADjb32Oz6VM2pZ5FTIfUEqNBkZDzvuJ58YUeyH7AzHnHD6XKN48PT0JDg4mODi4qKsiRIlR0rvIlPW5QH0WWus5Wuu2Wuu2VatWdfziAbWzPxAQUpDqCCFEqVLSA0ys9dk/lzJpx2JzKVMwPSaDV6bRFl7ljXwhhCjjSnqAOWV9rptLmbRmxqlcyhRM2DDoNwMCagPKeO43Q+6/CCEEJf8ezJ/W51ClVPkcRpK1y1TWucKGSUARQohslOgWjNb6LLAH8AaGZj6ulOoKhGDM8t/m3toJIUTZVqIDjNW71uf3lFKN0jKVUsHALGtyqtZa1voQQgg3KuldZGitw5VSszGWhdmrlPoNMAM9gIrAIoxFL4UQQrhRiQ8wAFrrsUqpzcA4oCvgARwC5gKzpfUihBDuVyoCDIDW+nvg+6KuhxBCCIOSdZUMSqkr5L7kTF6qAFedVB1RcsjPveySn72hrtY625nqEmCcRCm1S2vdtqjrIdxLfu5ll/zs81YaRpEJIYQohiTACCGEcAkJMM4zp6grIIqE/NzLLvnZ50HuwQghhHAJacEIIYRwCQkwQgghXEICTCEopR5WSm1SSsUopeKUUruUUuOUUvJ9LaWUUvOUUjqXx6GirqMoOKVUE6XUC0qpb5VSh5RSFuvPdUg+3iufB5mUmpn87qaUmgmMBRKBtdjXP/sE6KGUGqq1Ti3CKgrX2vL/7d1/7NVVHcfx50vAgNBUShRhgtkGoeaPclQOsJirYKWJ5rKRazrncrhck1h/sZWZ/aKl9ksb/bBaRmqLrfS74lv0S1AspExwfENUUPoxgUBM3/1xzh3X670f7v3e772X7uf12O7O7uece77v7/hy3vfzOedzPsCWOsef7nYgNqKuAa5r9UMeD+pzghkGSReT/ph2AHMiYnM+Pgn4FXARcC3w5Z4FaZ12e0Ss7HUQNuIeAT4HrAceBO4g7W/YkMeDxkp76tamZblcWvljAoiInaRvQACfKPOpsdn/o4i4PSJuiIgfRcTjTX7M40EDpfuF2yVpCnAOcAC4q7Y+IgaBJ4ETgNndjc7MusnjQTFfImvdWbnc1OARzQDrgJNy2991JSrrtvMlnQFMAHYCa4H7/WiI0vF4UMAJpnXTc1m08/K2mrbWfxbXOfYXSZdFxMauR2O94vGggC+RtW5CLvcWtNmTy6M6HIt138PAEmAW6W9hMrAQ+BPwRmBA0km9C8+6zONBAZ/BtE659B47JRQRK2oO7QVWS7ofGCRdZ19GWjVk/c/jQQGfwbRudy4nFLSp1O0uaGN9JCIOAJ/Jb9/Ty1isqzweFHCCad1QLk8uaDO1pq2VQ+Uufl8iK4+hXHo8qMMJpnUbcjlL0rgGbd5S09bKYWIu9xS2sn7i8aCAE0yLIuIJ4CHgSOCS2npJc4EppLt6f9/d6KzHLs3lup5GYV3j8aCYE8zwVK61f1bSqZWDko4Hbstvb/I9Ef1F0pmSFkoaVXN8tKTrSavLAL7U/eishzweNOAHjg2TpNtI20DsBwY4uLnd0cA9wKIybm7XzyRdCNwN/BN4DNhOWnp6Omm58kvAsoi4uWdBWlsknc3BpABp6flRwGbSvzsAETG75nMeD+pwgmmDpA8CHyUNMKNIk7zfAr5axm8r/U7SdNJOu+eSJnUnkpanbgd+A9waEQ/2LkJrl6R5pA0qC0WEao95PHglJxgzM+sIz8GYmVlHOMGYmVlHOMGYmVlHOMGYmVlHOMGYmVlHOMGYmVlHOMGYmVlHOMGYHYKkeZJC0ppex9IuSUvz7/KuNvo4W9JLkj4/krFZ/3GCsdKTNJQH3Wm9jqWTJJ0IfBL4dUT8fLj9RMRDwE+AJZLeMFLxWf9xgjE7tAeAmcDiXgfSpuWkfbWWj1BfYzi40aPZK3irGCs9SUOkvcWmR8RQb6PpDEkTSXumPQWcGiPwH1/SOuAs4JSI2NZuf9Z/fAZjpSXpCknBwacRbs2XyqL6klmjORhJ0/LxIUlHSLpe0iZJ+yRtl/RFSeNz22Mlrchtn5e0OW/x3yg2SbpM0n2SduXPbJP0zWFeyvsIMBb4Tr3kIukYSTfm+P9T9TuskbSsQZ/fJm3qePUw4rESGN3rAMx6aAtpkFwEvBpYxcufRtnKkym/DywE1uR+5wAfA2ZKuhz4A+ny1FrguFz/BUljI+LG6o4kjQF+CLwf2AesB3YCpwFXAhdLuiAi1rcQ34W5HKityEnwt6St6Z/JbfYCJ+Zjs6l/KazS1/tIcztmLxcRfvlV6hfpWekBTGtQPy/Xr6k5Pi0fD9LW7JOr6qYCu3LdRuAuYGxV/YJc9xwwvqbfm3LdIDClpu7aXLcFGN3k7zceOJBfY+vUL859/qy2T9IZyjsa9CvSM1ICmNTrf0e/Dr+XL5GZjYwlEfFU5U2kR+l+L789GbgmIvZX1a8G/kw6q3lz5bik40hPxtwDXBIR26t/SETcAqwGXg+8u8nYZpEm5LdWx1BlUi4HIuK/NT/vxYj4Zb1OIyKAv+a3ZzYZi5WIE4xZ+14A6g3CW3K5PiJ21anfnMvJVcfOB8YBgxHxTIOfN5jLtzYZ3/G5/EeD+gdyuVTShyQd02S/cPApj5MKW1kpeQ7GrH07ar/5Z5U5nO116qrrx1YdOyWXC/IChCKvazK+1+TyuXqVETEo6Wbg48B3gZD0KGm+aFVE/KKg70qfrSQlKwknGLP2HepxuK08LndULv9GWhhQ5I9N9vnvXB7dqEFELJX0NdKE/XnA24GrgKsk3QcsaJBEK33+q8lYrEScYMwOL0/kcmNEXDFCfVYutU0sahQRW4EV+YWk84AfABeQljl/o87HKn02upxnJeY5GLO0ugoOjy9cA6Q5nfktzoUU2QQ8D0yXNK7ZD0XEWmBlfvum2npJAmbktxvajNH6kBOMGTyZy5k9jQKIiJ3AraQ5jZ9KmlHbJt+0eaWkpibWI2If6XLaGOCcOv1dJGmOpCNqjo8D5ue3f6/T9QzgWGBTwYIEK7HD4RubWa/dTbrX5c4831CZs1gaEY1WXnXSDaSVZZcCj0h6GNhKWgwwlZQIj8zlzib7vId0c+d80uR9tbnAdcCzkjYAz5IWBryNdFPoo8DX6/RZST73NhmDlYwTjBncQpqsvpx0N/6r8vFP0Xhpb8dExAvAByTdSZr7OBc4A9gNPE3aNeBe4PEWul0JfBpYLGl5voelum4/aXL/NOC1pCS7hTQHc0dE7K7T54eBF6mffMy82aVZWeRVYlcD72x082QLfZ1OulF0VUQsGon4rP84wZiVhKQTgMeADRExt82+fgy8F5gVEZsP1d7KyZP8ZiURETtIl/3mtPtES9JGnF9xcrEiPoMxM7OO8BmMmZl1hBOMmZl1hBOMmZl1hBOMmZl1hBOMmZl1hBOMmZl1xP8A+Ervg7n5dOIAAAAASUVORK5CYII=\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(t,v_analytical(t,m,g,c),'-',label='analytical')\n", | |
| "plt.plot(t,v_numerical,'o-',label='numerical')\n", | |
| "plt.legend()\n", | |
| "plt.xlabel('time (s)')\n", | |
| "plt.ylabel('velocity (m/s)')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "**Note:** In the above plot, the numerical solution is given at discrete points connected by lines, while the analytical solution is drawn as a line. This is a helpful convention. We plot discrete data such as numerical solutions or measured data as points and lines while analytical solutions are drawn as lines. \n", | |
| "\n", | |
| "## Exercise\n", | |
| "\n", | |
| "Play with the values of `t` (defined above as `t=np.linspace(0,12,7)`). \n", | |
| "\n", | |
| "If you increase the number of time steps from 0 to 12 seconds what happens to v_analytical? to v_numerical?\n", | |
| "\n", | |
| "What happens when you decrease the number of time steps?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "# Errors in Numerical Modeling\n", | |
| "\n", | |
| "## 1 - Truncation\n", | |
| "## 2 - Roundoff " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "# 1- Truncation error\n", | |
| "## Freefall is example of \"truncation error\"\n", | |
| "### Truncation error results from approximating exact mathematical procedure\n", | |
| "\n", | |
| "We approximated the derivative as $\\frac{d v}{d t}\\approx\\frac{\\Delta v}{\\Delta t}$\n", | |
| "\n", | |
| "Can reduce error in two ways\n", | |
| "\n", | |
| "1. Decrease step size -> $\\Delta t$=`delta_time`\n", | |
| "\n", | |
| "2. Increase the accuracy of the approximation" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Truncation error as a Taylor series \n", | |
| "\n", | |
| "The freefall problem solution used a first-order Taylor series approximation\n", | |
| "\n", | |
| "Taylor series:\n", | |
| "$f(x)=f(a)+f'(a)(x-a)+\\frac{f''(a)}{2!}(x-a)^{2}+\\frac{f'''(a)}{3!}(x-a)^{3}+...$\n", | |
| "\n", | |
| "First-order approximation:\n", | |
| "$f(x_{i+1})=f(x_{i})+f'(x_{i})h$\n", | |
| "\n", | |
| "\n", | |
| "We can increase accuracy in a function by adding Taylor series terms:\n", | |
| "\n", | |
| "|Approximation | formula |\n", | |
| "|---|-----------------------------|\n", | |
| "|$0^{th}$-order | $f(x_{i+1})=f(x_{i})+R_{1}$ |\n", | |
| "|$1^{st}$-order | $f(x_{i+1})=f(x_{i})+f'(x_{i})h+R_{2}$ |\n", | |
| "|$2^{nd}$-order | $f(x_{i+1})=f(x_{i})+f'(x_{i})h+\\frac{f''(x_{i})}{2!}h^{2}+R_{3}$|\n", | |
| "|$n^{th}$-order | $f(x_{i+1})=f(x_{i})+f'(x_{i})h+\\frac{f''(x_{i})}{2!}h^{2}+...\\frac{f^{(n)}}{n!}h^{n}+R_{n}$|\n", | |
| "\n", | |
| "Where $R_{n}=O(h^{n+1})$ is the error associated with truncating the approximation at order $n$." | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "source": [ | |
| "In the .gif below, the error in the function is reduced by including higher-order terms in the Taylor series approximation. \n", | |
| "\n", | |
| "\n", | |
| "\n", | |
| "$n^{th}$-order approximation equivalent to \n", | |
| "an $n^{th}$-order polynomial. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "source": [ | |
| "# 2- Roundoff\n", | |
| "\n", | |
| "## Just storing a number in a computer requires rounding\n", | |
| "\n", | |
| "In our analytical solution, $v(t) = v_{terminal}\\tanh{\\left(\\frac{gt}{v_{terminal}}\\right)}$, we can solve for velocity, $v$ at any given time, $t$ by hand to avoid roundoff error, but this is typically more trouble than its worth. Roundoff error comes in two forms:\n", | |
| "\n", | |
| "1. digital representation of a number is rarely exact\n", | |
| "\n", | |
| "2. arithmetic (+,-,/,\\*) causes roundoff error" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "source": [ | |
| "1. digital representation of $\\pi$ \n", | |
| "\n", | |
| "[Consider the number $\\pi$](https://www.piday.org/million/). How many digits can a floating point number in a computer accurately represent?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 369, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "double precision 64 bit pi = 3.141592653589793115997963469\n", | |
| "single precision 32 bit pi = 3.141592741012573242187500000\n", | |
| "First 27 digits of pi = 3.141592653589793238462643383\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "pi=np.pi\n", | |
| "\n", | |
| "double=np.array([pi],dtype='float64')\n", | |
| "single=np.array([pi],dtype='float32')\n", | |
| "print('double precision 64 bit pi = {:1.27f}'.format(double[0])) # 64-bit\n", | |
| "print('single precision 32 bit pi = {:1.27f}'.format(single[0])) # 32-bit\n", | |
| "print('First 27 digits of pi = 3.141592653589793238462643383')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "In order to store the number in a computer you can only use so many bits, shown below is the [64-bit standard for floating point numbers](https://en.wikipedia.org/wiki/Double-precision_floating-point_format):\n", | |
| "\n", | |
| "<img src=\"../images/1236px-IEEE_754_Double_Floating_Point_Format.png\" style=\"width: 400px;\"/> \n", | |
| "\n", | |
| "Where the sign is either + or -, the exponent is a power of two as in, $2^{exponent}$, and the fraction (or base) is the binary representation of the number, $1+\\sum_{i=1}^{52}b_i2^{-i}$. We examine the floating point number representation to highlight that any number you store in a computer is an approximation of the real number you are trying to save. With 64-bit floating point numbers, these approximations are **extremely** good. \n", | |
| "\n", | |
| "2. Floating point arithmetic \n", | |
| "\n", | |
| "Each time you use an operation, e.g. `+ - / *` you lose some precision as well. \n", | |
| "\n", | |
| "Consider $\\pi$ again, but this time we will use a for loop to multiply $\\pi$ by a 1e-16 then divide by 1e-16, then multiply by 2e-16 and divide by 2e-16, and so on until we reach 10e-16. If we do these calculations by hand, we see that each step in the for loop returns $\\pi$, but due to floating point arithmetic errors we accumulate some error. \n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 374, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| " 0 operations 64 bit pi = 3.14159265358979311599796347\n", | |
| "\n", | |
| "20 operations 64 bit pi = 3.14159265358979089555191422\n", | |
| "\n", | |
| "First 26 digits of pi = 3.14159265358979323846264338\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "double=np.array([pi],dtype='float64')\n", | |
| "double_operated=double\n", | |
| "for i in range(0,10):\n", | |
| " double_operated=double_operated*(i+1)*1.0e-16\n", | |
| " double_operated=double_operated*1/(i+1)*1.0e16\n", | |
| "print(' 0 operations 64 bit pi = %1.26f\\n'%double) # 64-bit\n", | |
| "print('20 operations 64 bit pi = %1.26f\\n'%double_operated) # 64-bit after 1000 additions and 1 subtraction\n", | |
| "print('First 26 digits of pi = 3.14159265358979323846264338')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "In the previous block of code, we see $\\pi$ printed for 3 cases:\n", | |
| "\n", | |
| "1. the 64-bit representation of $\\pi$\n", | |
| "\n", | |
| "2. the value of $\\pi$ after it has gone through 20 math operations ($\\times (0..10)10^{-16}$, then $\\times 1/(0..10)10^{16}$)\n", | |
| "\n", | |
| "3. the actual value of $\\pi$ for the first 26 digits\n", | |
| "\n", | |
| "All three (1-3) have the same first 14 digits after the decimal, then we see a divergence between the actual value of $\\pi$ (3), and $\\pi$ as represented by floating point numbers. \n", | |
| "\n", | |
| "We can get an idea for computational limits using some built-in functions:\n", | |
| "\n", | |
| "- `np.info('float64').max`: the largest floating point 64-bit number the computer can represent\n", | |
| "\n", | |
| "- `np.info('float64').tiny`: the smallest non-negative 64-bit number the computer can represent\n", | |
| "\n", | |
| "- `np.info('float64').eps`: the smallest number that can be added to 1" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 20, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "realmax = 1.79769313486231570815e+308\n", | |
| "\n", | |
| "realmin = 2.22507385850720138309e-308\n", | |
| "\n", | |
| "maximum relative error = 2.22044604925031308085e-16\n", | |
| "\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "print('realmax = %1.20e\\n'%np.finfo('float64').max)\n", | |
| "print('realmin = %1.20e\\n'%np.finfo('float64').tiny)\n", | |
| "print('maximum relative error = %1.20e\\n'%np.finfo('float64').eps)\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "## Machine epsilon\n", | |
| "\n", | |
| "The smallest number that can be added to 1 and change the value in a computer is called \"machine epsilon\", $eps$. If your numerical results are supposed to return 0, but instead return $2eps$, have a drink and move on. You won't get any closer to your result. \n", | |
| "\n", | |
| "In the following example, we will add $eps/2$ 1,000$\\times$ to the variable s, set to 1. The result should be $s=1+500\\cdot eps$, but because $eps/2$ is smaller than floating point operations can track, we will get a different result depending upon how we do the addition.\n", | |
| "\n", | |
| "a. We make a `for`-loop and add $eps/2$ 1000 times in the loop\n", | |
| "\n", | |
| "b. We multiply $1000*eps/2$ and add it to the result" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 29, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "fragment" | |
| } | |
| }, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "summation 1+eps/2 over 1000 minus 1 = 0.0\n", | |
| "500.0 *eps= 1.11022302463e-13\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "s1=1;\n", | |
| "N=1000\n", | |
| "eps=np.finfo('float64').eps\n", | |
| "for i in range(1,N):\n", | |
| " s1+=eps/2;\n", | |
| "\n", | |
| "s2=1+500*eps\n", | |
| "print('summation 1+eps/2 over ',N,' minus 1 =',(s-1))\n", | |
| "print(N/2,'*eps=',(s2-1))" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Exercise\n", | |
| "\n", | |
| "1. Try adding $2eps$ to 1 and determine the result of the previous exercise. \n", | |
| "\n", | |
| "2. What is machine epsilon for a 32-bit floating point number?" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "slide" | |
| } | |
| }, | |
| "source": [ | |
| "# Freefall Model (revisited)\n", | |
| "\n", | |
| "In the following example, we judge the **convergence** of our solution with the new knowledge of truncation error and roundoff error. \n", | |
| "\n", | |
| "**The definition for convergence in mathematics is the limit of a sequence exists.** \n", | |
| "\n", | |
| "In the case of the Euler approximation, the sequence is smaller timesteps, $\\Delta t$, should converge to the analytical solution. \n", | |
| "\n", | |
| "Define time from 0 to 12 seconds with `N` timesteps \n", | |
| "function defined as `freefall`\n", | |
| "\n", | |
| "m=60 kg, c=0.25 kg/m" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "source": [ | |
| "## Freefall example\n", | |
| "\n", | |
| "Estimated the function with a $1^{st}$-order approximation, so \n", | |
| "\n", | |
| "$v(t_{i+1})=v(t_{i})+v'(t_{i})(t_{i+1}-t_{i})+R_{1}$\n", | |
| "\n", | |
| "$v'(t_{i})=\\frac{v(t_{i+1})-v(t_{i})}{t_{i+1}-t_{i}}-\\frac{R_{1}}{t_{i+1}-t_{i}}$\n", | |
| "\n", | |
| "$\\frac{R_{1}}{t_{i+1}-t_{i}}=\\frac{v''(\\xi)}{2!}(t_{i+1}-t_{i})$\n", | |
| "\n", | |
| "or the truncation error for a first-order Taylor series approximation is\n", | |
| "\n", | |
| "$R_{1}=O(\\Delta t^{2})$\n" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Computer model error = truncation + roundoff\n", | |
| "\n", | |
| "In the function `freefall(N)`, the speed of a 60-kg object is predicted in two ways:\n", | |
| "\n", | |
| "1. The analytical 64-bit representation, \n", | |
| "$v(t)=v_{terminal}\\tanh{\\left(\\frac{gt}{v_{terminal}}\\right)}$\n", | |
| "\n", | |
| "2. The numerical 32-bit$^{+}$ Euler approximation for `N`-steps from 0 to 2 seconds\n", | |
| "\n", | |
| "$^{+}$Here, we use a 32-bit representation to observe the transition from truncation error to floating point error in a reasonable number of steps. \n", | |
| "\n", | |
| "We can reduce truncation error by decreasing the timestep, $\\Delta t$. Here, we consider the speed from 0 to 2 seconds, so `N=3` means $\\Delta t$= 1 s and `N=21` means $\\Delta t$=0.1 s\n", | |
| "\n", | |
| "|N= | $\\Delta t$=|\n", | |
| "|---|---|\n", | |
| "|3 | 1 s|\n", | |
| "|21| 0.1 s|\n", | |
| "|201| 0.01 s|\n", | |
| "|??| 0.05 s|\n", | |
| "|?? | 0.001 s|\n", | |
| "\n", | |
| "What is N for 0.05 s and 0.001 s in the table above?\n", | |
| "\n", | |
| "Answer (0.05 s): <span style=\"color:white\"> 41 </span>\n", | |
| "\n", | |
| "Answer (0.001 s): <span style=\"color:white\"> 2001 </span>\n", | |
| "\n", | |
| "Highlight lines above for answer." | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 89, | |
| "metadata": { | |
| "slideshow": { | |
| "slide_type": "subslide" | |
| } | |
| }, | |
| "outputs": [], | |
| "source": [ | |
| "def freefall(N):\n", | |
| " ''' help file for freefall(N)\n", | |
| " computes the velocity as a function of time, t, for a\n", | |
| " 60-kg object with zero initial velocity and drag \n", | |
| " coefficient of 0.25 kg/s\n", | |
| " N is number of timesteps between 0 and 2 sec\n", | |
| " v_analytical is the 64-bit floating point \"true\" solution\n", | |
| " v_numerical is the 32-bit approximation of the velocity\n", | |
| " t is the timesteps between 0 and 2 sec, divided into N steps'''\n", | |
| " t=np.linspace(0,2,N)\n", | |
| " c=0.25\n", | |
| " m=60\n", | |
| " g=9.81\n", | |
| " v_terminal=np.sqrt(m*g/c)\n", | |
| "\n", | |
| " v_analytical = v_terminal*np.tanh(g*t/v_terminal);\n", | |
| " v_numerical=np.empty(len(t),dtype=np.float32)\n", | |
| " delta_time =np.diff(t)\n", | |
| " for i in range(0,len(t)-1):\n", | |
| " v_numerical[i+1]=v_numerical[i]+(g-c/m*v_numerical[i]**2)*delta_time[i];\n", | |
| " \n", | |
| " return v_analytical,v_numerical,t\n", | |
| " " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "We can visualize how the approximation approaches the exact solution with this method. The process of approaching the \"true\" solution is called **convergence**. \n", | |
| "\n", | |
| "First, we solve for `n=2` steps, so t=[0,2]. We can time the solution to get a sense of how long the computation will take for larger values of `n`. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 226, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "CPU times: user 288 µs, sys: 3 µs, total: 291 µs\n", | |
| "Wall time: 303 µs\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "%%time\n", | |
| "n=2\n", | |
| "\n", | |
| "v_analytical,v_numerical,t=freefall(n);" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "The block of code above assigned three variables from the function `freefall`. \n", | |
| "\n", | |
| "1. `v_analytical` = $v_{terminal}\\tanh{\\left(\\frac{gt}{v_{terminal}}\\right)}$\n", | |
| "\n", | |
| "2. `v_numerical` = Euler step integration of $\\frac{dv}{dt}= g - \\frac{c}{m}v^2$\n", | |
| "\n", | |
| "3. `t` = timesteps from 0..2 with `n` values, here t=np.array([0,2])\n", | |
| "\n", | |
| "All three variables have the same length, so we can plot them and visually compare `v_analytical` and `v_numerical`. This is the comparison method (1) from above. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 227, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "import matplotlib.pyplot as plt\n", | |
| "%matplotlib inline" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 228, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "<matplotlib.legend.Legend at 0x7f131d486c90>" | |
| ] | |
| }, | |
| "execution_count": 228, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEWCAYAAAB8LwAVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3dd5hTZfbA8e9hZugoXWkDqKgUqUMREEGkCKis6M+CbS1Y0N21FywsNhTsuIvIKrbV3bULCIKKIIIIAqICotQBpPc65fz+eO+QZEiGzJDkTjmf58lDcm47uRNycsv7vqKqGGOMMbmV8jsBY4wxhZMVCGOMMWFZgTDGGBOWFQhjjDFhWYEwxhgTlhUIY4wxYVmBMDEnIqkisltEkvzOpbATkZUicnaM19lJRJZ5f4P+YaafIiLzRWSXiPwlltvO7/pFpIGIqIgke6+nich1sc7JFIwVCFNg3pfbPu+LKOdRW1VXq2pFVc0qwDqvFpFvjjDPSO8LcJeILBGRKwv+LoqlYcAo72/wUZjpdwPTVLWSqr4Qh+3He/0mQaxAmKN1rvdFlPNYl9fM4hzt524PcC5wLHAV8LyIdDzKdRYn9YGfCzo9Bkd+R9q+KSKsQJiYi3Da4DERmQnsBU7wjhSWe0cBK0RkoIg0BkYDp3tHI9vDrV9VH1bVJaqararfATOA0yPkUl1ExovIdhHZKiIzcgqUiNQWkfdFZJOXw1+ClksSkftF5Hcvx3kiUs+b1lFEvheRHd6/HYOWmyYij4jITG+5z0WketD0K0RklYhsEZEhuXJtJyJzRWSniGwQkWfy2MfXi8hv3nv6RERqe/HfgROAT719WCbXcl8C3YBR3vSTRWSciPxTRCaKyB6gm4iU8Y7UVnu5jBaRckHr6SciC7z9+q2INM9j/X29U047RWSNiAyN9L5MIaOq9rBHgR7ASuDsMPEGgALJ3utpwGqgKZCM++W/EzjFm14LaOo9vxr4Jh85lAPWA70jTH8CV3RSvMcZgOB+HM0DHgJK475UlwO9vOXuAhYBp3jztwCqAVWBbcAV3nu51HtdLei9/g6c7OU2DRjuTWsC7Aa6AGWAZ4DMnH0IzAKu8J5XBDpEeE9nAZuB1t56XgSmH+nvEjR9GnBd0OtxwA6gk7dfygLPAZ9477cS8CnwhDd/a2Aj0B5Iwh3FrQTKRFh/V+A0b93NgQ1A/zw+K9dFyt0eiX3YEYQ5Wh95vyK3i0i48905xqnqz6qaiftSzAaaiUg5VV2vqgU9JTEaWAhMjjA9A1eA6qtqhqrOUPdN1BaooarDVPWgqi4HXgEu8Za7DnhAVZeqs1BVtwB9gWWq+qaqZqrqO8AS3CmvHK+p6q+qug/4L9DSi18IjFfV6ap6AHjQ2w/BuZ4kItVVdbeqzo7wngYCr6rqD9567sMddTWIZodF8LGqzlTVbOAAcD1wm6puVdVdwONB++Z64GVV/U5Vs1T1dW+ZDuFWrKrTVHWRuiO+H4F3gDOPIleTIFYgzNHqr6qVvcdhd8wEWZPzRFX3ABcDNwLrRWSCiJya3w2LyAigGfB/3pd+OCOA34DPvVNa93rx+kDtoOK2HbgfOM6bXg93JJBbbWBVrtgqoE7Q6z+Cnu/FHQ3kLJt7P2wJmvda3JHHEu/UVb8I7ykkB1Xd7a2nToT5o7Em6HkNoDwwL2jfTPLi4PbdHbn2XT0vr8OISHsR+co7lbcD93evHm5eU7hYgTCJEvIFrqqTVbUH7tf9Etyv98Pmi0RE/g6cA/RU1Z0RN6q6S1XvUNUTcL/ybxeR7rgvxBVBxa2yurtu+niLrgFODLPKdbgvyGCpwNoo0l6P+yLNeQ/lcaetcnJdpqqXAjWBJ4H3RKTCkXLw5qkWZQ6RBO/3zcA+3Gm/nH1zrKrmFLo1wGO59l1572gqnH/jTlfVU9VjcUd9chS5mgSxAmESTkSOE5HzvC+2A7jz8jm3xG4A6opI6TyWvw+4DOjhnfbJa1v9ROQkERHcdY8s7zEH2Cki94hIOe+idDMRaestOhZ4REQaidNcRKoBE4GTReQyEUkWkYtx1xbGR/HW3wP6iUhn7/0NI+j/oIhcLiI1vNM8ORfow90q/G/gzyLS0rsI/TjwnaqujCKHI/K2/wrwrIjU9HKrIyK9vFleAW70jgxERCp4F6IrRVhlJWCrqu4XkXa4v50pAqxAGD+UAu7A/RLeijsffbM37UvcLZJ/iMjmCMs/jvvVntMYbLeI3B9h3kbAVFwRmgX8wzsnnoU7omgJrMD9ah6Lu4AO7gLyf4HPcYXlX0A5ryD18/Lfgrvnv5+qRsr1EO86y2DcF/x63MXt9KBZegM/i8hu4HngElXdH2Y9X+CuX7zvredEAtcHYuUe3Km52SKyE7cPT/G2Pxd3HWKU9x5+w91cEMnNwDAR2YW7KeC/Mc7VxIlEPnVrjDGmJLMjCGOMMWFZgTDGGBOWFQhjjDFhWYEwxhgTVrLfCcRS9erVtUGDBn6nYYwxRca8efM2q2qNcNOKVYFo0KABc+fO9TsNY4wpMkQkd88Ah9gpJmOMMWHFrUCISD2v/5XFIvKziPzVi1cVkSniBnyZIiJVIix/lTfPMhG5Kl55GmOMCS+eRxCZwB2q2hjXy+NgEWkC3At8oaqNgC+81yFEpCrwMK474XbAw5EKiTHGmPiI2zUIVV2P6wYAVd0lIotxvU2ej+sfHuB1XP/v9+RavBcwRVW3AojIFFw3BJE6A4soIyOD9PR09u8/rMcC44OyZctSt25dUlJS/E7FGHMECblI7fVT3wr4DjjOKx6o6vqczsByqUNo98PpROjKWEQGAYMAUlNTD5uenp5OpUqVaNCgAa6/NuMXVWXLli2kp6fTsGFDv9MxxhxB3C9Si0hFXKdif8urW+bci4WJhe00SlXHqGqaqqbVqHH4nVr79++nWrVqVhwKARGhWrVqdjRnTIx8NH8tnYZ/ScN7J9Bp+Jd8NP9oenw/XFwLhIik4IrD26r6gRfeICK1vOm1cEMX5pZOUL/5QF1cz58FzaOgi5oYs7+FMbHx0fy13PfBItZu34cCa7fv474PFsW0SMTzLibBdZG8WFWDB1//BDeGLd6/H4dZfDLQU0SqeBenexJ5SEljjClxRkxeyr6M0OFC9mVkMWLy0phtI55HEJ1wA7ufJSILvEcfYDjQQ0SWAT2814hImoiMBfAuTj8CfO89huVcsC5q1qxZQ7du3WjcuDFNmzbl+eefDzvf0KFDqVOnDi1btjz02L59e9h5czRo0IDNm484DEG+TJs2jW+//Tam6zTGxN667fvyFS+IeN7F9A2RhxXsHmb+ubiB4nNevwq8Gp/sIvto/lpGTF7Kuu37qF25HHf1OoX+rQo+1G9ycjJPP/00rVu3ZteuXbRp04YePXrQpEmTw+a97bbbuPPOO48m/TxlZWWRlJSU5zzTpk2jYsWKdOzYMW55GGOOXu3K5VgbphjUrlwuZtuwltRB4nFOr1atWrRu3RqASpUq0bhxY9aujX5948aN45Zbbjn0ul+/fkybNu2w+d566y3atWtHy5YtueGGG8jKcoeeFStW5KGHHqJ9+/bMmjUrZJkXXniBJk2a0Lx5cy655BJWrlzJ6NGjefbZZ2nZsiUzZsxg06ZNDBgwgLZt29K2bVtmzpwJuCOeK664grPOOotGjRrxyituSOn169fTpUsXWrZsSbNmzZgxY0a+9pcxJjp39TqFcimhP/jKpSRxV69TYraNYtUX09HK65ze0RxF5Fi5ciXz58+nffv2Yac/++yzvPXWWwBUqVKFr776Kqr1Ll68mP/85z/MnDmTlJQUbr75Zt5++22uvPJK9uzZQ7NmzRg2bNhhyw0fPpwVK1ZQpkwZtm/fTuXKlbnxxhupWLHioSOZyy67jNtuu43OnTuzevVqevXqxeLFiwH48ccfmT17Nnv27KFVq1b07duXd955h169ejFkyBCysrLYu3dvQXaVMeYIcr6TYnnGIzcrEEHieU5v9+7dDBgwgOeee45jjjkm7DwFPcX0xRdfMG/ePNq2bQvAvn37qFnTNS9JSkpiwIABYZdr3rw5AwcOpH///vTv3z/sPFOnTuWXX3459Hrnzp3s2rULgPPPP59y5cpRrlw5unXrxpw5c2jbti3XXHMNGRkZ9O/fn5YtW+b7/RhjotO/VZ2YFoTc7BRTkEjn7o72nF5GRgYDBgxg4MCBXHDBBflaNjk5mezs7EOvw7UhUFWuuuoqFixYwIIFC1i6dClDhw4FXMvlSNcdJkyYwODBg5k3bx5t2rQhMzPzsHmys7OZNWvWoXWvXbuWSpUqAYffsioidOnShenTp1OnTh2uuOIK3njjjXy9X2NM4WEFIkg8zumpKtdeey2NGzfm9ttvz/fyDRo0YMGCBWRnZ7NmzRrmzJlz2Dzdu3fnvffeY+NG16Rk69atrFoVsQdfgEPr69atG0899RTbt29n9+7dVKpU6dARAkDPnj0ZNWrUodcLFiw49Pzjjz9m//79bNmyhWnTptG2bVtWrVpFzZo1uf7667n22mv54Ycf8v2ejTGFg51iChKPc3ozZ87kzTff5LTTTjt0uuXxxx+nT58+h80bfA0C4KOPPqJTp040bNiQ0047jWbNmh264B2sSZMmPProo/Ts2ZPs7GxSUlJ46aWXqF+/fsS8srKyuPzyy9mxYweqym233UblypU599xzufDCC/n444958cUXeeGFFxg8eDDNmzcnMzOTLl26MHr0aADatWtH3759Wb16NQ8++CC1a9fm9ddfZ8SIEaSkpFCxYkU7gjCmCBPVsD1YFElpaWmae8CgxYsX07hxY58yKr6GDh0acjE7P+xvYkwM7dkC816FzndAqfyfFBKReaqaFm6aHUEYY0xRlJ0NC96CKQ/Bvm1QoSa0ie3QOVYgTIHkXAQ3xvhgw88w/nZYMzsQm/IQNO0PZY+N2WasQBhjTFFxcA9MGw6z/wHZQXcdVk6FPiNjWhzACoQxxhQNSybCZ3fDjqChckqlQMdboctdULp8zDdpBcIYYwqz7Wvgs3tg6YTQeP1O0PcZqHlq3DZtBcIYYwqjrAx3KmnacMgI6rKmfDXo+Si0uBTiPL6KNZQrInJ32hdpnnXrAuMqXXfddSHdZERr2rRp9OvXL9/LGWNiZPVseLmLu/AcXBxaXwm3zIWWl8W9OIAdQRQr48aNo1mzZtSuXRuAsWPH+pyRMSZf9m51RWH+m6Hxmk2h3zOQ2iGh6dgRRAL079+fNm3a0LRpU8aMGQO4briHDBlCixYt6NChAxs2bADg008/pX379rRq1Yqzzz77UDzHrl27aNiwIRkZGYDrPK9Bgwb873//Y+7cuQwcOJCWLVuyb98+unbtSk7DwUmTJtG6dWtatGhB9+5uOI45c+bQsWNHWrVqRceOHVm6NHYjURlj8kEV5r8NL7YJLQ4p5aHHI3DD1wkvDhDHIwgReRXoB2xU1WZe7D9ATsdGlYHtqnpYd58ishLYBWQBmZFa+eXb0NjeAha67h0RJ7366qtUrVqVffv20bZtWwYMGMCePXvo0KEDjz32GHfffTevvPIKDzzwAJ07d2b27NmICGPHjuWpp57i6aefPrSuSpUq0bVrVyZMmED//v159913GTBgABdddBEvvfQSI0eOJC0tdHdt2rSJ66+/nunTp9OwYUO2bnWD85166qlMnz6d5ORkpk6dyv3338/7778fn/1jjAlv4xKYcDusmhkaP7Uf9B4Olev5kxfxPcU0DhgFHOqMR1UvznkuIk8Dkb9VoZuqxnY8TZ+88MILfPjhh4AbgnTZsmWULl360Hn+Nm3aMGXKFADS09O5+OKLWb9+PQcPHqRhw4aHre+6667jqaeeon///rz22muHBuuJZPbs2XTp0uXQuqpWrQrAjh07uOqqq1i2bBkicuioxBiTAAf3wvSn4NsXQ9s0HFsPznkKTj28v7ZEi9spJlWdDoQdR1pcP9H/B7wTr+0XFtOmTWPq1KnMmjWLhQsX0qpVK/bv309KSsqh7rKTkpIOdbV96623csstt7Bo0SJefvnlsN17d+rUiZUrV/L111+TlZVFs2bN8sxBVQ/rmhvgwQcfpFu3bvz00098+umnYbdljImDXyfDP9rDN88GikOpZOj0Nxj8XaEoDuDfReozgA2quizCdAU+FxEFXlbVMTHZah6ngeJlx44dVKlShfLly7NkyRJmz559xPnr1HG9x77++usR57vyyiu59NJLefDBBw/FcnfVneP0009n8ODBrFix4tAppqpVq4Zsa9y4cQV4d8aYfNmR7to0LBkfGq/XAfo9C8cdPla9n/y6SH0peR89dFLV1sA5wGAR6RJpRhEZJCJzRWTupk2bYp3nUevduzeZmZk0b96cBx98kA4d8r7QNHToUC666CLOOOMMqlevHnG+gQMHsm3bNi699NJDsauvvpobb7zx0EXqHDVq1GDMmDFccMEFtGjRgosvdmf67r77bu677z46dep0aAxrY0wcZGXCt6NgVLvQ4lCuCpw3Cv78WaErDhDn7r5FpAEwPucitRdLBtYCbVQ1PYp1DAV2q+rII81bkrr7fu+99/j444958803jzxzIVNc/ybGhLVmjutYb8Oi0Hiry+HsYVChmj95eQpbd99nA0siFQcRqQCUUtVd3vOewLBEJljY3XrrrXz22WdMnDjR71SMMZHs3Qpf/B3mjQuN12js2jTU7+hLWvkRz9tc3wG6AtVFJB14WFX/BVxCrtNLIlIbGKuqfYDjgA+9i6rJwL9VdVK88iyKXnzxRb9TMMZEogo//gcmD4G9QTdiJpeDrvdAh8GQXNq//PIhbgVCVS+NEL86TGwd0Md7vhxoEeNcwt7FYxKvOI1gaMxhNv3q2jSsnBEaP7m3u3W1SuRhgAujYt/VRtmyZdmyZQvVqlWzIuEzVWXLli2ULVvW71SMia2MfTB9JMx8HrKD2hMdU8dr09A3IX0nxVqxLxB169YlPT2dwniHU0lUtmxZ6tat63caxsTOsikw8U7YtjIQkyTocBN0vQ/KVPQttaNV7AtESkpK2NbIxhhzVHaug0n3wi8fh8brtnNtGo7PuwFrUVDsC4QxxsRUViZ8/wp8+Sgc3B2Il60MPf4Ora6EUsWjH1QrEMYYE630eTD+r/BHrjYNLS6DHsOgYg1/8ooTKxDGGHMk+7bDF8Ng7qu4noA81U92w342PMO31OLJCoQxxkSiCov+59o07NkYiCeXhS53Qce/FJk2DQVhBcIYY8LZvMy1aVgxPTR+Ug/oMwKqFv+bX6xAGGNMsIz98M0zrivurIOBeKVacM6T0Pi8ItmmoSCsQBhjTI7fvoAJd8C2FYGYlIL2N0K3+6FMJf9y84EVCGOM2bkeJt8PP38QGq/TxrVpqBXT3n+KDCsQxpiSKzsLvv8XfPkIHNgZiJc9Fro/DG2uhlJJvqXnNysQxpiSae0PMP42WL8gNN78Yuj5KFSs6U9ehYgVCGNMybJ/h2sFPecVQto0VDvJtWk44UzfUitsrEAYY0oGVfjpfXetYfeGQDypjGvT0OkvkFzGv/wKISsQxpjib8vv7u6k5V+Fxk88C/qMhGon+pNXIWcFwhhTfGUegG+egxlPQ9aBQLzi8dD7CWj6pxLTpqEg4tbloIi8KiIbReSnoNhQEVkrIgu8R58Iy/YWkaUi8puI3BuvHI0xxdjvX8E/TodpjweKQ06bhlu+h2YXWHE4gngeQYwDRgFv5Io/q6ojIy0kIknAS0APIB34XkQ+UdVf4pWoMaYY2bUBPh/i+lAKVruVa9NQu5U/eRVB8RyTerqINCjAou2A37yxqRGRd4HzASsQxpjIsrNg3mswdRgc2BGIlzkGuj8EadeU6DYNBeHHNYhbRORKYC5wh6puyzW9DrAm6HU60D7SykRkEDAIIDU1NcapGmOKhHULXJuGdT+ExptdCL0eg0rH+5NXEZfoYY/+CZwItATWA0+HmSfcSUENE3MTVMeoapqqptWoUbwG6zDGHMH+nfDZvfBKt9DiUPVEuOIjuPBfVhyOQkKPIFT10M3HIvIKMD7MbOlAvaDXdYF1cU7NGFOUqMIvH8Gk+2DX+kA8qTSccQd0+huklPUvv2IioQVCRGqpas5f80/AT2Fm+x5oJCINgbXAJcBlCUrRGFPYbV0OE++C36aGxk/oCn2ehuon+ZFVsRS3AiEi7wBdgeoikg48DHQVkZa4U0YrgRu8eWsDY1W1j6pmisgtwGQgCXhVVX+OV57GmCIi8wDMfAFmjITM/YF4hZquTUOzAXbbaoyJasTT+0VOWlqazp071+80jDGxtmI6jL8dtiwLCgq0vQ7OegDKVfYttaJOROapalq4adaS2hhTeO3eBJ8/AD++Gxqv1cK1aajTxp+8SggrEMaYwic7G34YB1OHut5Xc5SuBN0fdEcO1qYh7qxAGGMKlz8WuTYN6d+Hxpv+CXo9AcfU8ievEsgKhDGmcDiwC756Ar4bDZoViFdpAH2fhpPO9i21ksoKhDHGX6qw+FP47B7YFdTkqVQKdL4NzrgdUsr5l18JZgXCGOOfbSth4t2wbHJovMEZbnS3Gif7kpZxrEAYYxIv8yDMehG+HgGZ+wLxCjWg1+Nw2kXWpqEQsAJhjEmslTPdRejNS4OCAml/dr2ulqviW2omlBUIY0xi7NkMUx6CBW+Hxo8/Dfo9B3XDttUyPrICYYyJr+xsmP+mKw77twfipStCtyHQbhAk2VdRYWR/FWNM/Gz42Z1OWvNdaLzxedB7OBxbx5+8TFSsQBhjYu/Abvh6OMz6R2ibhsr1oc9IOLmnf7mZqFmBMMbE1pIJ7tbVnemBWKkU6PQXOONOKF3ev9xMvliBMMbExvbVrjD8+llovH5n6PcM1DjFn7xMgVmBMMYcnawMmPUSfP0kZOwNxMtXg56PQYtLrE1DEWUFwhhTcKtmwYTbYeMvofHWV8HZQ6F8VT+yMjESzxHlXgX6ARtVtZkXGwGcCxwEfgf+rKrbwyy7EtgFZAGZkQazMMb4ZM8WmPoQzH8rNF6zqRunIbW9P3mZmCoVx3WPA3rnik0Bmqlqc+BX4L48lu+mqi2tOBhTiGRnu6IwKi20OKRUgJ6Pwg1fW3EoRuJ2BKGq00WkQa7Y50EvZwMXxmv7xpgY27jYDfu5+tvQ+Kn9XJuGyvX8ycvEzRELhIjUBDoBtYF9wE/AXFXNPsptXwP8J8I0BT4XEQVeVtUxeeQ3CBgEkJqaepQpGWMOc3APfP0UzBoF2ZmB+LGp0OcpOOUc/3IzcRWxQIhIN+BeoCowH9gIlAX6AyeKyHvA06q6M78bFZEhQCbwdoRZOqnqOq84TRGRJao6PdyMXvEYA5CWlqb5zcUYk4elk2DiXbBjdSBWKhlOvwXOvBtKV/AvNxN3eR1B9AGuV9XVuSeISDLuAnQP4P38bFBErvKW7a6qYb/QVXWd9+9GEfkQaAeELRDGmDjYke4G8FkyPjSeerq7CF2zsT95mYSKWCBU9a48pmUCH+V3YyLSG7gHOFNV90aYpwJQSlV3ec97AsPyuy1jTAFkZbghP796AjL2BOLlqkLPR6DFZVAqnve2mMLkiH9pEfmriBwjzr9E5AcROWJHKiLyDjALOEVE0kXkWmAUUAl32miBiIz25q0tIhO9RY8DvhGRhcAcYIKqTirg+zPGRGvNHBjTFT5/ILQ4tLoCbp0HrS634lDCRHMX0zWq+ryI9AJqAH8GXgM+z2shVb00TPhfEeZdhzulhaouB1pEkZcxJhb2boWpQ+GH10PjNRq700n1T/clLeO/aApEThv5PsBrqrpQxNrNG1PkqcLCd+HzIbB3SyCeUh7OvAdOHwxJKf7lZ3wXTYGYJyKfAw2B+0SkEnC0t7gaY/y0aalr07Dqm9D4yee4W1cr2y3jJu/bXJO9i9HXAi2B5aq6V0Sq4U4zGWOKmoN7YcZImPkCZGcE4sfUdYXh1L7+5WYKnbyOIGaLSDowCZiU02eSqm4BtuSxnDGmMPr1c5h4J2xfFYhJkjuVdOY9UKaif7mZQimv21zTRKQ+cA7wnIjUAb4BPgO+VtUDCcrRGHM0dqyFSffC4k9C4/XaQ99n4Phm/uRlCr08r0Go6ipgNDBaRFKAM3Ad8D0qIptU1Y5HjSmssjJhzhj46jE4uDsQL1cFegyDlnbbqslb1J31qWoG8KX3wDuiMMYURulzYfzf4I9FofGWA11xqFDdn7xMkRJNZ339gEeABkAS7rZXVdVj4puaMSbf9m2DL4bB3NdwfV56qp/ihv1s0Nm31EzRE80RxHPABcCiSH0nGWN8pgo//te1adizKRBPLuc61Tv9Fkgu7V9+pkiKpkCsAX6y4mBMIbV5GYy/DVbOCI036gl9RkCVBr6kZYq+aArE3cBEEfkaOHTnkqo+E7esjDFHlrEPZjwDM5+DrIOBeKXacM6T0PhcsE4PzFGIpkA8BuzGjQVhx6jGFAa/TYUJd8K2FYGYJEGHm6DrvVCmkn+5mWIjmgJRVVWP2HurMSYBdq6HyffBzx+GxuukuY71ajX3Jy9TLEVTIKaKSM9c40kbYxIpOwvmvAJfPgoHdwXiZY+Fs4dC66utTYOJuWgKxGDgbhE5AGRgt7kak1hr57mL0OsXhsabX+IG8alY05+8TLF3xAKhqnYy0xg/7Nvujhi+H0tIm4ZqjVybhoZdfEvNlAwRj0lFpEFeC3ojzNU9wjyvishGEfkpKFZVRKaIyDLv3yoRlr3Km2eZN461MSWDKix6D0a1he9f4VBxSC4LZz0AN8204mASIq+TliNE5H0RuVJEmopITRFJFZGzROQRYCZwpJHLx+H6bgp2L/CFqjYCvvBehxCRqsDDQHugHfBwpEJiTLGy5Xd4sz+8fy3s2RiIn3Q23DwLutwFyWX8y8+UKHn15nqRiDQBBgLXALWAvcBiYCLwmKruz2vlqjo9zJHI+UBX7/nrwDTgnlzz9AKmqOpWABGZgis07xzpDRlTJGXsd+0ZZjwDWUEdJVeqBb2fgCb9rU2DSbgj9eb6CzAkxts8TlXXe+tfLyLhrrDVwbXgzpHuxQ4jIueJveEAABmKSURBVIOAQQCpqTYKlimCfv8SJtwBW5cHYlIK2t0A3e6HsnY/iPFH1L25Jli4n0phu/pQ1THAGIC0tDTrDsQUHbs2wOT74af3QuO1W7s2DbVb+pOXMR4/CsQGEanlHT3UAjaGmSedwGkogLq4U1HGFH3ZWTD3Vdfr6oGdgXiZY6H7g5B2DZRK8i8/Yzx+FIhPgKuA4d6/H4eZZzLweNCF6Z7AfYlJz5g4WjfftWlYNz80ftpF0PMxqHScP3kZE0bEAiEirfNaUFV/ONLKReQd3JFAdW9864dxheG/InItsBq4yJs3DbhRVa9T1a3enVLfe6salnPB2pgiaf8O+PIxd9uqZgfiVU+Evk/Did38y82YCCRSL94i8pX3tCyQBizEXRtoDnynqoVu5JG0tDSdO3eu32kYE6Dq+k2adB/s/iMQTyoDZ9wBnf4KKWX9y8+UeCIyT1XTwk3L6zbXbt7C7wKDVHWR97oZcGc8EjWmWNm63PW4+vsXofETurmjhmon+pOXMVGK5hrEqTnFAUBVfxIRu73CmEgyD8DM52H6yNA2DRWPc20aml5gbRpMkRBNgVgsImOBt3C3ml6OayxnjMlt+deuTcOWZUFBgXbXu24yyh7rW2rG5Fc0BeLPwE3AX73X04F/xi0jY4qi3Rvh8wfgx/+Exmu1dG0a6uR5z4cxhVI0vbnuF5HRwERVXZqAnIwpOrKzYd5r8MXf3Z1KOUpXgu4PQdtrrU2DKbKOWCBE5DxgBG640Ybe9YdhqnpevJMzplBb/6Nr07A2151zTS+AXo/DMbX8ycuYGInmFNPDuB5VpwGo6oIjdQVuTLF2YBd89QR898/QNg1VGrq7k07q7l9uxsRQNAUiU1V3iN11YUo6VVj8CXx2L+xaF4gnlYbOt7lHSjn/8jMmxqIpED+JyGVAkog0Av4CfBvftIwpZLathIl3wbJcQ7M3PBP6PgPVT/IlLWPiKZoCcSuuy+8DwL9x/SQ9Gs+kjCk0Mg/Cty/A9BGQGTT8SYUa0OsJOO1Ca9Ngiq1o7mLaCwwRkcdVdU8CcjKmcFj5DYy/HTYH37wnrrfV7g9CORvk0BRv0dzF1BEYC1QEUkWkBXCDqt4c7+SM8cXuTTDlQViYawDD40+Dfs9D3Tb+5GVMgkVziulZ3BCgnwCo6kIRsRHTTfGTnQ3z34ApD8P+7YF46YquFXTb6yGpsI6xZUzsRfVpV9U1ue5iyopPOsb45I+fXJuG9Dmh8Sb9Xf9Jx9T2Jy9jfBRNgVjjnWZSESmNu4vJ+mIyxcOB3TDtCZj9T9Cg3z2V67s2DY16+JebMT6LpkDcCDwP1AHW4u5iGhzPpIyJO1VYMgE+uxt2rg3ES6W4MRq63GltGkyJF81dTJuBgbHaoIicAgT3aHYC8JCqPhc0T1fcUKQrvNAHqjosVjmYEm7bKlcYfp0UGm9whjtqqHGKP3kZU8hEcxfTCbgjiA647r5nAbep6vKCbNDr8K+lt+4k3FHJh2FmnaGq/QqyDWPCysqAWaNg2pOQuS8QL18dej0GzS+2Ng3GBInmFNO/gZeAP3mvLwHeAdrHYPvdgd9VdVUM1mVMZKu+dW0aNuW6fNbmauj+MJSv6ktaxhRm0RQIUdU3g16/JSK3xGj7OcUmnNNFZCGwDrhTVX8Om5zIIGAQQGpqaozSMsXGni0w5SFY8FZo/LhmbpyGeu38ycuYIkBUNe8ZRIYD24F3caeYLgbK4I4qUNWtBdqwuyNqHdBUVTfkmnYMkK2qu0WkD/C8qjY60jrT0tJ07ty5R5rNlATZ2bDgbdfgbd+2QDylAnS7H9rfaG0ajAFEZJ6qpoWbFs3/kIu9f2/IFb8GVzBOKGBe5wA/5C4OAKq6M+j5RBH5h4hU9y6YG5O3Db/AhNth9azQeONzofdwOLauP3kZU8REcxdTwzht+1IinF4SkeOBDaqqItIOKAVsiVMeprg4uAe+fhJmvQTZmYH4sanQZwSc0tu/3IwpgqK5i+kiYJKq7hKRB4DWwCOqOr+gGxWR8kAPgo5KRORGAFUdDVwI3CQimcA+4BI90rkwU7ItmehuXd2xJhArlQwdb4Uud0Pp8v7lZkwRFc0ppgdV9X8i0hnXJ9NIYDRHcReT10NstVyx0UHPRwGjCrp+U4JsXwOf3QNLJ4TGUztCv2egZmN/8jKmGIimQOT0P9AX+KeqfiwiQ+OXkjFRyMpw3WNMewIy9gbi5apCz0eh5WXWpsGYoxRNgVgrIi8DZwNPikgZ3DUBY/yxerbrWG/jL6Hx1lfC2X+3Ng3GxEg0BeL/gN7ASFXdLiK1gLvim5YxYezdClMfhh/eCI3XbOLaNKR28CcvY4qpaEeU+yDo9XpgfTyTMiaEqhu85/MHYG/QzWwp5aHrvdDhZkhK8S8/Y4opaylkCreNS1ybhlUzQ+On9IVznoTK9fzJy5gSwAqEKZwO7oXpI+DbF0LbNBxTF/o8Baf29S83Y0oIKxCm8Pl1Mky8E7avDsRKJbtTSWfeA2Uq+pebMSWIFQhTeOxYC5PugcWfhsbrdXBtGo5r6k9expRQViCM/7IyYc7L8NXjcHB3IF6uCvQYBi0vh1J2Z7UxiWYFwvhrzfeuTcOGRaHxlpe74lChWvjljDFxZwXC+GPfNpj6d5g3DtcpsKfGqdD3GWjQya/MjDEeKxAmsVThx//C5Pthb1Dv7cnloOs90GEwJJf2Lz9jzCFWIEzibPrVtWlYOSM03qiX6467Sn1/8jLGhGUFwsRfxj6Y8TR88xxkZwTix9Rxjd1O7Wcd6xlTCFmBMPG1bCpMvAO2rQzEJAk63ARd77M2DcYUYlYgTHzsXAeT7oNfPgqN123rOtY7/jR/8jLGRM23AiEiK4FduPEmMnMPmi0iAjwP9AH2Aler6g+JztPkU1YmfD8WvnwUDu4KxMtWhh5/h1ZXWpsGY4oIv48guqnq5gjTzgEaeY/2wD85ilHsTAKkz4Pxf4M/fgyNt7gUejwCFWv4k5cxpkD8LhB5OR94wxuLeraIVBaRWl5346Yw2bcdvnwEvv8XIW0aqp/s2jQ0PMO31IwxBedngVDgcxFR4GVVHZNreh0gaAR60r2YFYjCQhUWvefaNOzZGIgnl4Uud0HHv1ibBmOKMD8LRCdVXSciNYEpIrJEVacHTQ9336PmDojIIGAQQGpqanwyNYfb/Jtr07Di69D4ST1cm4aqDf3JyxgTM74VCFVd5/27UUQ+BNoBwQUiHQgeDaYusC7MesYAYwDS0tIOKyAmxjL2wzfPwDfPQtbBQLxSLdemofF51qbBmGLClwIhIhWAUqq6y3veExiWa7ZPgFtE5F3cxekddv3BZ7994cZp2Lo8EJNS0P5G16ah7DH+5WaMiTm/jiCOAz50d7KSDPxbVSeJyI0AqjoamIi7xfU33G2uf/YpV7PrD3ed4af3Q+N12rg2DbVa+JOXMSaufCkQqrocOOxbxSsMOc8VGJzIvEwu2VnuzqQvH4EDOwPxMsfC2Q9Dm6uhVJJv6Rlj4qsw3+Zq/LRuPnz6N1i/IDR+2v9Br8egYk1/8jLGJIwVCBNq/w7XCvr7saDZgXi1k6Dv03BCV78yM8YkmBUI46jCzx+4/pN2bwjEk8pAlzuh018huYx/+RljEs4KhIEtv7u7k37/MjR+4lnQZyRUO9GfvIwxvrICUZJlHnBjNMx4GrIOBOIVj4fej0PTC6xNgzElmBWIkmr5NJhwB2z5LRCTUtD2ejhrCJQ91rfUjDGFgxWIkmb3Rpg8BBb9NzReu5Vr01C7lT95GWMKHSsQJUV2Fsx7DaYOgwM7AvEyx0D3hyDtGmvTYIwJYQWiJFi/EMbfBmvnhcabDYBej0Ol4/3JyxhTqFmBKM7274SvHoc5L4e2aah6gmvTcOJZ/uVmjCn0rEAUR6rwy8cw6V7YFdS/YVJp6Hw7dL4NUsr6l58xpkiwAlHcbF3h2jT8NjU0fkJX6PM0VD/Jj6yMMUWQFYjiIvMAfPsCTB8JmfsD8Qo1ofcT7nqDtWkwxuSDFYjiYMUMN7rb5l+DggJtr4WzHoRylX1LzRhTdFmBKMp2b4LPH4Af3w2NH98czn3OjddgjDEFZAWiKMrOhh9eh6lDYf/2QLx0JTjrAWh7HSTZn9YYc3TsW6So+WORa9OQ/n1ovEl/d63hmNr+5GWMKXYSXiBEpB7wBnA8kA2MUdXnc83TFfgYWOGFPlDV3GNWlywHdsO0J2D2P0GzAvEqDdzdSY3O9i01Y0zx5McRRCZwh6r+ICKVgHkiMkVVf8k13wxV7edDfoWLKiwZD5/dAzvXBuKlUqDz3+CMOyClnH/5GWOKrYQXCFVdD6z3nu8SkcVAHSB3gTDbVsFnd8Ovk0LjDc6Avs9AjZP9ycsYUyL4eg1CRBoArYDvwkw+XUQWAuuAO1X15wjrGAQMAkhNTY1PoomWeRBmjYKvn4LMfYF4+equ76Tm/2dtGowxcedbgRCRisD7wN9UdWeuyT8A9VV1t4j0AT4CGoVbj6qOAcYApKWlaRxTToyVM12bhk1LgoICba6Gsx+GclX8yswYU8L4UiBEJAVXHN5W1Q9yTw8uGKo6UUT+ISLVVXVzIvNMqD2bYcpDsODt0Phxp7lxGuq19ScvY0yJ5cddTAL8C1isqs9EmOd4YIOqqoi0A0oBWxKYZuJkZ8OCt1xx2LctEC9dEbrdD+1usDYNxhhf+PHN0wm4AlgkIgu82P1AKoCqjgYuBG4SkUxgH3CJqhb900e5bfgZxt8Oa2aHxhufB72Hw7F1/MnLGGPw5y6mb4A8r7Cq6ihgVGIy8sHBPTBtOMz+B2RnBuKVU6HPSDi5l3+5GWOMx85dJNqSie7W1R1rArFSKdDpL3DGnVC6vH+5GWNMECsQibJ9tWvstnRiaLx+J9emoeap/uRljDERWIGIt6wMdypp2nDI2BuIl68GPR+FFpdamwZjTKFkBSKeVs92HettzNVIvPVVcPZQKF/Vj6yMMSYqViDiYe9Wd9vq/DdD4zWbQr9nILWDP3kZY0w+WIGIJVVY8G83iM++rYF4Snnoeh90uAmSUvzLzxhj8sEKRKxsXOzaNKz+NjR+aj/XpqFyPX/yMsaYArICcbQO7oXpT8G3L4a2aTi2HpzzFJzax7/cjDHmKFiBOBpLJ8Fnd7lbWHOUSobTb4Ez74bSFfzLzRhjjpIViILYke7aNCwZHxpPPd21aTiuiT95GWNMDFmByI+sDPhuNHz1BGTsCcTLVYEej0DLgVCqlH/5GWNMDFmBiNaaOa5Nw4afQuOtLoezh0GFav7kZYwxcWIF4kj2boUv/g7zxoXGazR2bRrqd/QlLWOMiTcrEJGowsJ3XZuGvUHjFCWXg673QIfBkFzav/yMMSbOrECEs2kpTLgDVs4IjZ/c2926WqW+P3kZY0wCWYEIlrEPpo+Emc9DdkYgfkwdr01DX+tYzxhTYvg1JnVv4HkgCRirqsNzTS8DvAG0wQ01erGqroxHLh/NX8uIyUtptHMWj5V5nTq6ISiRJDj9ZjjzXihTMR6bN8aYQsuPMamTgJeAHkA68L2IfKKqwV2eXgtsU9WTROQS4Eng4ljn8tH8tTz7wdfcz2v0LT0Hggc1rdsO+j0LxzeL9WaNMaZI8OOm/XbAb6q6XFUPAu8C5+ea53zgde/5e0B3kdif23ll0hw+LXU7fZPmHIpt1woMT74JrplsxcEYU6L5USDqAEHjbZLuxcLOo6qZwA4gbEMDERkkInNFZO6mTZvylcgvO0ozMav9odfvZXXhrANP8/LuM6zBmzGmxPPjGkS4IwEtwDwuqDoGGAOQlpYWdp5Ialcux/Dtl9Kw1B88m3khs7NdFxl1KpfLz2qMMaZY8uNncjoQ3Pd1XWBdpHlEJBk4FthKjN3V6xQOpFTm4oMPHSoO5VKSuKvXKbHelDHGFDl+FIjvgUYi0lBESgOXAJ/kmucT4Crv+YXAl6qar6ODaPRvVYcnLjiNOpXLIbgjhycuOI3+rXKf8TLGmJIn4aeYVDVTRG4BJuNuc31VVX8WkWHAXFX9BPgX8KaI/IY7crgkXvn0b1XHCoIxxoThSzsIVZ0ITMwVeyjo+X7gokTnZYwxJsBu1THGGBOWFQhjjDFhWYEwxhgTlhUIY4wxYUkc7h71jYhsAlYVcPHqwOYjzpV4llf+WF75Y3nlT3HMq76q1gg3oVgViKMhInNVNc3vPHKzvPLH8sofyyt/SlpedorJGGNMWFYgjDHGhGUFImCM3wlEYHnlj+WVP5ZX/pSovOwahDHGmLDsCMIYY0xYViCMMcaEVewLhIj0FpGlIvKbiNwbZnoZEfmPN/07EWkQNO0+L75URHolOK/bReQXEflRRL4QkfpB07JEZIH3yN1VerzzulpENgVt/7qgaVeJyDLvcVXuZeOc17NBOf0qItuDpsVzf70qIhtF5KcI00VEXvDy/lFEWgdNi+f+OlJeA718fhSRb0WkRdC0lSKyyNtfcxOcV1cR2RH093ooaFqen4E453VXUE4/eZ+pqt60eO6veiLylYgsFpGfReSvYeaJ32dMVYvtA9ed+O/ACUBpYCHQJNc8NwOjveeXAP/xnjfx5i8DNPTWk5TAvLoB5b3nN+Xk5b3e7eP+uhoYFWbZqsBy798q3vMqicor1/y34rqRj+v+8tbdBWgN/BRheh/gM9woiR2A7+K9v6LMq2PO9oBzcvLyXq8Eqvu0v7oC44/2MxDrvHLNey5ujJpE7K9aQGvveSXg1zD/J+P2GSvuRxDtgN9UdbmqHgTeBc7PNc/5wOve8/eA7iIiXvxdVT2gqiuA37z1JSQvVf1KVfd6L2fjRt6Lt2j2VyS9gCmqulVVtwFTgN4+5XUp8E6Mtp0nVZ1O3qMdng+8oc5soLKI1CK+++uIeanqt952IXGfr2j2VyRH89mMdV6J/HytV9UfvOe7gMVA7gFs4vYZK+4Fog6wJuh1Oofv3EPzqGomsAOoFuWy8cwr2LW4Xwg5yorIXBGZLSL9Y5RTfvIa4B3KviciOcPHFor95Z2Kawh8GRSO1/6KRqTc47m/8iv350uBz0VknogM8iGf00VkoYh8JiJNvVih2F8iUh73Jft+UDgh+0vc6e9WwHe5JsXtM+bLgEEJJGFiue/rjTRPNMsWVNTrFpHLgTTgzKBwqqquE5ETgC9FZJGq/p6gvD4F3lHVAyJyI+7o66wol41nXjkuAd5T1aygWLz2VzT8+HxFTUS64QpE56BwJ29/1QSmiMgS7xd2IvyA6xtot4j0AT4CGlFI9hfu9NJMVQ0+2oj7/hKRirii9DdV3Zl7cphFYvIZK+5HEOlAvaDXdYF1keYRkWTgWNyhZjTLxjMvRORsYAhwnqoeyImr6jrv3+XANNyvioTkpapbgnJ5BWgT7bLxzCvIJeQ6/I/j/opGpNzjub+iIiLNgbHA+aq6JScetL82Ah8Su1OrR6SqO1V1t/d8IpAiItUpBPvLk9fnKy77S0RScMXhbVX9IMws8fuMxePCSmF54I6QluNOOeRc2Gqaa57BhF6k/q/3vCmhF6mXE7uL1NHk1Qp3Ua5RrngVoIz3vDqwjBhdrIsyr1pBz/8EzNbABbEVXn5VvOdVE5WXN98puAuGkoj9FbSNBkS+6NqX0AuIc+K9v6LMKxV3Xa1jrngFoFLQ82+B3gnM6/icvx/ui3a1t++i+gzEKy9ves6PxwqJ2l/ee38DeC6PeeL2GYvZzi2sD9wV/l9xX7ZDvNgw3K9ygLLA/7z/LHOAE4KWHeIttxQ4J8F5TQU2AAu8xydevCOwyPsPsgi4NsF5PQH87G3/K+DUoGWv8fbjb8CfE5mX93ooMDzXcvHeX+8A64EM3C+2a4EbgRu96QK85OW9CEhL0P46Ul5jgW1Bn6+5XvwEb18t9P7OQxKc1y1Bn6/ZBBWwcJ+BROXlzXM17saV4OXivb86404L/Rj0t+qTqM+YdbVhjDEmrOJ+DcIYY0wBWYEwxhgTlhUIY4wxYVmBMMYYE5YVCGOMMWFZgTAmAhGpLCI3B72uLSLvxWlb/YN7Lg0z/TQRGRePbRsTid3makwEXt8341W1WQK29S2uTcfmPOaZClyjqqvjnY8xYEcQxuRlOHCi18//CBFpkDNegLhxMT4SkU9FZIWI3CJuDI/5XqeAOWMFnCgik7yO3GaIyKm5NyIiJwMHcoqDiFzkjTmwUESC+/T5FNfa35iEsAJhTGT3Ar+raktVvSvM9GbAZbguIR4D9qpqK2AWcKU3zxjgVlVtA9wJ/CPMejrhOqnL8RDQS1VbAOcFxecCZxzF+zEmX4p7b67GxNNX6vro3yUiO3C/8MF1d9Dc64GzI/A/N8QI4Pr2yq0WsCno9UxgnIj8FwjunG0jUDuG+RuTJysQxhTcgaDn2UGvs3H/t0oB21W15RHWsw/XERwAqnqjiLTHdcK2QERaqutttaw3rzEJYaeYjIlsF26YxwJR12//ChG5CA6NHdwizKyLgZNyXojIiar6nao+BGwm0GXzyUDYMZONiQcrEMZE4P1qn+ldMB5RwNUMBK4VkZzePsMNkzkdaCWB81AjRGSRd0F8Oq6nUHDjlE8oYB7G5Jvd5mpMISAizwOfqurUCNPLAF8DndUNjWtM3NkRhDGFw+NA+TympwL3WnEwiWRHEMYYY8KyIwhjjDFhWYEwxhgTlhUIY4wxYVmBMMYYE5YVCGOMMWH9Pz+7VJSzcPGwAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.plot(t,v_numerical,'o',label=str(n)+' Euler steps')\n", | |
| "plt.plot(t,v_analytical,label='analytical')\n", | |
| "plt.title('First 2 seconds of freefall')\n", | |
| "plt.xlabel('time (s)')\n", | |
| "plt.ylabel('speed (m/s)')\n", | |
| "plt.legend()" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Exercise\n", | |
| "\n", | |
| "Try adjusting `n` in the code above to watch the solution converge. You should notice the Euler approximation becomes almost indistinguishable from the analytical solution as `n` increases. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## Convergence of a numerical model\n", | |
| "\n", | |
| "You should see that the more time steps we use, the closer the Euler approximation resembles the analytical solution. This is true only to a point, due to **roundoff error**. In our `freefall` function, the numerical result is saved as a 32-bit floating point array. Our best approximation of the freefall function is the 64-bit analytical equation `v_terminal*np.tanh(g*t/v_terminal)`.$^{+}$ \n", | |
| "\n", | |
| "In the next plot, we consider the relative error for the velocity at t=2 s, as a function of `N`. \n", | |
| "\n", | |
| "$^+$ Note: In practice, there is no reason to restrict the precision of our floating point numbers. The function was written this way to highlight the effect of roundoff error without significant computational resources. You would need significantly more timesteps to observe floating point error with 64-bit floating point numbers. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 232, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "N=100\n", | |
| "error=np.zeros(N) # initialize an N-valued array of relative errors\n", | |
| "n=np.logspace(2,5,N) # create an array from 10^2 to 10^5 with N values\n", | |
| "for i in range(0,N):\n", | |
| " v_an,v_num,t=freefall(n[i]) # return the analytical and numerical solutions to our equation\n", | |
| " error[i]=(v_num[-1]-v_an[-1])/v_an[-1] #calculate relative error in velocity at final time t=2 s\n", | |
| "\n", | |
| " " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 233, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "data": { | |
| "text/plain": [ | |
| "Text(0.5, 1.0, 'Truncation and roundoff error \\naccumulation in log-log plot')" | |
| ] | |
| }, | |
| "execution_count": 233, | |
| "metadata": {}, | |
| "output_type": "execute_result" | |
| }, | |
| { | |
| "data": { | |
| "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEpCAYAAACN9mVQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5xcdX3/8debEGQBycpFbZZggNAoiiWan/ZHtAVqBSuXSNFApS0q0Npii22poVpRa00qtkXFSkElRRS5GFMCSKoEykVUEhO5Gq5BshGQy4bbKiF8+sf5DpkMM2fO7M7s3N7Px2Mf2TlzzpnP7NnMZ8/38vkqIjAzM6tlq3YHYGZmnc2JwszMcjlRmJlZLicKMzPL5URhZma5nCjMzCyXE4V1NUn/KOmsdseRR9LWkkLS9HbHUqnR2CQtlPSopHXp8VGS1kl6StK+rYzV2keeR2GSnip7uB3wa2BTevxnEfGNiY/qxSS9DfhKRExvdyyNkLQ1sBHYIyLWtjmcLTQSm6Q9gNuA3SPikbTtfuAvIuLyVsdq7bN1uwOw9ouIHUrfS1oLHB8R36+1v6StI+K5iYit30jaCiAinm93LFW8Cni4LElsBUwjSx7jJmlSRGyqt63OOfy72QJuerK6JH1a0oWSLpD0JHCspPMlfaJsn7elJFN6vE7S30i6RdKGdOxLyp4/UtJqSU9IulvS29P24yXdIelJSfdIOj5tnwIsBXZPzRxPSXp5im1R2XnnSrpN0oik5ZJmFo2p4j3vLenq1MzyiKSvpxiKvr/5kh6UNAz8aZ2f7/WS/knSjcDT6T3uJukySY9JukvS+8v2H+/PvmZskgbT+X8paa2kU5U5BPhu2c//XOAJQMBtktbUeG/7SPp+eh8/k/SHFe/jS5KulPQ08NYa26rGlM5xvKRrJX1B0mPAx/J+1jZGEeEvf73wBawF3lax7dPAs8BhZH9cDADnA58o2+dtwNqyx+uAHwKvBHYG7iS7UwHYHxgBfi+dbxowMz13GLAn2QfQQcAo8Ppqr1EW26L0/WuAp9Jxk4F/SK87uV5MVX4Ov5ni2wZ4OXAD8LmC7+9Q4BfAPsD2wEVAANNrvNb16ef+mhT31un1vghsC7wBeAT43bT/eH72ubEB3wQWAy9N1+Fu4E9rvM7Wdd7XS4Fh4E/Svm8EHi271ucDjwP/P/0evKTGtryYjgeeAz4ITAIG2v1/qBe/fEdhRV0fEUsj4vmIGC14zBkR8WBEPApcBuyXtn8AOCcirkrneyAi1gCk17g3MsuBq4C3Fny9o4FLI2J5RGwEFgI7Am8uENMWIuLOFN+zEfEw8O/A7xZ8f+8BvhoRt0fE08AnCsT+tYi4I8U9DXgTMD8ifhURPwHOBf64wHnGHJukyen5+RHxZETcm953I69b7nDgzog4LyKei4iVwBLgqLJ9vhMRN6bfg19XbgOeLxDTzyPiyxGxqYHfTWuAE4UV9cAYjnmw7PtngFJfyDTgnmoHSDpU0o9SU8UI8HZgl4KvNxW4v/QgfdCsA4YKxFQZxyslXSRpWNITwKIqcdQ611S2/HndT33l+08FHkkf5OXnGKK4scT2crK/yu+veL6R1y33KmBOagYcSddzHvAbZftU+70q31YkprH8bloDnCisqMrhcU+TjZAqeWUD53oA2Ktyo6QB4BJgAfCKiBgE/oesGapaDJXWk304lc63FbAbWfNHo/6FbPTXvhGxI3BcWRz1/IIsGZbsXuCY8ve2HthF0vYV5yi9j/H87PNie5hstNurKp4fy88Psut8VUQMln3tEBEnle1T7ZqWbysSk4dutpgThY3VauCdkl4m6TeAv2rg2K8Cx0s6UNJWqeN2Jll79DbAL4FNkg4l6ycoeYjsA/SlNc57EXC4pANSM8opwJPAjxp7a0DWHv40sEHSNODvGjj2IuD9kl6dPuxPa+SFI+I+YAXwGUkvkbQf8D6gNEx5PD/7mrGlZq9L0uvuoGw47IfJ+g3G4lLgtZL+SNLk9PWm8gEG9bQgJhsDJwobq0XAHWTNAFcC3yp6YET8ADgB+AKwAbgamBYRI2QfAt8BHiNry76s7LhbgW8Da1NTxssrznsb2SieL5Mlm0OAw9OHTaNOI+sn2ED2gfftBt7fUuBLwP+SdSR/bwyvPw/Ym6wJ6RLgHyLi6vTcIsb+s68X21+QDVy4L+3zX8B5Y4ifiNgAHAwcS3Yn8yDZ3WLVkWY5mhaTjY0n3JmZWS7fUZiZWS4nCjMzy+VEYWZmuZwozMwslxOF9TRJx0m6fhzHf1dSbq2mMZ73LEn/OMZjr1GqgdXkmKYrKzk+ocVCx3uNrPVcPdYsSYX2ZkTEsaVtEfGOVrxWRPx5K87b6yRdA5wfEV9pdyz9xHcUZmaWy4nC6kplqe9RVvr7dknvqnj+BG0uDX67pDek7dMkLU7loR+VdGba/glJ55cdv0WTR2pa+bSkH6SS1ksl7SzpG8rKkt+ktCJbteaSvKYZSZ+X9EA6z0pJb03bDyGrNjsvveZPK8+VZpF/TNL9kh6WdJ5S6fGyOP5U0s+VlSb/aM7PdJGkT6fvD1BWGvxv03l/Iel9Ba9NzZjS83+SnntU2WqAa5UtAFXk3FMlXaqs7tbdkk4oe25A0n9Jejxd+79XWvWuxrlC0l9Jujf9bE5XWnujyr77p2u8If27f9r+z2QFIs9M1+jMIu/Dxs+Jwoq4h+w/6BTgk8D5ykpHIOndZBVI/4SsUuvhwKOSJpHNqr4fmE5WxK3wDGKySrB/nI7bC7iRrILqTmSzkhsqi1HmJrJKqjuRla++WNK2EXEl8BngwlSP6LeqHHtc+jqQrNz1DkDlh9VbgJlkpUc+Luk1BeN6JdnPd4isuu6XJL2swHE1Y5K0D/AfwHvJCvGVzl/UBWRFFaeSzZL/jKRSSZXTyK7rnsDvk82+ruddwGyysulHAO+v3EHSTsDlZLP2dwb+Dbhc0s4R8VHgOuCkKjWjrIWcKKyuiLg4ItanUtAXAneRlbeAbD2Az0bETak0+N0RcX96fipwSkQ8ncplN9JheW5E3JPKQHwXuCcivh/Z6mUXA7PG+F7Oj4hHU9nrfyUrJ1G09tB7gX9LZdCfAk4Fjq7o/P1kRIxGxE+BnwLVEk41G4FPRcTGiLiCbF2NInHlxXQUsDQiro+IZ4GPU7CAnrL6Vm8BPpKu3WrgK2wu7/0e4DMR8XhErCP7YK/nXyLisYj4OXAGcEyVfd4J3BURX0/X6ALgZ2TrlFibOFFYXan5YrU2l4p+HZtLbtcqGT4NuD/GvizlQ2Xfj1Z5XLU8eD2peeeO1KwxQvZX9pjKmKfvtwZeUbatUBnzKh6t+FkVPTYvpi1KikfEM2QLBwHZWullX5UVbqcCj0XEkxXnHip7vry8d5FS35XlzacWeD+Vr2tt4ERhuSS9CjgHOAnYOZX+vpXNJberlgxP23dX9aGW4ymTXe1cFDlf6o/4CNlfwy9L72UDYyxjTlbu+jm2TGITLS+mX5CVWQdeKOO+c+lxar4pff28ynl30paVesvLe29xbrYsXV5LZXnz9QXeT+XrujhdGzhRWD3bk/3n/CVA6mR9XdnzXwH+TtIblZmRksuPyT5MFkraXtK2kuakY1YDvyNp99TxeupYg4uIX5J9iBwraZKytaWrJS7ISoc/l97L1pI+TtavUvIQML1WJytZm/2HJe0haQc292mM9a6pGfJiugQ4LHUOb0PWv1RoTY2IeAD4AbAgXbvXk/WdlEqdXwScqqzU+RDZHxL1nJL2nwb8NXBhlX2uAH5TWWnyrSXNI1u2tVRF+CGyfhGbQE4Ulisibgf+lawz+SFgX7L1nEvPXwz8M1nH8JNkS13uFBGbyNqVZwA/J+sUnZeO+R7Zh8TNwErKSomP0Qlka088CryW7AOummVk/R13kjVn/Iotm0MuTv8+KuknVY7/GvB14Fqykte/Aj40ztjHq2ZMqez6h8gGEfyC7Po8TLYgUxHHkHVYrycr/X5aunYAnyK7pvcB3ydLSvXO+99k13s1WYf1Vyt3SEu3Hgr8Ldn1/Hvg0Ih4JO3yeeCoNNqqSL+INYHLjJv1iXTHMQLsnRZHaua5PwgcHRGV64qXno/0unc383VtYviOwqyHSTpM0nbKVrP7HHALsLYJ5/0NSXPSPI6ZZHcA3xnvea0zOVGY9bYjyJqO1pOtmHd0NKcZYRvgP8mas5aTNSv9RxPOax3ITU9mZpbLdxRmZpbLicLMzHL1ZJnxXXbZJaZPn97uMMzMusrKlSsfiYhdK7f3VKKQdBhw2IwZM1ixYkW7wzEz6yqSKsunAD3W9BQRSyPixClTptTf2czMCumpRGFmZs3nRGFmZrmcKMzMLFfPdmY3asmqYU5ftob1I6NMHRzglINnMneWS+CbmfXUHcVYO7OXrBrm1MW3MDwySgDDI6OcuvgWlqwarnusmVmv66lEMVanL1vD6MZNW2wb3biJ05etaVNEZmadw4kCWD8yWnX78MgocxYu952FmfU1Jwpg6uBAzefcDGVm/a6nEkWqvX/2hg0bGjrulINnMjB5Us3nRzdu4uQLV/vuwsz6Uk8lirF2Zs+dNcSCI/dlKOfOAnx3YWb9qacSxXjMnTXEDfMPqpss3MltZv3GiaJCvWYocCe3mfWXnppw1wylSXanL1vDcI3RULC5Gar8GDOzXtRTdxRj7cyuVGqGOmPefnU7ud0MZWa9rqcSRbPLjBfp5HYzlJn1up5KFK1QpJPbo6HMrJc5URRUZK6Fm6HMrBe5M7ugIp3ctUqBmJl1M99RNKBeM1ReKRAzs27lRDEG1ZqhBiZP4pSDZ7YpIjOz1umppqfxLFzUiPJmKC90ZGa9ThHR7hiabvbs2bFixYp2h2Fm1lUkrYyI2ZXb3fRkZma5nCjMzCyXE4WZmeVyojAzs1xOFGZmlqunhsd2kyWrhj281sy6ghNFGyxZNcypi29hdOMmwGtbmFln66mmp2atR9Fqpy9b80KSKHFRQTPrVD2VKJq9HkWr1Coe6LUtzKwT9VSi6BZ5xQO9toWZdRonijbw2hZm1k3cmd0GXtvCzLqJ7yjaxGtbmFm3cKJos7y1LZasGmbOwuXsMf9yd3KbWdu46anNaq1tAXiuhZl1BCeKDjB31tCLPvznLFxec66FE4WZTSQ3PXWoWp3Z7uQ2s4nmRNGhanVmu5PbzCZaTyWKbinhUYQ7uc2sU/RUouiWEh5FzJ01xIIj92VocAABQ4MDLDhyXyDr5B4eGSXwTG4zaz13Zncwd3KbWSfoqTuKfuCCgmY20ZwouowLCprZRHOi6DIuKGhmE82JosuUd3LX4mYoM2smJ4ouVK+gILgZysyax4miixVphjr5wtW+uzCzcfHw2C5WZF0LcEFBMxsf31F0uSLNUOBObjMbOyeKHlGvGQrcyW1mY+Ompx7hZigzaxXfUfSQUjPUGfP281wLM2saJ4oeVHSuhavPmlkRPZUoeqnM+HgV6eR29VkzK6KnEkUvlRlvliKd3J5vYWZ53Jnd48o7udenNSxqcUe3mVXTU3cUVl2pGeq+he/0fAsza5gTRZ/xfAsza5SbnvqM51uYWaN8R9GHPN/CzBrhO4o+VuTuotbSq2bWP3xH0efqzbfIW3rVzPqDE4UB1Tu5ByZP4pSDZ7Jk1TBzFi73TG6zPuWmJwNePN9i6uAApxw8E4BTF9/C6MZNgDu5zfqRE4W9YO6soRd9+M9ZuPyFJFFS6uR2ojDrD256sly1OrM918KsfzhRWK68zmwXFDTrD04UlqveTG7PtTDrfU4Ulqvo2hZuhjLrXU4UVleRtS3cDGXWu5worDA3Q5n1Jw+PtcJc8sOsP+XeUUiaJOn0iQrGOp9Lfpj1n9xEERGbgDdK0gTFY10ir+SHmfWWIk1Pq4D/lnQx8HRpY0QsbllU1vFqlfzwbG2z3lMkUewEPAocVLYtACeKPlet5IeZ9Z66iSIi3jcRgeSRtD1wLXBaRFzW7njMzPpJ3eGxknaT9B1JD0t6SNK3Je1W5OSSvpaOu7Vi+yGS1ki6W9L8Aqf6CHBRkdc0M7PmKtL0dC7wTeDd6fGxadvvFzh2EXAmcF5pg6RJwJfS8euAmyRdCkwCFlQc/37g9cDtwLYFXs862JJVw+7TMOtCRRLFrhFxbtnjRZJOLnLyiLhW0vSKzW8C7o6IewEkfQs4IiIWAIdWnkPSgcD2wD7AqKQrIuL5KvudCJwIsPvuuxcJzybQklXDXtfCrEsVmZn9iKRj05yKSZKOJevcHqsh4IGyx+vStqoi4qMRcTLZXc051ZJE2u/siJgdEbN33XXXcYRnrXD6sjU117Uws85WJFG8H3gP8CDwC+CotG2sqs3JiHoHRcQid2R3L69rYZ3Ey/s2JrfpKfUn/GFEHN7E11wHTCt7vBuwvonntw40dXCgZtkPN0PZRHIzaOOKzMw+osmveROwt6Q9JG0DHA1c2owTSzpM0tkbNmxoxumsiYoUFDz5wtX+685azs2gjSvS9HSDpDMlvVXSG0pfRU4u6QLgRmCmpHWSPhARzwEnAcuAO4CLIuK2Mb+DMhGxNCJOnDJlSjNOZ01UZF0LcLlya71azaAuaFlbkVFP+6d/P1W2LdhypnZVEXFMje1XAFcUeG3rIaWZ3HMWLq/ZDAWb/7pzM4C1Qq1mUBe0rK1e9ditgC9HxIEVX3WThFkt9ZqhILuzcEejtYILWjauXh/F82TNRGZNU7QZKnBTlDVf+e+fgKHBARYcua/vYHMoIn9kqqR/BEaBC9myeuxjrQ2tcZIOAw6bMWPGCXfddVe7w7ECKkeg1DI0OMAN830ja9ZKklZGxOwXbS+QKO6rsjkiYs9mBddss2fPjhUrVrQ7DCuovLRH3m/jkMt+mLVUrURRpHrsHq0JySxTXq48r6Pb493N2qNI9djtJH1M0tnp8d6SXlSTyawZisy38Hh3s4lVZB7FucCzbB4muw74dMsisr5WpKPb493NJlaRRLFXRHwW2AgQEaNUr9fUdp6Z3RvmzhrihvkH1UwWHu9uNrGKJIpnJQ2QCvdJ2gv4dUujGiPPzO4teePdXdTNbOIUmZl9GnAlME3SN4A5wHGtDMoMNndYVy52BLiom9kEqjs8FkDSzsBvkzU5/TAiHml1YOPh4bG9LW9klIfQmo3dmIfHAkTEo8DlTY/KbAzyOrN9d2HWfEX6KLqGO7P7Q73ObA+hNWuunkoU7szuD0WLCrqT26w5CiUKSW+R9L70/a6SPFvb2sZrW5hNrCIzs08DPgKcmjZNBs5vZVBm9ZTmWpwxbz/P5DZrsSJ3FO8CDidVjo2I9cBLWxmUWVFF7i7cDGU2PoUm3EU2hrY04W771oZk1ph6M7nBzVBm41EkUVwk6T+BQUknAN8HzmltWGaNc0FBs9aomygi4nPAJcC3gZnAxyPii60ObCw8PLa/uRnKrDWKLFz0YeDiiFg3MSGNn2dmW97sbchqRnn5S7Mt1ZqZXaTpaUdgmaTrJP2lpFc0Pzyz5nIzlFnzFGl6+mREvBb4S2Aq8L+Svt/yyMzGwetamDVPoVpPycPAg8CjwMtbE45Z85SWWK3VDDV1cGCL9bqnuqCgWVVFJtx9UNI1wFXALsAJEfH6Vgdm1iy11rU48NW7curiWxgeGSXwEFqzWor0UbwKODkiXhsRp0XE7a0OyqyZypuhRFaKfMGR+3L1z375wpoWJe67MHuxmk1PknaMiCeAz6bHO5U/HxGPtTg2s6YpNUOV+/CFq6vu674Lsy3l9VF8EzgUWEk2K7t8newA9mxhXGYtN3VwoGbfhZltVrPpKSIOTf/uERF7pn9LXx2ZJDzhzhqRtya3mW1WpDP7qiLbOoHXo7BG1Oq78Kgnsy3l9VFsC2wH7CLpZWxuetqRbD6FWder1ndhZlvK66P4M+BksqSwks2J4gngSy2Oy6ytPL/CbLOaiSIiPg98XtKHOrUIoFkrLFk1zKmLb3lh6GxpfgXgZGF9qe7M7Ij4oqTXAfsA25ZtP6+VgZm1y+nL1tScX+FEYf2obqJIS6EeQJYorgDeAVwPOFFYT6o1j2L9yKibpKwvFZmZfRTwe8CDEfE+4LeAl7Q0KrM2qjWPYsrAZJf8sL5UJFGMRsTzwHOSdiQrDtiR8yjMmqHW/AoJl/ywvlQkUayQNEi2/OlK4CfAj1salVkb1ZpfMfLMxqr7D4+Mssf8y71ynvWsuivcbbGzNB3YMSJublVAzeAV7qwV6q2aB145z7pbwyvcSXpD5RewE7B1+r7juISHtVK9VfMga4o6+cLVvruwnlLzjkLS1TnHRUQc1JqQxs93FNYq5aOe6t2L++7Cuk2tO4qGmp66hROFTYQiTVGQ9XF4GK11g4abnsoO3E7SxySdnR7vLenQVgRp1k2KNEWBh9Fa9ysy6ulc4Flg//R4HfDplkVk1iXKR0fV42G01s2KJIq9IuKzwEaAiBhly0WMzPrW3FlD3DD/IM6Yt1/duwuvnGfdqm4JD+BZSQNkq9ohaS/g1y2NyqzLlPofTl+2pma/hVfOs25V5I7iNOBKYJqkbwBXAX/f0qjMulDe3YVXzrNulntHIUnAz4Ajgd8ma3L664h4ZAJiM+tK5XcXlcUDXVTQulHd4bFpuNQbJyiepvDwWOtEletcgOdaWGcZ8/BY4IeS/l8LYjLrK3nrXJh1siKJ4kDgRkn3SLpZ0i2SOrrWk1knqjXqaXhk1CU/rKMVGfX0jpZHYdYHpg4O1BwR5eVWrZPVvaOIiPurfU1EcGa9pN5MbjdDWacqckdhZk1QZK6FJ+VZJyrSR9E1XGbcOl1prkWtsh+elGedqKcSRUQsjYgTp0yZ0u5QzHLVWm7Vk/KsE7npyawN8iblmXUaJwqzNpk7a8iJwbpCTzU9mZlZ8zlRmJlZLjc9mXWw8iKCUwYmI8HIMxvdp2ETyonCrENVFhEcGd34wnOeyW0TyU1PZh2qWhHBcp7JbRPFicKsQxWZpe2Z3DYRnCjMOlSRWdqeyW0TwYnCrEPVKyLomdw2UdyZbdahKmdve9STtYsThVkHKzJ72+twW6vVXTO7G3nNbOsX1dbhFhDAkJOGNWg8a2abWYeqNoS29Kdfaa6Fl1i18XKiMOti9YbHeq6FNYMThVkXKzI81nMtbLycKMy6WL0htLA5mSxZNcychcvZY/7lzFm43E1SVphHPZl1scp1uEsd2SWluRaVnd6uFWWNcKIw63LlQ2hrDZWds3D5izq9S/0XThRWjxOFWQ+pNe+iVj+F+y+sCPdRmPWBWp3erhVlRXR8opB0gKTrJJ0l6YB2x2PWjap1ertWlBXV0kQh6WuSHpZ0a8X2QyStkXS3pPl1ThPAU8C2wLpWxWrWy+bOGmLBkfsyNDiAgMGByWw7eSs+fOFqj4CyulpawkPS75B9yJ8XEa9L2yYBdwK/T/bBfxNwDDAJWFBxivcDj0TE85JeAfxbRLy33uu6hIdZbdXKfgxMnsSCI/d1x3afq1XCo6Wd2RFxraTpFZvfBNwdEfemwL4FHBERC4BDc073OPCSWk9KOhE4EWD33XcfR9Rmva1a2Q+PgLI87eijGAIeKHu8Lm2rStKRkv4T+DpwZq39IuLsiJgdEbN33XXXpgVr1mtqjXQaHhl1M1QT9dIEx3YMj1WVbTXbvyJiMbC4deGY9ZepgwMM5yQLT8Qbv7FMcOzkcvHtuKNYB0wre7wbsL4NcZj1pXplP1xIcPzymveqKSWW4ZFRgs6r/NuORHETsLekPSRtAxwNXNqME0s6TNLZGzZsaMbpzHpS+QioWtwMNT6NTnBsNLFMtFYPj70AuBGYKWmdpA9ExHPAScAy4A7gooi4rRmvFxFLI+LEKVOmNON0Zj1r7qwhbph/UN1k0Ul/1XaTRic4dvrM+VaPejqmxvYrgCta+dpmVt8pB8980VDZcqMbN3Hyhav5xKW3eb3uBlT7ueZNcKzVb9QpM+c7fma2mbVOkWYogJHRjTz+zMaObD/vRJUTHIcGB3LnqXT6zHmvmW1mAMxZuLzmaKhqhgYHuGH+QS2MqL90wqintky4m2iSDgMOmzFjRrtDMes69ZqhKnVK+3mvqFX5txP0VNOTO7PNxq5oM1RJgEdG9YmeShRmNj6l0VBnzNuv7hKr4P6KfuFEYWYvUq3a7Mu2m1x139LIKN9d9K6e6qMws+ap1ma+x/zLa9bbcfmP3tVTdxSemW3WWvXG9XfSbGJrnp5KFO7MNmutenWiwKOhepGbnsyssFKT0unL1tScc9Eps4mteXrqjsLMWi9vZFQnzSa25vEdhZmNSfndRSeuoWDN01OJwjOzzSZWrdnEnVCOwprHtZ7MrKkqV3eDbFnLIKsP5aTRuWrVenIfhZk1VbVFeEp/jnomd3dyojCzpqo3PNZzLbqPE4WZNVWR4bFearW7OFGYWVMVmZQHbobqJk4UZtZUleXKlbOvm6G6g4fHmlnTlQ+bLQ2VrTWT2yU/Ol9P3VG41pNZ5ynN5K61IFK3lPxYsmqYOQuXs8f8y/uuf6WnEoWZda5qfRfdUvKjNDdkeGSUoP/6V5wozGxCVC6GNDQ4wIIj9+2KyXfV5ob0U/9KT/VRmFlnq1XyAzq77EetfpR+6V/xHYWZtV2nN+3U6kfplv6V8XKiMLO26/SmnW7uX2kGNz2ZWdvlNe2UN0lNGZiMBCPPbJzQ5ql+L6nuRGFmbTd1cKDqPIspA5O3qEQ7MrrxhedKzVPAhCWLfkkMlXqq6UnSYZLO3rBhQ7tDMbMG1GrakXhRk1S5Tmqe6mU9lSg84c6sO9UaOjvyzMa6x/bLyKN2ctOTmXWEak07eaU/Svpl5FE79dQdhZn1lnqVaPtp5FE7+Y7CzDpW5Wijdo166ndOFGbW0fp5tFGncNOTmZnlcqIwM7NcThRmZpbLicLMzHI5UZiZWa6eGvVUWjMbeELSXRVPTwEqa3tU27YL8EhrIqyrWjwTcZ6i+9fbL+/5Ws91+nVp1zUpesx49unWawLNuS6tuiZF9mvV/5XxXpNXVd0aEX3xBZxdcNuKTopxIs5TdP96++U9X+u5Tr8u7bomRY8Zzz7dekSpMKEAAAb3SURBVE2adV1adU2K7Neq/yutuib91PS0tOC2dmpWPI2ep+j+9fbLe77Wc51+Xdp1TYoeM559uvWaQHPiadU1KbJfV/1fUcpClkhaERGz2x2HbcnXpfP4mnSeVl2TfrqjKOrsdgdgVfm6dB5fk87TkmviOwozM8vlOwozM8vlRGFmZrmcKMzMLJcTRR2S5ko6R9J/S3p7u+MxkPQaSWdJukTSB9sdj20maXtJKyUd2u5YDCQdIOm69P/lgLGepy8ThaSvSXpY0q0V2w+RtEbS3ZLmA0TEkog4ATgOmNeGcPtCg9fkjoj4c+A9gIdntlAj1yX5CHDRxEbZXxq8JgE8BWwLrBvra/ZlogAWAYeUb5A0CfgS8A5gH+AYSfuU7fKx9Ly1xiIauCaSDgeuB66a2DD7ziIKXhdJbwNuBx6a6CD7zCKK/1+5LiLeQZbAPznWF+zLRBER1wKPVWx+E3B3RNwbEc8C3wKOUOZfgO9GxE8mOtZ+0cg1SftfGhH7A++d2Ej7S4PX5UDgt4E/Ak6Q1JefL63WyDWJiOfT848DLxnra/ZUUcBxGgIeKHu8Dngz8CHgbcAUSTMi4qx2BNenql6T1NZ6JNkv/hVtiKvfVb0uEXESgKTjgEfKPqSs9Wr9XzkSOBgYBM4c68mdKDZTlW0REV8AvjDRwRhQ+5pcA1wzsaFYmarX5YVvIhZNXCiW1Pq/shhYPN6T+9Zws3XAtLLHuwHr2xSLZXxNOpOvS+dp6TVxotjsJmBvSXtI2gY4Gri0zTH1O1+TzuTr0nlaek36MlFIugC4EZgpaZ2kD0TEc8BJwDLgDuCiiLitnXH2E1+TzuTr0nnacU1cFNDMzHL15R2FmZkV50RhZma5nCjMzCyXE4WZmeVyojAzs1xOFGZmlsuJwvqCpGsktbwkuaS/knSHpG9UbN9P0h+UPT68ojx3s+OYW1H9uBWvsVbSt8seHyVpUStf09rDicKsDkmN1ET7C+APIqKyqu1+wAuJIlW/XdiM+GqYS1ZuutVmS3rtBLyOtZEThXUMSdPTX+PnSLpN0v9IGkjPvXBHIGkXSWvT98dJWiJpqaT7JJ0k6W8krZL0Q0k7lb3EsZJ+IOlWSW9Kx2+fFoK5KR1zRNl5L5a0FPifKrH+TTrPrZJOTtvOAvYELpX04bJ9twE+BcyTtFrSvHT+M9PziyR9WdLVku6V9LsppjvK/0KX9HZJN0r6SYpth7R9oaTbJd0s6XOS9gcOB05Pr7dX+rpS2epz10l6ddlrn5W23am0Mp2k10r6cTr+Zkl717hsnwP+oZHrbF0oIvzlr474AqYDzwH7pccXAcem768BZqfvdwHWpu+PA+4GXgrsCmwA/jw99+/AyWXHn5O+/x3g1vT9Z8peYxC4E9g+nXcdsFOVON8I3JL22wG4DZiVnlsL7FLlmOOAM6s9JluI5ltkFUCPAJ4A9iX7Q24l2d3ILsC1wPbpmI8AHwd2AtawucrCYNk5jyp7vauAvdP3bwaWl+13ZXqtvdN73hb4IvDetM82wECV97QWeAVZyYgZwFHAonb/Hvmr+V8uM26d5r6IWJ2+X0mWPOq5OiKeBJ6UtAFYmrbfAry+bL8LIFv4RdKOkgaBtwOHS/q7tM+2wO7p++9FROUCMQBvAb4TEU8DSFoMvBVYVeQN1rA0IkLSLcBDEXFLOvdtZD+D3ciakm6QBNmH941kSeVXwFckXQ5cVnnidOexP3BxOha2XMTmosjWjrhL0r3Aq9O5PyppN2BxRNxVI+5NwOnAqcB3x/jercM5UVin+XXZ95uAgfT9c2xuKt0255jnyx4/z5a/45WFzYLsr/g/jIg15U9IejPwdI0Yq9X+H6/ymCvfz9ZkP4vvRcQxLwoma0b7PbKKoScBB1XsshUwEhH71XjtF/1cIuKbkn4EvBNYJun4iFhe4/ivkyUKFwbsUe6jsG6xlqzJB7ImjrGYByDpLcCGiNhAVm3zQ0p/akuaVeA81wJzJW0naXvgXcB1dY55kqx5bKx+CMyRNCPFuZ2k30x3C1Mi4grgZLJmqi1eLyKeAO6T9O50rCT9Vtm53y1pK0l7kfWxrJG0J3BvZAt3XcqWd2ZbiIiNpGa+cbw/62BOFNYtPgd8UNIPyNrrx+LxdPxZwAfStn8CJgM3S7o1Pc4V2drpi4AfAz8CvhIR9Zqdrgb2KXVmNxp4RPySrF/jAkk3kyWOV5Mlg8vStv8FSp3o3wJOSR30e5GtLf4BST8l+8v/iLLTr0nHfpesf+dXZEn1Vkmr0+ucVyfEr+IWip7lMuNmfSyNqrosIi5pdyzWuXxHYWZmuXxHYWZmuXxHYWZmuZwozMwslxOFmZnlcqIwM7NcThRmZpbLicLMzHL9HxFBxpFJKiBxAAAAAElFTkSuQmCC\n", | |
| "text/plain": [ | |
| "<Figure size 432x288 with 1 Axes>" | |
| ] | |
| }, | |
| "metadata": { | |
| "needs_background": "light" | |
| }, | |
| "output_type": "display_data" | |
| } | |
| ], | |
| "source": [ | |
| "plt.loglog(n,error,'o')\n", | |
| "plt.xlabel('number of timesteps N')\n", | |
| "plt.ylabel('relative error')\n", | |
| "plt.title('Truncation and roundoff error \\naccumulation in log-log plot')" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "In the above plot \"Truncation and roundoff error accumulation in log-log plot\", we see that around N=20,000 steps we stop decreasing the error with more steps. This is because we are approaching the limit of how precise we can store a number using a 32-bit floating point number. \n", | |
| "\n", | |
| "In any computational solution, there will be some point of similar diminishing in terms of accuracy (error) and computational time (in this case, number of timesteps). If you were to attempt a solution for N=1 billion, the solution could take $\\approx$(1 billion)*(200 $\\mu$s\\[cpu time for n=2\\])$\\approx$55 hours, but would not increase the accuracy of the solution. " | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "## What we've learned\n", | |
| "\n", | |
| "* Numerical integration with the Euler approximation\n", | |
| "* The source of truncation errors\n", | |
| "* The source of roundoff errors\n", | |
| "* How to time a numerical solution or a function\n", | |
| "* How to compare solutions\n", | |
| "* The definition of absolute error and relative error\n", | |
| "* How a numerical solution converges" | |
| ] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "# Problems\n", | |
| "\n", | |
| "1. The growth of populations of organisms has many engineering and scientific applications. One of the simplest\n", | |
| "models assumes that the rate of change of the population p is proportional to the existing population at any time t:\n", | |
| "\n", | |
| " $\\frac{dp}{dt} = k_g p$\n", | |
| "\n", | |
| " where $t$ is time in years, and $k_g$ is growth rate in \\[1/years\\]. \n", | |
| " \n", | |
| " The world population has been increasing dramatically, let's make a prediction based upon the [following data](https://worldpopulationhistory.org/map/2020/mercator/1/0/25/) saved in [world_population_1900-2020.csv](../data/world_population_1900-2020.csv):\n", | |
| " \n", | |
| " |year| world population |\n", | |
| " |---|---|\n", | |
| " |1900|1,578,000,000|\n", | |
| " |1950|2,526,000,000|\n", | |
| " |2000|6,127,000,000|\n", | |
| " |2020|7,795,482,000|\n", | |
| " \n", | |
| " a. Calculate the average population growth, $\\frac{\\Delta p}{\\Delta t}$, from 1900-1950, 1950-2000, and 2000-2020\n", | |
| " \n", | |
| " b. Determine the average growth rates. $k_g$, from 1900-1950, 1950-2000, and 2000-2020\n", | |
| " \n", | |
| " c. Use a growth rate of $k_g=0.013$ [1/years] and compare the analytical solution (use initial condition p(1900) = 1578000000) to the Euler integration for time steps of 20 years from 1900 to 2020 (Hint: use method (1)- plot the two solutions together with the given data) \n", | |
| " \n", | |
| " d. Discussion question: If you decrease the time steps further and the solution converges, will it converge to the actual world population? Why or why not? \n", | |
| "\n", | |
| "**Note: We have used a new function `np.loadtxt` here. Use the `help` or `?` to learn about what this function does and how the arguments can change the output. In the next module, we will go into more details on how to load data, plot data, and present trends.**" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 12, | |
| "metadata": {}, | |
| "outputs": [ | |
| { | |
| "name": "stdout", | |
| "output_type": "stream", | |
| "text": [ | |
| "years= [1900. 1950. 2000. 2020.]\n", | |
| "population = [1.578000e+09 2.526000e+09 6.127000e+09 7.795482e+09]\n" | |
| ] | |
| } | |
| ], | |
| "source": [ | |
| "import numpy as np\n", | |
| "year,pop=np.loadtxt('../data/world_population_1900-2020.csv',skiprows=1,delimiter=',',unpack=True)\n", | |
| "print('years=',year)\n", | |
| "print('population =', pop)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| }, | |
| { | |
| "cell_type": "markdown", | |
| "metadata": {}, | |
| "source": [ | |
| "2. In the freefall example we used smaller time steps to decrease the **truncation error** in our Euler approximation. Another way to decrease approximation error is to continue expanding the Taylor series. Consider the function f(x)\n", | |
| "\n", | |
| " $f(x)=e^x = 1+x+\\frac{x^2}{2!}+\\frac{x^3}{3!}+\\frac{x^4}{4!}+...$\n", | |
| "\n", | |
| " We can approximate $e^x$ as $1+x$ (first order), $1+x+x^2/2$ (second order), and so on each higher order results in smaller error. \n", | |
| " \n", | |
| " a. Use the given `exptaylor` function to approximate the value of exp(1) with a second-order Taylor series expansion. What is the relative error compared to `np.exp(1)`?\n", | |
| " \n", | |
| " b. Time the solution for a second-order Taylor series and a tenth-order Taylor series. How long would a 100,000-order series take (approximate this, you don't have to run it)\n", | |
| " \n", | |
| " c. Plot the relative error as a function of the Taylor series expansion order from first order upwards. (Hint: use method (4) in the comparison methods from the \"Truncation and roundoff error accumulation in log-log plot\" figure)" | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": 1, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [ | |
| "from math import factorial\n", | |
| "def exptaylor(x,n):\n", | |
| " '''Taylor series expansion about x=0 for the function e^x\n", | |
| " the full expansion follows the function\n", | |
| " e^x = 1+ x + x**2/2! + x**3/3! + x**4/4! + x**5/5! +...'''\n", | |
| " if n<1:\n", | |
| " print('lowest order expansion is 0 where e^x = 1')\n", | |
| " return 1\n", | |
| " else:\n", | |
| " ex = 1+x # define the first-order taylor series result\n", | |
| " for i in range(1,n):\n", | |
| " ex+=x**(i+1)/factorial(i+1) # add the nth-order result for each step in loop\n", | |
| " return ex\n", | |
| " " | |
| ] | |
| }, | |
| { | |
| "cell_type": "code", | |
| "execution_count": null, | |
| "metadata": {}, | |
| "outputs": [], | |
| "source": [] | |
| } | |
| ], | |
| "metadata": { | |
| "celltoolbar": "Slideshow", | |
| "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.5" | |
| } | |
| }, | |
| "nbformat": 4, | |
| "nbformat_minor": 4 | |
| } |