[1]_Module_1_Integrals_and_Derivatives.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G0BcK0H_G1Iu"
      },
      "source": [
        "# **AP155 Lab Assignment** \n",
        "## Module 1: Integrals and Derivatives\n",
        "\n",
        "_Instructions_: Answer each problem as completely as you can. Discuss **all** your answers as clearly and concisely as possible. \n",
        "\n",
        "_Scoring Criteria_: 50% - *correctness of code*; 50% - *discussion of the code and results*. Maximum score is **100 points** (Problems 1 to 5 are worth 20, 20, 10, 30, and 20 points, respectively).\n",
        "\n",
        "_Due Date_: **5PM of November 3 (Thu)**\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "PeHZ9kP8JaVR"
      },
      "source": [
        "### Student Information \n",
        "\n",
        "_Full Name (Last Name, First Name)_: Jonabel Eleanor B. Baldres \\\n",
        "_Student No._: 2020-07279 \\\n",
        "_Section_: THU-HX-3"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Hc183co-J4eI"
      },
      "source": [
        "### Submission Information\n",
        "\n",
        "_Date and Time Submitted (most recent upload)_: 27 November 2022 \n",
        "\n",
        "**HONOR PLEDGE** I affirm that I have upheld the highest principles of honesty and integrity in my academic work and that this lab assignment is my own work.\n",
        "\n",
        "**Sign here with your full name:** sgd. Jonabel Eleanor Baldres"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "X1UGsRpW-eQs"
      },
      "source": [
        "### Grading Information (c/o Lab Instructor)\n",
        "\n",
        "TOTAL SCORE: **96**/100\n",
        "\n",
        "Score breakdown:\n",
        "* Problem 1 - 20/20\n",
        "* Problem 2 - 20/20\n",
        "* Problem 3 - 8/10\n",
        "* Problem 4 - 28/30\n",
        "* Problem 5 - 20/20\n",
        "### The notebook is very neat. Keep up the good work! -Sir Rene \n",
        "_Date and Time Scored (MM/DD/YYYY HH:MM AM/PM):_ "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "hMqWOsbIMoUh"
      },
      "source": [
        "# importing all necessary functions in the problems\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "from math import sin, sqrt, tanh, cosh, inf\n"
      ],
      "execution_count": 972,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LxWuqiCFMsEh"
      },
      "source": [
        "### PROBLEM 1 (20 pts.) \n",
        "**Integration using the trapezoidal rule and Simpson's rule**\n",
        "\n",
        "Calculate the integral of the function $ 2x^4 + 3x^2 + 4x + 5$ from $x=0$ to $x=5$.\n",
        "\n",
        "1. Do the numerical calculation using the trapezoidal rule with 10, 100, and 1,000 slices. Describe each step of the calculation clearly and concisely. *(8 pts.; 4 pts. for the code & 4 pts. for the discussion text)*\n",
        "\n",
        "2. What are the percentage errors of your answers in (1)? *(2 pts.)*\n",
        "\n",
        "3. Repeat the calculation using Simpson's rule with 10, 100, and 1,000 slices. Describe each step of the calculation clearly and concisely. *(8 pts.;  4 pts. for the code & 4 pts. for the discussion text)*\n",
        "\n",
        "4. What are the percentage errors of your answers in (3)? How do they compare with the results from (1) for the same number of slices? *(2 pts.)*"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def f(x):\n",
        "  ''' defining the function provided that needs to be integrated'''\n",
        "  return 2*x**4 + 3*x**2 + 4*x + 5"
      ],
      "metadata": {
        "id": "8pAti1JPGzsy"
      },
      "execution_count": 973,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def int_f(x):\n",
        "  ''' defining the integral of the function '''\n",
        "  return 2*x**5/5 + x**3 + 2*x**2 + 5*x\n"
      ],
      "metadata": {
        "id": "ySEdfabVG7kZ"
      },
      "execution_count": 974,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "print(\"The actual value of the integral is \" + str(int_f(5))) \n",
        "#since the value of the function is 0 when x = 0, the author decided to just find the integral when x = 5\n",
        "# this will be used in obtaining the percentage errors in the next steps"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "8LcBINg9HTIc",
        "outputId": "4dfde00a-9930-494e-fbc8-6f4745eb2c06"
      },
      "execution_count": 975,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The actual value of the integral is 1450.0\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Using the **Trapezoidal Rule** \\\n",
        "\n",
        "Here, we created a function trapezoidal to calculate the integral of the given function. We create an argument *N* which is the variable for the number of steps needed in the numerical calculation using the Trapezoidal Rule. Note that the Trapezoidal rule is given by the equation\n",
        "\\begin{align}\n",
        "        I(a,b) = h \\left[ \\frac{1}{2} f(a) + \\frac{1}{2} f(b) + \\sum_{k=1}^{N-1} f(a+kh) \\right]\n",
        "    \\end{align}\n",
        "where h is the width of the slices, a is the lower limit of integration, b is the upper limit of integration, k is the position of slice, and $f(n)$ is the function integrated. "
      ],
      "metadata": {
        "id": "MiJzqAqiHy4h"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "a = 0 # lower limit of integration\n",
        "b = 5 # upper limit of integration\n",
        "\n",
        "def trapezoidal(N):\n",
        "  ''' this function provides the integral using the trapezoidal method'''\n",
        "  h = (b-a)/N # width of the slices \n",
        "  s = 0.5*f(a) + 0.5*f(b) \n",
        "\n",
        "  for k in range(1,N): \n",
        "    s+= f(a+k*h) # this adds a value to the variable s and the new value becomes the value of s\n",
        "  return h*s # returning h*s means the functions shows the product of s with h that results to the integral of the given function given N"
      ],
      "metadata": {
        "id": "f6NfSN4yHx2y"
      },
      "execution_count": 976,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "for i in [10,100,1000]: # used for loop to yield the integral of the function when N= 10,100, and 1000\n",
        "  print(\"When N=\" + str(i) +  \",the integral is \" + str(trapezoidal(i)))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "jmZ-InuaIRpv",
        "outputId": "75ea5b10-9505-495f-e65b-e1fa99c4a0e5"
      },
      "execution_count": 977,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "When N=10,the integral is 1471.4375\n",
            "When N=100,the integral is 1450.2145812500003\n",
            "When N=1000,the integral is 1450.0021458331234\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def trapezoidal_error(N):\n",
        "  '''calculates the numerical error produced by the trapezoidal rule method'''\n",
        "  return (trapezoidal(N) - int_f(5)) / int_f(5) * 100"
      ],
      "metadata": {
        "id": "etfbe2kbJp1Y"
      },
      "execution_count": 978,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "for k in [10,100,1000]: # shows the error when the integral has N number of steps\n",
        "  print(\"The error when the integral has \" + str(k) + \" steps is \" + str(trapezoidal_error(k)) )"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "95QtATXGLJEK",
        "outputId": "0c151d59-0927-493d-a84d-d467793c2922"
      },
      "execution_count": 979,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The error when the integral has 10 steps is 1.478448275862069\n",
            "The error when the integral has 100 steps is 0.014798706896569988\n",
            "The error when the integral has 1000 steps is 0.00014798849127240326\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Using the **Simpsons Rule** \\\n",
        "In this problem, we define the function simpsons with argument N for the number of slices to produce the integral of a given function. We note that the Simpson's rule is given by the equation:\n",
        "\\begin{align}\n",
        "        \\frac{1}{3} h \\left[ f(a) + f(b) + 4 \\sum_{k \\ odd} f(a+kh) + 2 \\sum_{k \\ even} f(a+kh) \\right]\n",
        "    \\end{align}\n",
        "where h is the width of the slices, a and b is the lower and upper limit of integration respectively, k is the position of the slice, and *f(n)* is the function"
      ],
      "metadata": {
        "id": "aqCML6ISMO4-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#using the same lower and upper integration bounds a and b\n",
        "\n",
        "def simpsons(N):\n",
        "  ''' the function is used to create the code for the Simpson's Rule'''\n",
        "  h = (b-a) / N # width of the slice\n",
        "  c = h/3 \n",
        "  s = f(a) + f(b) # first two equation that needs to be added together\n",
        "\n",
        "  for k in range(2,N,2): # when the position of the slice is even, we add the equation below to the initial value of s and the sum becomes the new s\n",
        "    s += 2 * f(a+k*h)\n",
        "\n",
        "  for k in range(1,N,2): # when the position of the slice is odd, we add the equation below to the initial value of s and the sum becomes the new s\n",
        "    s += 4 * f(a+k*h) \n",
        "\n",
        "  return c*s #calling the function returns the product of c and s which is the value of the integral using the Simpson's rule\n"
      ],
      "metadata": {
        "id": "KXgRYXb4MSRJ"
      },
      "execution_count": 980,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "for k in [10,100,1000]: # used for loop to yield the integral of the function when N= 10,100, and 1000\n",
        "  print(\"When N=\" + str(k) +  \",the integral is \" + str(simpsons(k)))"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "4ljQrffAMy0o",
        "outputId": "75731338-6669-4f52-fd60-14a28b3f2425"
      },
      "execution_count": 981,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "When N=10,the integral is 1450.0833333333333\n",
            "When N=100,the integral is 1450.0000083333334\n",
            "When N=1000,the integral is 1450.0000000008336\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def simpsons_error(N): # returns the value of the percentage error when Simpsons rule is used\n",
        "  return (simpsons(N) - int_f(5)) / int_f(5) * 100"
      ],
      "metadata": {
        "id": "i_PVN2VzNMsE"
      },
      "execution_count": 982,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "for l in [10,100,1000]:\n",
        "  print(\"The error when the integral has \" + str(l) + \" steps is \" + str(simpsons_error(l)) )"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "3_XtuqZmNSZa",
        "outputId": "f8a01a3d-fde0-4c92-93ff-9284950c0a42"
      },
      "execution_count": 983,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The error when the integral has 10 steps is 0.005747126436776382\n",
            "The error when the integral has 100 steps is 5.747126500676226e-07\n",
            "The error when the integral has 1000 steps is 5.748633752930267e-11\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "In summary, the percentage errors when Simpsons Rule was used were smaller compared to Trapezoidal rule. This means that Simpsons Rule produces more accurate integral value than Trapezoidal rule. "
      ],
      "metadata": {
        "id": "-XJDuH3HOeSl"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kBX9EyMLOnXj"
      },
      "source": [
        "### PROBLEM 2 (20 pts.)\n",
        "\n",
        "**Integration using the adaptive trapezoidal rule method and Romberg integration**\n",
        "\n",
        "Calculate the integral $ I = \\int_0^1 \\sin^2 \\sqrt{100x} \\mbox{d}x $.\n",
        "\n",
        "1. Do the numerical calculation using the adaptive trapezoidal rule method to an approximate accuracy of $\\epsilon = 10^{-6}$. Start with one single integration slice and work up from there to two, four, eight, and so forth. Save and print out the number of slices, its estimate of the integral, and its estimate of the error on the integral, for each value of the number of slices $N$, until the target accuracy is reached. *(8 pts.; 4 pts. for the code & 4 pts. for the discussion text)*\n",
        "\n",
        "2. Repeat the calculation using the Romberg integration technique. Save and print out the same variables as in (1). *(8 pts.; 4 pts. for the code & 4 pts. for the discussion text)*\n",
        "\n",
        "3. Compare the results from both calculations by creating a plot of error estimates vs. $N$. Annotate your chart completely (include proper axis labels and legend). Describe how the trends for the two techniques compare. *(4 pts.; 2 pts. for the chart & 2 pts. for the discussion text)*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "GABjwlXCR97i"
      },
      "source": [
        "a = 0 # lower integration bound\n",
        "b =1 #upper integration bound\n",
        "data = [] # to store the values "
      ],
      "execution_count": 984,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "Using the **Adaptive Trapezoidal Method** \\\n",
        "Here, we take note of the following equation\n",
        "\n",
        "\\begin{align}\n",
        "    I_i = \\frac{1}{2} I_{i - 1} + h_i \\sum_{k \\ odd} f(a+kh_i)\n",
        "  \\end{align}\n",
        "\n",
        "which is the adaptive trapezoidal method where $I_i$ is integral of the ith step of the calculation, $h_i$ is the width of the slice, and $k$ is the position of the slice. "
      ],
      "metadata": {
        "id": "aGyqZz2yEFFd"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def g(x):\n",
        "  '''returns the function'''\n",
        "  return (sin(10*sqrt(x)) ** 2)\n",
        "\n",
        "def integrate(function, N):\n",
        "  h = (b-a)/N # width of the slice\n",
        "  s = 0 #starting point\n",
        "  for i in range(1,N,2): #when i is odd, we add the equation below to the initial value of s and the sum becomes the new s\n",
        "    s += function(a+i*h) \n",
        "\n",
        "  return s * h #calling the function returns the product of s and h\n",
        "\n",
        "def adaptive_trapezoidal(function):\n",
        "  N= 1 #starting number of slice\n",
        "  initial = 0 # initial value of integral\n",
        "  h = (b-a) / N\n",
        "  final = h * (function(a) + function(b)) \n",
        "  error= abs(final - initial)/3 #calculates the error \n",
        "\n",
        "  while error > 1.e-6: #while error > 1.e-6, the program continues printing the integral until the desired value of error is achieved\n",
        "    initial = final\n",
        "    final = (0.5 * initial) + integrate(function, N)\n",
        "    error = abs(final - initial)/3\n",
        "    N*= 2 # doubles the number of steps \n",
        "    data.append([N/2, error, final]) #adds the value of final, error, and N/2 to the list to store the value\n",
        "\n",
        "  return data\n",
        "\n",
        "data = adaptive_trapezoidal(g)"
      ],
      "metadata": {
        "id": "TeKwVpFzPX3h"
      },
      "execution_count": 985,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "The table below shows the number of slices (4096) it takes to approach the desired accuracy and the value of its integral.  From the program, it stated it took 1 second to fully show the output. "
      ],
      "metadata": {
        "id": "tji3Wv4bPzpm"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "df = pd.DataFrame(data, columns =['Number of slices', 'Error', 'Value of Integral'])\n",
        "df.index += 1\n",
        "df"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 452
        },
        "id": "ic8kCLIeMkm0",
        "outputId": "b835c4e4-544f-4edc-c469-e6206e182426"
      },
      "execution_count": 986,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "    Number of slices         Error  Value of Integral\n",
              "1                1.0  4.932649e-02           0.147979\n",
              "2                2.0  5.908414e-02           0.325232\n",
              "3                4.0  6.235031e-02           0.512283\n",
              "4                8.0  3.642847e-02           0.402997\n",
              "5               16.0  9.035307e-03           0.430103\n",
              "6               32.0  6.103765e-03           0.448415\n",
              "7               64.0  1.832755e-03           0.453913\n",
              "8              128.0  4.785244e-04           0.455349\n",
              "9              256.0  1.209207e-04           0.455711\n",
              "10             512.0  3.031107e-05           0.455802\n",
              "11            1024.0  7.582827e-06           0.455825\n",
              "12            2048.0  1.896023e-06           0.455831\n",
              "13            4096.0  4.740255e-07           0.455832"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-cbb4e8c5-5f99-4046-8133-fa97fe9fecb0\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Number of slices</th>\n",
              "      <th>Error</th>\n",
              "      <th>Value of Integral</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1.0</td>\n",
              "      <td>4.932649e-02</td>\n",
              "      <td>0.147979</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>2.0</td>\n",
              "      <td>5.908414e-02</td>\n",
              "      <td>0.325232</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>4.0</td>\n",
              "      <td>6.235031e-02</td>\n",
              "      <td>0.512283</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>8.0</td>\n",
              "      <td>3.642847e-02</td>\n",
              "      <td>0.402997</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>16.0</td>\n",
              "      <td>9.035307e-03</td>\n",
              "      <td>0.430103</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>32.0</td>\n",
              "      <td>6.103765e-03</td>\n",
              "      <td>0.448415</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>64.0</td>\n",
              "      <td>1.832755e-03</td>\n",
              "      <td>0.453913</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>8</th>\n",
              "      <td>128.0</td>\n",
              "      <td>4.785244e-04</td>\n",
              "      <td>0.455349</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>9</th>\n",
              "      <td>256.0</td>\n",
              "      <td>1.209207e-04</td>\n",
              "      <td>0.455711</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>10</th>\n",
              "      <td>512.0</td>\n",
              "      <td>3.031107e-05</td>\n",
              "      <td>0.455802</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>11</th>\n",
              "      <td>1024.0</td>\n",
              "      <td>7.582827e-06</td>\n",
              "      <td>0.455825</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>12</th>\n",
              "      <td>2048.0</td>\n",
              "      <td>1.896023e-06</td>\n",
              "      <td>0.455831</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>13</th>\n",
              "      <td>4096.0</td>\n",
              "      <td>4.740255e-07</td>\n",
              "      <td>0.455832</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-cbb4e8c5-5f99-4046-8133-fa97fe9fecb0')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-cbb4e8c5-5f99-4046-8133-fa97fe9fecb0 button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-cbb4e8c5-5f99-4046-8133-fa97fe9fecb0');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ]
          },
          "metadata": {},
          "execution_count": 986
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Using the **Romberg Integration** \\\n",
        "The Romberg integration is given by the formula\n",
        "\n",
        "\\begin{align}\n",
        "   I = R_{i,m+1} + \\text{error}\n",
        "  \\end{align}\n",
        "where\n",
        "\\begin{align}\n",
        "   R_{i,m+1} = R_{i,m} + \\frac{1}{4^m -1 } (R_{i,m} - R_{i-1,m})\n",
        "  \\end{align}\n",
        "and\n",
        "\\begin{align}\n",
        "  \\text{error} = \\frac{1}{4^m -1} (R_{i,m} - R_{i-1,m})\n",
        "  \\end{align}"
      ],
      "metadata": {
        "id": "wwW8IA5RQtVR"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "#adding important constants\n",
        "N = 1 # initial number of slice\n",
        "a = 0 #lower integration limit\n",
        "b = 1 # upper intergration limit\n",
        "h = (b-a)/N #width of slice\n",
        "error = 10**(-6) #approximate accuracy\n",
        "err = [inf] #list for creating the errors until desired approximate error is reached\n",
        "slices = [] #list for storing the number of slices required to reach the desired integral given error"
      ],
      "metadata": {
        "id": "aeEo3-7gRNtn"
      },
      "execution_count": 987,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "#defining all the needed functions\n",
        "def TrapI(N):\n",
        "    '''trapezoidal rule that will be used for romberg integration'''\n",
        "    h = (b-a)/N\n",
        "    s = 0.5*g(a) + 0.5*g(b)\n",
        "    \n",
        "    for k in range(1,N):\n",
        "        s += g(a+k*h)\n",
        "    return(h*s)\n",
        "\n",
        "I1 = 0\n",
        "I2 = TrapI(N)\n",
        "\n",
        "while abs(1/3*(I2-I1)) > error: #continues the assigning of new values of I1 and I2 until it reaches the accuracy needed\n",
        "    h_i = (b-a)/N\n",
        "    I1 = I2\n",
        "    I2 = 1/2*I1\n",
        "    for k in range(1,N,2):\n",
        "        I2 += h_i*g(a+k*h_i)\n",
        "    N = 2*N\n",
        "\n",
        "def Romberg():\n",
        "    '''Romberg integration with no argument needed'''\n",
        "    size = 7\n",
        "    R = np.zeros((size,size))\n",
        "\n",
        "    for i in range(0,size):\n",
        "        N = 2**i\n",
        "        slices.append(N)\n",
        "        R[i,0] = TrapI(N) # using the trapezoidal function above to calculate R[i,0]\n",
        "        for k in range(0,i): \n",
        "            m = k + 2\n",
        "            R[i,k+1] = 1.0/(4**(m-1)-1)*(4**(m-1)*R[i,k] - R[i-1,k]) #calculating the next trapezoidal estimate\n",
        "        if (i > 0):\n",
        "          if (np.inf >= error):\n",
        "            err.append(abs(R[i,k+1] - R[i,k]))\n",
        "\n",
        "    return R, err, slices #return the value R, err, and slices\n",
        "\n",
        "R, err, slices = Romberg() # assigns the value of Romberg into R, err, and slices"
      ],
      "metadata": {
        "id": "L6alcgbkR10_"
      },
      "execution_count": 988,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "values = [] #list for the R[i,i] values which is the final integration per step\n",
        "for i in range(7):\n",
        "  values.append(R[i,i])"
      ],
      "metadata": {
        "id": "VDROFXGHSgCB"
      },
      "execution_count": 989,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "df = pd.DataFrame(\n",
        "    {'Number of slices': slices,\n",
        "     'Error': err,\n",
        "     'Value of the integral': values\n",
        "    })\n",
        "df.index += 1\n",
        "df #displays the integral given the number of slices and its error"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 266
        },
        "id": "7yRcsXUCSfk-",
        "outputId": "b38d08e3-8e5a-450b-ed85-1583d15d15ed"
      },
      "execution_count": 990,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "   Number of slices         Error  Value of the integral\n",
              "1                 1           inf               0.147979\n",
              "2                 2  5.908414e-02               0.384316\n",
              "3                 4  1.268781e-02               0.587321\n",
              "4                 8  3.724174e-03               0.348974\n",
              "5                16  3.782418e-04               0.445804\n",
              "6                32  9.739968e-06               0.455777\n",
              "7                64  1.342828e-08               0.455832"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-359e8027-3c5b-4dce-8b67-fc2e6be9d620\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Number of slices</th>\n",
              "      <th>Error</th>\n",
              "      <th>Value of the integral</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>1</td>\n",
              "      <td>inf</td>\n",
              "      <td>0.147979</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>2</td>\n",
              "      <td>5.908414e-02</td>\n",
              "      <td>0.384316</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>4</td>\n",
              "      <td>1.268781e-02</td>\n",
              "      <td>0.587321</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>8</td>\n",
              "      <td>3.724174e-03</td>\n",
              "      <td>0.348974</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>16</td>\n",
              "      <td>3.782418e-04</td>\n",
              "      <td>0.445804</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>6</th>\n",
              "      <td>32</td>\n",
              "      <td>9.739968e-06</td>\n",
              "      <td>0.455777</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>7</th>\n",
              "      <td>64</td>\n",
              "      <td>1.342828e-08</td>\n",
              "      <td>0.455832</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-359e8027-3c5b-4dce-8b67-fc2e6be9d620')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-359e8027-3c5b-4dce-8b67-fc2e6be9d620 button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-359e8027-3c5b-4dce-8b67-fc2e6be9d620');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ]
          },
          "metadata": {},
          "execution_count": 990
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "#plotting the values \n",
        "x = []\n",
        "y = []\n",
        "for i in range(0,len(data)):\n",
        "     x.append(data[i][0])\n",
        "\n",
        "for i in range(0,len(data)):\n",
        "     y.append(data[i][1])\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "adtrap = ax.plot(x,y,'salmon', marker='o', label='Adaptive Trapezoidal Rule (ATR)')\n",
        "romberg = ax.plot(slices, err, 'teal', marker='o',label='Romberg Integration (RI)')\n",
        "ax.legend()\n",
        "ax.set_yscale('log')\n",
        "ax.set_xscale('log')\n",
        "ax.set_xlabel(\"no. of slices N\")\n",
        "ax.set_ylabel(\"error estimate |err|\")\n",
        "ax.axhline([10**(-6)], color = 'tan', linestyle = '--')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 302
        },
        "id": "W19eCo-pTw9o",
        "outputId": "b83b11fa-092f-4045-84c3-d6f0dbcb8aea"
      },
      "execution_count": 991,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "<matplotlib.lines.Line2D at 0x7f7b1b9ab6d0>"
            ]
          },
          "metadata": {},
          "execution_count": 991
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEMCAYAAAA1VZrrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd3iUVfbA8e+ZSYAkQOiiFCkm1IQgAUFESkBAEpqgFEFFYWXtu7p2V1356a4NCxYURBECClKCoCgdQWkiIEhvQaQJAUJJZub+/pgkJmFCJslMZjI5n+eZJ847bzkzhty57733HDHGoJRSSl2OxdcBKKWU8n/aWCillMqXNhZKKaXypY2FUkqpfGljoZRSKl/aWCillMqXNhZKKaXypY2FUkqpfAX5OgB3iEgD4Gkg3BgzIL/9q1WrZurVq+f1uJRSKpCsX7/+uDGmuqvXvN5YiMhEIB44aoxpnm17D+AtwAp8bIx5Ja9zGGP2AHeLyAx3rlmvXj3WrVtXtMCVUqqUEZH9eb1WHD2LScC7wGfZArIC44BuQDKwVkTm4mw4Xs51/AhjzNFiiFMppVQevN5YGGOWi0i9XJvbALsyegyIyDSgjzHmZZy9kAITkVHAKIC6desWOl6llFKX8tUAdy3gYLbnyRnbXBKRqiLyAdBSRJ50tY8xZrwxJtYYE1u9ustbbkoppQqpRAxwG2NOAPfmt5+IJAAJ11xzjfeDUiVWeno6ycnJXLhwwdehKOUT5cqVo3bt2gQHB7t9jK8ai0NAnWzPa2dsUx5g37wex6IFkHISwitjieuJNaqVr8PyG8nJyVSoUIF69eohIr4OR6liZYzhxIkTJCcnU79+fbeP81VjsRaIEJH6OBuJQcCQop7UGJMEJMXGxo4s6rm8ZfLXc3hm/UYOGkMdEV5qFcOwXn08dn775vU4kr6E9HTnhpSTzuegDUaGCxcuaEOhSi0RoWrVqhw7dqxAx3l9zEJEEoHVQCMRSRaRu40xNuB+4FtgG/CFMeZXD1wrQUTGp6SkFPVUXjH56zncu+5nDhgwCAcM3LvuZyZ/PcftcxiHA5Oehjl/DnP2NObUn5gTxzBHDuP4/SCOb5P+aigypac7exoqizYUqjQrzO9/ccyGGpzH9vnAfA9fq0g9C0/fvjE2G6SewZw9A2fP8Mz6nzlHzv9J5xCeWbeBwcd/B7sNbDZMxk/s9oyf2f7bOAoXTMpJHIcOIFfWRiy6cN8fzJ49m379+rFt2zYaN27scp9OnTrx2muvERsbW6jzR0ZG0rRpUwCee+45brzxRrp27VromO+77z5++OEH0tLS2Lt3L40aNQLgmWeeYcCAfNfLetXcuXPZunUrTzzxxCWvlS9fnrNnz172+Lz2sVqtREVFYbPZqF+/PpMnT6ZSpUp5nmfp0qW89tprzJs3z+3YDx8+zMiRI3Mc8/DDD/Pll19y8OBBfv31V4YNGwbAgQMHCA8PJzw8nGrVqvHxxx/TpEkTGjVqRFpaGrGxsUyYMIHg4GA2b97M66+/zqRJk9yOJS8lYoDbXUUZ4Hb39o2x2yD1LJzNaASyGoOzmNTTzp9nTzv3uXA+xzUOUsHltQ8izgYhKAjKlkOCgsAa5HxuDUKs1qz/zvEzKAjJtc0+d7rz2q7e48dvQUgo0iACS4NGSMNIJLxygT+r0sZbY0CJiYnccMMNJCYm8sILL3gg0pxmz55NfHx8VmPx4osvFvmc48aNA2Dfvn3Ex8ezcePGHK/bbDaCgnzzZ6V379707t3b4+cNCQnJep933HEH48aN4+mnn/boNd544w1GjvzrO67D4WDWrFnUqVOHZcuW0blz56wY7rzzTuLj47Ma53379tGwYUM2btyI3W6nW7dufPHFFwwdOpSoqCiSk5M5cOBAkZcUBFRjUZSehWPRAte3b5K+xGxYk9UwcP6c6xOULQdh5ZHyFZErrkTCKkD5Clk/KV+BOhM+4YCLkuc1RQga8UBBQ3bJdO+ds9EDCA5GuiVgKReCY88OzO7t2H/9xflatRp/NRz1GiJlynokjkDhrTGgs2fPsnLlSpYsWUJCQkJWY3H+/HnuuusufvnlFxo3bsz583994Rg9ejRr167l/PnzDBgwIOuYevXqceutt7JgwQJCQkKYOnUqR48eZe7cuSxbtoyXXnqJmTNn8p///If4+HjKly/PhAkT+PJL5/vI/k144cKF/Pvf/+bixYs0bNiQTz75hPLly1/2vSxdupRnn32WypUr89tvv7Fjxw769u3LwYMHuXDhAg899BCjRo0CnN/eR44cycKFC6lZsybTpk2jevXq7N69m/vuu49jx44RGhrKRx99ROPGjYmJicm6zvbt2/nmm2+IiopixIgR7Nmzh9DQUMaPH090dDSTJk1i3bp1vPvuu+zdu5chQ4Zw9uxZ+vT5a0ww8/nJkydJT0/npZdeyvF6ftq1a8emTZuAnL2+48ePExsby759+3Lsn5qaygMPPMCWLVtIT0/n+eefd3m9mTNn8tJLL+X4TJs1a8Ztt91GYmIinTt3dis+q9VKmzZtOHTor/lCCQkJTJs2jX/9619uv09XAqqxKNLU2ZSTrrenp4PdhlSvgdRr6PzDH1YByWgAJPO5G1PQXmoVw73rct+KMpwSC3O3b6d3Rpe+KDL/gOX1TdgSdS3GGDh2BMfu7Zg9O3Bs+BHWrACLFalTz9lwNGyE1KwV8Les7N/Mxvzxe56vm+T9zl5fdunpOOZ8gVn/k8tjpOZVWHv0vex158yZQ48ePYiMjKRq1aqsX7+eVq1a8f777xMaGsq2bdvYtGkT1157bdYxY8aMoUqVKtjtduLi4ti0aRPR0dEAhIeHs3nzZj777DMefvhh5s2bR+/evXN8A83UtWtXRo0aRWpqKmFhYUyfPp1BgwZx/PhxXnrpJb7//nvCwsL473//yxtvvMFzzz132fcCsGHDBrZs2ZI1u2bixIlUqVKF8+fP07p1a2655RaqVq1KamoqsbGxvPnmm7z44ou88MILvPvuu4waNYoPPviAiIgIfvrpJ/7+97+zePHirG/TSUlJ/O9//+P666/nH//4By1btmT27NksXryY4cOHX9LDeeihhxg9ejTDhw/P6g2Bc8rorFmzqFixIsePH6dt27b07t3brXv4drudRYsWcffdd+e7b6YxY8bQpUsXJk6cyKlTp2jTpg1du3YlLCwsa5+9e/dSuXJlypb964taYmIigwcPpk+fPjz11FOkp6e7Nc31woUL/PTTT7z11ltZ22JjY3nllVe0sciuSGMW4ZVdNxjhlT32rT9z1lP22VAPN2/GlOMn6TNtGk936MALnTphLeIfaGtUq8t+6xURqFETa42a0K4jxpaOObgPs3s7jt07MIsXwOIF+d6yKhVTdHM3FPltd1NiYiIPPfQQAIMGDSIxMZFWrVqxfPlyHnzwQQCio6OzGgOAL774gvHjx2Oz2Th8+DBbt27Nen3w4MFZPx955JHLXjsoKIgePXqQlJTEgAED+Prrr/nf//7HsmXL2Lp1K+3btwcgLS2Ndu3aufV+2rRpk2Ma5ttvv82sWbMAOHjwIDt37qRq1apYLBZuu+02AG6//Xb69+/P2bNnWbVqFQMHDsw6/uLFi1n/vXPnTh577DGWLFlCcHAwK1euZObMmQB06dKFEydOcPr06Rzx/PDDD1n7DBs2jMcffxxwTht96qmnWL58ORaLhUOHDnHkyBFq1qyZ53s7f/48MTExHDp0iCZNmtCtWze3PhOAhQsXMnfuXF577TXA+cf8wIEDNGnSJGufw4cPk30hcVpaGvPnz+eNN96gQoUKXHfddXz77bfEx+ed3GL37t3ExMSwd+9eevXqleP3pkaNGvz+e95fiNwVUI1FUVjierq8fWOJ6+nR6wzr1eeSqbKjbTbunz+fMStWsObQIabecgvVQkM9et3LkaBgpH4E1I/A2hVM6hnMnp2XvWXlSD2NWTC7xE/Rza8HkD72pby/RNz590Jd888//2Tx4sVs3rwZEcFutyMivPrqq3kes3fvXl577TXWrl1L5cqVufPOO3MsKsz+zdidb8mDBg3i3XffpUqVKsTGxlKhQgWMMXTr1o3ExMQCv6fs35SXLl3K999/z+rVqwkNDaVTp055LoAUERwOB5UqVbqkdwDO20a33norH330EVdeeWWBYnL1OUyZMoVjx46xfv16goODqVevXr6LMzPHLM6dO0f37t0ZN24cDz74IEFBQTgczgkneZ3DGMPMmTOzJgLkdf7sx3/77becOnWKqKgoAM6dO0dISMhlG4vMMYvjx4/Tvn175s6dmzV+c+HCBUJCQi77Ht0RUPcYijJ11hrVCkvCQGcPA5zflBMGFssfvnJBQXzcuzcfJSSwfP9+Wo0fzzoPfBMoLAmrgCXqWoL6DCLokecIGv0Ylpt6I5Wq4NjwI/bECZi5X5aKKbqWuJ6Qu/tfxC8RM2bMYNiwYezfv599+/Zx8OBB6tevz4oVK7jxxhuZOnUqAFu2bMm6P3769GnCwsIIDw/nyJEjLFiQ83OePn161s/M3kCFChU4c+aMyxg6duzIhg0b+Oijjxg0aBAAbdu25YcffmDXrl2A8377jh07Cvz+UlJSqFy5MqGhofz222/8+OOPWa85HA5mzHAmj546dSo33HADFStWpH79+lljKMYYfvnF+QVlxIgR3HXXXXTo0CHrHB06dGDKlCmAs2GqVq0aFStWzBFD+/btmTZtGkDWvpmx1ahRg+DgYJYsWcL+/XkmWb1EaGgob7/9Nq+//jo2m4169eqxfv16gKz3lFv37t155513nLd+gZ9//vmSfSIjI3OMdSQmJvLxxx+zb98+9u3bx969e/nuu+84dy6P8dJsqlWrxiuvvMLLL/+Vj3XHjh00b978Mke5J6AaC2NMkjFmVHh4eKGOt0a1IvjhZwj+9+sEP/xMsX9Dvufaa1k5YgQA7SdO5OMNG4r1+q6ICFKjJtZ2HQkaOpKgx/+Dddjf8j4gr7GfEsobXyISExPp169fjm233HILiYmJjB49mrNnz9KkSROee+45WrVyXqdFixa0bNmSxo0bM2TIkKxbRZlOnjxJdHQ0b731Fm+++Sbg7D28+uqrtGzZkt27d+d8X1Yr8fHxLFiwIOsba/Xq1Zk0aRKDBw8mOjqadu3a8dtvvxX4/fXo0QObzUaTJk144oknaNu2bdZrYWFhrFmzhubNm7N48eKs8ZApU6YwYcIEWrRoQbNmzZgzZw779+9nxowZTJw4kZiYGGJiYli3bh3PP/8869evJzo6mieeeIJPP/30khjeeustxo0bR1RUVI7B3qFDh7Ju3TqioqL47LPP8pyynJeWLVsSHR1NYmIijz76KO+//z4tW7bk+PHjLvd/9tlnSU9PJzo6mmbNmvHss89esk9YWBgNGzZk165dnDt3jm+++YZevXrleP2GG24gKSnJrRj79u3LuXPnWLFiBQBLlizJcb7CkswWL5DExsaaklzP4vi5cwz96isW7t7NiJgY3r35ZkIKkMOlOOR5e6ZiJYIfufQfhD/Ztm1bjnvGJV1m/ZZq1ar5OpR8ubPeoTSaNWsW69evzzEjyhMuXrxIx44dWbly5SVTml39OxCR9cYYl4t6AqpnESiqhYYyf8gQnunQgYkbN3LDJ5+w79QpX4eVg8vbMwB2O+b3g5duV0rlqV+/fnijuueBAwd45ZVXPLL2JaB6Ftmmzo7cuXOnr8PxiKTt2xk2axZWi4Wp/fvT3Y8y6uaeDSVR12J+WQepZ7B07onl+k5+OfU20HoWShVGQXsWAdVYZCrpt6Fy2/Xnn/SfPp0tR4/yQqdOPH3jjVj8NLeROX8O+7wZmK2/IFc3xNpvsN+tEtfGQim9DRWQrqlShR/vuYeh0dE8t3QpvRMTOXn+fP4H+oCEhGIdMAxrn0GYw8nY3n8Nx5ZLZ4AopUoWbSxKiNDgYD7r25d3e/Zk4e7dxH70Eb/88Yevw3JJRLDEtCbob/9Aql+Bfebn2GZNxVzUYkNKlVTaWJQgIsJ9bdqw7M47uWCz0W7CBCZnzEf3R1KlGta77sPS8SbM5g3YPngdx4G9vg5LKVUIAdVY+Hs9C09pV6cOG0aNok2tWgyfPZv7vv6aNLvd12G5JBYr1k7dsd51P4hgnzQO++IFGD+Nt7hYrVZiYmJo3rw5CQkJnPLAbLd9+/Z5ZPGVp67zf//3f16PBWDjxo3Mn/9XtYO5c+fyyiuveOTc58+fp2PHjtjtdvbt20dISAgxMTE0bdqU4cOHk56xMHXp0qVZ61XmzZvnVj6tkiagGouiLsorSa4oX57vhw/n0XbteG/dOjpOmkRyrvw4/sRSp57ztlSLWBwrvsf+ybuYEwWr1OUrUzZvpt7YsVheeIF6Y8cyZfPmIp8zM4XEli1bqFKlSo5kd75isxUt31VunmwsLhdb7said+/eLmtaFMbEiRPp378/VqsV+CutxubNm0lOTuaLL7645JhevXqRlJTk1orrkiSgGovSJshi4dWbbuLLgQPZcvQo1374IUv2+u9tHilbjqA+g7AOGI45cQzbh2/g2PAT/jwjb8rmzYxKSmJ/SgoG2J+SwqikJI80GJnatWuXtcp448aNtG3blujoaPr168fJk86Fj506deKRRx4hNjaWJk2asHbtWvr3709ERATPPPNM1rlsNhtDhw6lSZMmDBgwIOsP1vr16+nYsSOtWrWie/fuHD58OOu8Dz/8MLGxsbz11lusXbuW6OhoYmJieOyxx/LtQUyaNIn+/fvTo0cPIiIisjKbPvHEE1kJ+IYOHQrA559/Tps2bYiJieFvf/sb9oze5YQJE4iMjKRNmzaMHDmS+++/H3DWbbj33nu57rrr+Ne//sWaNWto164dLVu25Prrr2f79u2kpaXx3HPPMX36dGJiYpg+fTqTJk3KOse+ffvo0qUL0dHRxMXFceDAgaxzP/jgg1x//fU0aNAgz3QdU6ZMcZlS3FUq8EwiQqdOnQpU/Kgk0ESCAWBA06Y0r1GD/tOn03XyZG5t2pRVyckcTEmhbng4Y+LiGJqRlMwfWJq1QOpcjX12IvakL5Cd27AmDERCw/I/2MMe/uYbNl5mosCPyclczHXL7Fx6OnfPmcNHGXmBcoupWZOxPXq4df3caa+HDx/OO++8Q8eOHXnuued44YUXGDt2LABlypRh3bp1vPXWW/Tp04f169dTpUoVGjZsmJVpdvv27UyYMIH27dszYsQI3nvvPR566CEeeOAB5syZQ/Xq1Zk+fTpPP/00EydOBJxZTjOnmjdv3pyPPvqIdu3auf3tfOPGjfz888+ULVuWRo0a8cADD/DKK6/w7rvvZiUH3LZtG9OnT+eHH34gODiYv//970yZMoWuXbvyn//8hw0bNlChQgW6dOlCixYtss6dnJzMqlWrsFqtnD59mhUrVhAUFMT333/PU089xcyZM3nxxRez6lgAOarCPfDAA9xxxx3ccccdTJw4kQcffJDZs2cDzmyvK1eu5LfffqN3796XpHJPS0tjz549LhfLuUoFnl1sbCwrVqzg1ltvdeszLAm0ZxEgGlerxk/33EPslVcy7ddfOeDFb8KeIBUrYR32Nyzd4jE7tjqn2O7e7uuwLpG7ochvu7syv3XXrFmTI0eO0K1bN1JSUjh16hQdO3YEnFXZli9fnnVMZhbRqKgomjVrxpVXXknZsmVp0KABBw86V83XqVMnK2/U7bffzsqVK9m+fTtbtmyhW7duxMTE8NJLL5GcnJx13syU4adOneLMmTNZiQiHDBni1nuJi4sjPDyccuXK0bRpU5fJ+RYtWsT69etp3bo1MTExLFq0iD179rBmzRo6duxIlSpVCA4OzpGmHGDgwIFZt4BSUlIYOHAgzZs355FHHuHXX3/NN7bVq1dnvY9hw4axcuXKrNf69u2LxWKhadOmHDly5JJjjx8/fkn51MxU4FdccQVXXnlljlTg2XkqLbg/0Z5FAKlQtix/pKZesv1cejpPL1rkV70LABEL1us7Y6kfge2rqdg/H49peyOWuJuRoOLJhZVfD6De2LHsdzFh4urwcJbeeWehr+sq7fUdd9xx2WMyi+NYLJYchXIsFkvWPf3cablFBGMMzZo1Y/Xq1S7Pmz29eGFkj8VqtbocXzDGcMcdd+TIhgpkfcvPS/bYnn32WTp37sysWbPYt28fnTp18ljcrm6F5k4dDpdPBZ6dp9KC+5MS0bMQkb4i8pGITBeRm3wdjz87mMdMsAN+PENMrqxN0KiHsbRuj+PH5dg+egtz9LCvwwJgTFwcoblyYIUGBzMmLs4j58+e9josLIzKlStnZQudPHlyVi/DXQcOHMhqFDJTgDdq1Ihjx45lbU9PT3f5rbxSpUpUqFCBn35yVgDMTPFdWMHBwVmzheLi4pgxYwZHjx4FnDU99u/fT+vWrVm2bBknT57EZrNlFSxyJSUlhVq1agE5bzVdLhX79ddfnyNVefZU5/mpXLkydrvdZa0KV6nAs/NUWnB/4vXGQkQmishREdmSa3sPEdkuIrtE5LI3R40xs40xI4F7gdu8GW9JVzePmWAiwuRffvHbwWQJLoP15v5Yh9wDqWewjR+L/cfl2DatI33sS6S/8E/Sx76EfbPrcQJvGRoVxfiEBK4OD0dw9ijGJyR4tJeWPe31p59+ymOPPUZ0dDQbN24s8BTMRo0aMW7cOJo0acLJkycZPXo0ZcqUYcaMGTz++OO0aNGCmJgYVq1a5fL4CRMmMHLkSGJiYkhNTaUoMwtHjRpFdHQ0Q4cOpWnTprz00kvcdNNNREdH061bNw4fPkytWrV46qmnaNOmDe3bt6devXp5XvNf//oXTz75JC1btszRe+ncuTNbt27NGuDO7p133uGTTz4hOjqayZMn5znGkJebbropx62r7HKnAs/OU2nB/YnXc0OJyI3AWeAzY0zzjG1WYAfQDUgG1gKDASuQu6keYYw5mnHc68AUY8xlCz0EWm6ogsicvXMuW2GickFB1Cpfnt2nTtGlfn3e79WLyKpVfRjl5ZnUM9jnTMfs3AYikP13NDi4yPUkNDdU3s6ePUv58uUBeOWVVzh8+HCB/8AW9po2m41+/foxYsSIS+p9+MqGDRt48803mTx5stvHHDlyhCFDhrBo0SIvRlZ0fpcbyhizHPgz1+Y2wC5jzB5jTBowDehjjNlsjInP9TgqTv8FFuTVUIjIKBFZJyLrjh0rGfP3vcHVN+GPe/dmx4MP8n6vXqz//Xei3n+fF5Yu5aKH59V7ioRVwDr4bigXkrOhgICsxudPvv7666zFgitWrMgxLddbnn/++axr1q9fn759L1/qtjhde+21dO7cOWuarzsOHDjA66+/7sWofKNYss6KSD1gXraexQCghzHmnoznw4DrjDH353H8g8AdOHsgG40xH+SxX8ClKPe0P86e5ZFvv2Xali00qlqV93v1onP9+r4Oy6X0F/6Z52uWNjcgV9ZGrqoN1WogFqvb59WehVIF71mUiNlQxpi3gbfd2C8JSIqNjR3p/ahKpprly5N4yy3c2aIFf58/ny6ffcbwFi14rVs3qhdxVozHhVd2XY3PasXx8xpYk3EvOSgYqXlVVuMhV9aB6gVrQJRSl+erxuIQUCfb89oZ24okW8+iqKcKeN2vuYYto0fz0vLlvLpqFfN27ODVbt24KybmkumXvmKJ64kj6UvINv6SOWZhadYSThzDHE7G/J6MOZyM45d1sPYH537ZG5DMHkj1K7IaEGNM1vs051IxZ1LAbgerFakQ7pMFgkoVl8LcUfLVbaggnAPccTgbibXAEGNM/qts3FCaB7gLY+uxY/xt3jxWHjhAh7p1+SA+nqbVq/s6LODSanyWuJ55Dm4b43A2IBmNh/NxCNIuOncICkKuuIoDzVtToU5dqlavjthsmJSTOcdGRJDwytpgqIBkjOHEiROcOXOG+rluQfu0Up6IJAKdgGrAEeDfxpgJInIzMBbnDKiJxpgxHriWjlkUksMYPvn5Zx777jvOpqXxr/btebpDB0Jc1dkuQZwNyPEcPZD0E8f4vXkrLlasDOTRi7JYkIqBn5BSlU7lypWjdu3aBOf6961lVZXbjqam8ujChUzetImGlSvzXq9e3NSwoa/D8ihjHPDnCczvB7F/NSXP/YL/HXgzWpS6nFJTVrW01LPwphphYXzWrx+Lhg/HarHQ/fPPGTJzJn+cPevr0DxGxIJUrY4l6lrnILorwcGYM/6b8l2p4qY9C5WnCzYb/125kv9buZKQoCBe6dqVUa1aYfGTAXBPsG9ef+kgusXiHMMoUxZL5x5YWrdHLAH1vUopl0rNbSgds/CO7cePM/rrr1mybx/tatfmw/h4oq64wtdheYyrQXTLVXWxz/8Ks2cH1KyFtdctWGpf7etQlfKqUtNYZNKehecZY/h80yb+sXAhJ8+f5x/t2tG4WjVeXLaMA35aN6OojDGYrZuwfzsbzpzB0qqtMyNuSKivQ1PKK7SxUB5z4tw5Hv/+eyb8/DMCZP/tCQ0O9niSPX9gLl7AsfRbHD+thJBQrN3ikRaxfrMeRSlP0QFu5TFVQ0P5uHdvrggLI/fXjMy6GYFGypbD2r0PQaMeQapUxT5nGvZJ7/lNGnWlikNANRbGmCRjzKiipFVW7jnqosgS+HfdjKKSmldhHXE/1oRbMcf+wPbhG9i/S8JkLvpTKoCViNxQyv/UDQ93WUEur3oagULEglx7HdK4Gfbvv8axaimOLRux9uiLNG6ut6ZUwAqonoUqPq4qyAHc17q1D6IpfhJanqDet2EdcT+EhGD/YhL2xAmYkyd8HZpSXhFQjYWOWRSf3HUzalWoQMUyZXhv3TqOBNACvvxY6tQnaNQjWG7qjdm/B9t7/8O+/DuMn9YKUaqwdDaU8pi1hw7RcdIkoq64giV33OGy5xHIzOkU7N/OwWz9BapWx3pzfywNIn0dllJuKzWzoZRvta5Vi6m33MLaQ4cYPmsWjgD8InI5UjGcoIHDsQ4dCcZgn/whthmTnenPlSrhtGehPO7N1av5x8KFPHb99fyvWzdfh+MTxpaOY+ViHJlO1jgAACAASURBVCsXg9WKpXMPTEgIZsm3bqVbV8oXSnylPFWyPNy2Lbv+/JNXV62iYeXK/C3W5e9eQJOgYKydumOJboV9/lc4vp0D2Zcxppx05qQCbTBUiRBQt6F0gNs/iAhv9ezJzRER3Dd/Pt/s2uXrkHxGqlRz3pYKDYPcyxjT0505qZQqAQKqsdBFef4jyGJh2i23EHXFFdz65ZdsOnLE1yH5jIjAOdeLGF3WGFfKDwVUY6H8S4WyZZk3eDAVy5al19Sp/H7mjK9D8p286maULavTbFWJoI2F8qpaFSsyb8gQTl24QPzUqZxNS/N1SD5hiesJuacSiwUuXsT20ZuYw8m+CUwpN2ljobwupmZNpg8YwC9HjjBk5kzsDoevQyp21qhWWBIG/tXDCK+Mpd8grEPugXPnsH38FvZl32Ecdt8GqlQeSsTUWRFpAjwEVAMWGWPev9z+OnXWP723di33zZ/Pg23a8FbPnr4Ox2+Y8+ewL/gKs/ln5Ko6WPsORqoHTnEpVXL4dFGeiEwUkaMisiXX9h4isl1EdonIE5c7hzFmmzHmXuBWoL0341Xe8/fWrflH27a8vWYNb//0k6/D8RsSEkpQ/9uxDhyOOXnCmc129TKMKX09MOW/imOdxSTgXeCzzA0iYgXGAd2AZGCtiMwFrMDLuY4fYYw5KiK9gdHA5GKIWXnJ/7p1Y8+pUzz8zTfUq1SJ3o0a+Tokv2Fp2gKp2wD7vC9xLJyL2b4Fa59BSOWqvg5NqeK5DSUi9YB5xpjmGc/bAc8bY7pnPH8SwBiTu6Fwda6vjTG9LreP3obyb6lpaXT69FO2HjvG8jvvpNVVV/k6JL9ijMH8sg77N7PB4cDavTdybVtNf668zh9zQ9UCDmZ7npyxzSUR6SQib4vIh8D8PPYZJSLrRGTdsWPHPBut8qiwMmVIGjyYaqGhxCcmBnTBpMIQESwxrQka/ShS+2rs82Zgn/ox5rR+Tsp3LttYiMheEdmT7WfuR+b2B70ZpDFmqTHmQWPM34wx4/LYZzzwArChTJky3gxHeUDN8uWZP2QI59LTiZ86ldMXtdpcbhJeGeuwUVh69sPs243t/VdxbFpPSZiUogLPZRsLY0x9Y0yDbD9zPzK3v13A6x4C6mR7XjtjW5HoCu6SpVmNGsy89Va2HT/OwC+/JN2u00ZzE7FgbXMDQff+E6lWA/usqdi//AyTWnpqhij/kO9tKBGxisgSD193LRAhIvVFpAwwCJhb1JNqbqiSp2uDBnzQqxcLd+/m/vnz9VtzHqRqdax33Y8lrhdmx6/OXsZvW/I/UCkPybexMMbYAYeIFOrruogkAquBRiKSLCJ3G2NswP3At8A24AtjzK+FOb8q+e6+9lqevOEGxm/YwGurVvk6HL8lFgvWG7oQNPIRqFAR+/RPsM1OxFw47+vQVCng1mwoEZkDtAS+A7IyohljvDpWUVg6G6rkcRjDkJkzmf7rr3w5cCADmjb1dUh+zdhtOJZ/j2PFIqhQAWufQVqVTxWZJ+pZfJXx8GsikgAkXHPNNb4ORRWQRYRJffty8PRphs2aRe2KFWlbu7avw/JbYg3C2rkHEtkU++xE7JM/xMRej7mqFmbZ91pgSXlcvj2LjAV03xtjOhdPSEWnPYuS61hqKm0nTODMxYv8eM89NKicR7ZWlcWkp+NYPB/Hj8svfTE4GEvCQG0wlFuKtM6iqGMWShVE9bAw5g8Zgs3hoNfUqZw8r/fj8yPBwVi794Gw8pe+qAWWlIe4exvqLLBZRPx6zEJvQwWGRtWqMeu22+g2eTLtJ04kNT2dgykp1A0PZ0xcHEOjonwdon/KazqtFlhSHuDuCu6vgGeB5cD6bA+/oussAkfHevW4+9pr2Xb8OAdSUjDA/pQURiUlMWXzZl+H55/yKrBULkSTEqoic6uxMMZ8CnwB/GiM+TTz4d3QCk7XWQSWBTt3XrLtXHo6Ty9a5INo/J/rAksCF847B8C1h6GKwK3GIuP2zkbgm4znMRlZYv2K9iwCS145ozSXlGuuCixJ30FYEwZikg9g++B1HJs3+DZIVWK5O2bxPNAGWApgjNkoIg28FJNSANQND2e/i4ahrn4ZyJM1qpXLmU9S7xpnqpCvpuDYsRXrzf2RkNDiD1CVWO6OWaQbY3L/q9WboMqrxsTFEZrrtkq5oCDGxMX5KKKSS6pUw3rXfVg698Bs/QXbB6/h2HvpbT6l8uJuY/GriAwBrCISISLvAH6Xl0HHLALL0KgoxickcHV4OIJz4V5ElSo6G6qQxGLFemM3rHc/CMFlsH/2AfZv52Bs6b4OTZUA7qb7CAWeBm4CBGdOp/8YYy54N7zC0UV5gemN1av558KFfD9sGHEN9C5oUZj0NBzfzcOx9geoXpOg/kORmlqEqrS73KK8YqmUV9y0sQhMF2w2It95hyvKl2fNPfdo5TgPcOz6DfucaXDuHJYuPbC064RYfFUTTflaoXNDichwN6+x0RizqcCRKVUA5YKC+E/nztw5Zw4ztm5lYLNmvg6pxLNc0xgZ/Sj2eTNwfP81Zuc2rH0HI5Wq+Do05Wfy+wpR381HJS/GqFSW26OjaV6jBk8tXqzFkjxEQstjHXgH1j6DMIcPOafY/rJOa4uoHALqNlS2dB8jd7pY0KUCw9c7dhCfmMh7N9/M6NatfR1OQDEnT2CfnYg5sBdp2gJrr1uQ0DBfh6WKSZESCZYkuiivdLg5IoIbr76aF5Yt42xamq/DCShSuSrWO/6OJe5mzG9bsL3/Go7d230dlvIDAdVYqNJBRPhv164cSU3lzdWrfR1OwHFW5Isj6J4HoVwI9s/HY18wC5OuU2xLM20sVInUtnZt+jVuzKurVnEsNTX/A1SByZW1CRr1CJbrOuBYsxLb+Dcwh5N9HZbyEXfXWQgwFGhgjHlRROoCNY0xa7wdYGHo1NnS4bfjx2n23ns80KYNY3v08HU4Ac2xe7tzim3qWSydumMqVsIs+UYr8gUYT4xZvAe0AwZnPD8DjPNAbEoVWuNq1bi7ZUveW7uWvSc1o6o3WRo2Imj0Y0iTKByLF2DmTPurTkbKSRxJX2Lf7HdVC5QHudtYXGeMuQ+4AGCMOQmU8VpULohImIisE5H44ryu8m//7tiRIIuFZ5cs8XUoAU9CQrHeMgxCQiH3HQmtyBfw3E4kmFGL2wCISHXcTCQoIhNF5KiIbMm1vYeIbBeRXSLyhBunehxnTQ2lstSqWJGHrruOKZs3s/GPP3wdTsATETh/zvWLWi8joLnbWLwNzAJqiMgYYCXwspvHTgJy3FDOaHjGAT2BpsBgEWkqIlEiMi/Xo4aIdAO2AkfdvKYqRR6/4QYqlyvHk1oUqXjkVZEv1EUNcBUw3K2UNwX4F84G4jDQ1xjj1rd8Y8xy4M9cm9sAu4wxe4wxacA0oI8xZrMxJj7X4yjQCWgLDAFGisglcYvIqIzbVOuOHTvmTmgqQFQqV46nO3Tgm127WLx3r6/DCXguK/IBnDuLfcFsnWIboNytlDfZGPObMWacMeZdY8w2EZlchOvWAg5me56csc0lY8zTxpiHganAR8ZFQWFjzHhjTKwxJrZ69epFCE2VRPe1aUOdihV5/PvvNU2Fl7msyNf7towptiuwffSmTrENQO5WysuRsS3jNlKxz5Mzxky63OvZ0n0UT0DKb2iSweKVZ0W+iCbYZ0/D9vHbmsU2wFz2/6KIPCkiZ4BoETktImcynh8F5hThuoeAOtme187YplShaZJB33NOsX0UadQUx/dfY//sA4wOfAeEyzYWxpiXjTEVgFeNMRWNMRUyHlWNMU8W4bprgQgRqS8iZYBBwNwinC8zXs0NVYpZLRZejotj159/MuHnn30dTqkloWEZWWxvwxxOduaX2rzB12GpInJ3gPtJEaksIm1E5MbMhzvHikgisBpoJCLJInK3McYG3I+z4t424AtjzK+FfRPZrqVlVUu5XhERdKhbl+eXLtUkgz4kIlhi2hB07z+RGjWxfzUF28zPMXlNu1V+z910H/cAD+G8XbQR58yk1caYLt4Nr3A03UfptvrgQa6fOJH/dO7MMze69Z1GeZFx2HGsXIxj2UIoXwFr3yFY6uu4oj/yRLqPh4DWwH5jTGegJXDKQ/F5jPYsFEC7OnXo17gx//vhB00y6AfEYsV6YzesIx6A4DLYP/sA+8IkjM3m69BUAbjbWFwwxlwAEJGyxpjfgEbeC6twdMxCZRrTpQup6emMWbHC16GoDJZadZ1ZbFu1xbF6KbaPx2KOHvZ1WMpN7jYWySJSCZgNfCcic4D93gurcLRnoTI1qV6dETExmmTQz0iZsljjB2AdfDecOYNt/FjsPy7HxdIp5WcKXFZVRDoC4cA3Gauv/Y6OWSiAQ6dPc8077zCgaVMm9+vn63BULubsGexJX2B2bEUaRGDtMxipqHcFfMkjZVUzZkNF40xPngw091B8SnlFrYoVefi665iyaZMmGfRDUr4C1kEjsPQagDm4H9v7r+LY+ouvw1J5cHc21H+AO4E9/JVt1vjbbKhsK7hH7ty509fhKD9w6sIFGrz1FtfVrs2CoUN9HY7KgzlxDPtXUzC/H0RaxGLt2Q8pW87XYZU6l+tZuNtYbAei/PW2U256G0pl99qqVTz23XcsHj6czvXr+zoclQdjt+NY/h2OFd87801Ft8L8sk6r8RUjT9yG2gJU8lxIShWf+zXJYIkgVivWzj2w3nUfXLyIWf6dVuPzI+42Fi8DP4vItyIyN/PhzcCU8pRyQUG82Lkza3//nZnbtvk6HJUPS536rlOgazU+n3I36+ynwH+BzbhZIc8XNOusysuw6GheW7WKpxYtok+jRgRbrb4OSV3O6TzW/GpSQp9xt2dxzhjztjFmiTFmWebDq5EVgi7KU3mxWiy80rUrOzXJYMmQVzW+oCDM2TPFG4sC3G8sVojIyyLSTkSuzXx4NTKlPKxXRAQ31K3LC8uWkapJBv2ay2p8FivY7c4sttuLnHdUFZC7jUVLnMkD/w94PePxmreCUsobRIT/du3KH2fPMvbHH30djroMV9X4LH1vI+jeR6FCRezTJmKfNwOTdtG3gZYiBV7BXRLo1Fl1Of2mT2fRnj3seeghqoWG+jocVUDGZsOxZAGOVcugSlWs/YdiqVXX12EFhEJPnRWR2zN+/sPVwxvBFoXmhlLu+L/MJIPLl/s6FFUIEhSEtVsC1uH3gi0d+4R3sC//DuPQ6ojelN9tqLCMnxVcPMp7Ma5C0QFu5Y6sJIPr1rHvlN9l2ldustS/hqB7H0WaReNY8g32Se9hTp7wdVgBy90V3O2NMT/kt81f6G0olZ/MJIMDmzblM00yWOI5Nm/A/vVMMMaZKqRFLCLi67BKHE+s4H7HzW1KlQi1Klbkoeuu4/NNm/hFkwyWeJaoa50lXK+shX3ONOwzPsOc08JXnpTfmEU7EfknUD3XeMXzgK5qUiXa4+3bU6lcOZ5ctMjXoSgPkEpVsA4fjaVrL8xvv2L74DUcu7f7OqyAkV/PogzOsYkgco5XnAYGeDc0pbyrckgIT3XowIJdu1iyd6+vw1EeIBYL1vZdCLrnQSgbgv3z8di/mYOxpfs6tBLP3TGLq40x+zP+2wKUN8ac9nZw2a7fCfgP8CswzRiz9HL765iFctcFm41ar7/O2fR00u126oaHMyYujqFRUb4OTRWRSU/H8V0SjrU/QPWaBPUfitS8ytdh+TVPjFm8LCIVRSQMZwbarSLymJsXnygiR0VkS67tPURku4jsEpEn8jmNAc4C5XAWXlLKI2Zu28bZ9HTS7HYMsD8lhVFJSUzZvNnXoakikuBgrDf3xzrkHjifiu3jsdhXLdUSroXkbmPRNKMn0RdYANQHhrl57CSgR/YNImIFxgE9gabAYBFpKiJRIjIv16MGsMIY0xN4HHjBzesqla+nFy0izZ5zfv659HSe1nGMgGGJaOKcYhvRBMd3Sdg/+xCjCQkLzN3GIlhEgnE2FnONMek4v+3nyxizHPgz1+Y2wC5jzJ6MgkrTgD7GmM3GmPhcj6Pmr68CJ4Gyrq4jIqNEZJ2IrDt27Jibb0uVdgfyWMCZ13ZVMklYeay33ok14VbMoQPYPngdxxZNKFkQ7qYo/xDYB/wCLBeRq3EOchdWLeBgtufJwHV57Swi/YHuOAswvetqH2PMeBE5DCSUKVNGy2kpt9QND2e/i4ahri7sDDgiglx7HVKvIfavpmCf+Tn2NSudac9Pp2g1vny41bPISE9eyxhzs3HaD3T2cmzZr/+VMeZvxpjbLje4rSu4VUGNiYsjNFd209DgYMbExfkoIuVtUqUa1hH3Q+PmcHCfs6EArcaXD7caCxG5QkQmiMiCjOdNgTuKcN1DQJ1sz2tnbCsSzQ2lCmpoVBTjExK4OtsXjGc6dNDZUAFOLFY47OJPjlbjy5O7YxaTgG+BzHlnO4CHi3DdtUCEiNQXkTLAIKDIZVq1Z6EKY2hUFPsefpg//vlPAOwBmIlZuZDXILcOfrvkbmNRzRjzBRklVY0xNsCtFI8ikgisBhqJSLKI3J1x/P04G6BtwBfGmCJXM9GehSqKK8qXp02tWszbscPXoajikFc1PsC+9gcCsXxDUbjbWKSKSFUyZkCJSFvArb/IxpjBxpgrjTHBxpjaxpgJGdvnG2MijTENjTFjChX9pdfSnoUqkoTISNYcOsSRs2d9HYryMpfV+IKCoUZNHPO/wp44QUu4ZuNuY/EPnLeJGorID8BnwANei6qQtGehiio+MhIDzN+509ehKC9zWY2v90CC7n0US89+mL07M0q4brn8iUoJtyvliUgQ0AgQYHvGWgu/pOk+VGEZY6g7diytr7qKr267zdfhKB8yx/7A9tUU+ON3LNe2xdK9N1LG5TKvgOGJdB8YY2zGmF+NMVv8taHQnoUqKhEhPiKChbt3c9Fm83U4yoekek2C7n4Iy/WdcWz4CduHb+A4dMDXYfmM241FSaBjFsoTEho1IjU9naX79vk6FOVjzhKu8VjvuBdstlJdwjXfxkKc6uS3n1KBokv9+oQGB5Oks6JUBku9awga/SjSvEWpLeGab2NhnIMa84shliLT21DKE8oFBdG1QQPm7dih0ydVFikXQlD/27H2H4o5+oczv9TGNaXmd8Td21AbRKS1VyPxAL0NpTwlITKS/SkpbDl61NehKD/zVwnX2tjnTC81JVzdbSyuA1aLyG4R2SQim0VkkzcDU8qXekVEAOgCPeWSs4TrvaWqhKu7jUV3oCHQBUgA4jN++hW9DaU85coKFYi96iodt1B5Km0lXN3NOrsfZ3rwhIxHpcwyq/5Eb0MpT4qPiODH5GSOpQb+LQZVeHJlbYJGPYKldXscPy3HNn4s5sjvvg7L49zNOvsQMAWokfH4XET8bgW3Up6U0KiRruZWbrmkhOtHY7GvXhZQJVzdWsGdMT7RzhiTmvE8DFhtjIn2cnyFoiu4lScYY6j95ptcX6cOXw4c6OtwVAlhUs9in/cl5rctSP1rsPYdjGP/bmfq85STfl1k6XIruN2tlCfkzDJrz9imVMDKXM2duGULaXY7ZaxWX4ekSoDMEq7m5zXYv5mN7Z1XwDggs9Z7RpElwC8bjLy4O8D9CfCTiDwvIs8DPwITvBaVUn4ioVEjzqSlsUxXc6sCEBEs115H0L3/zNlQZCqBRZbcWcFtwdk43AX8mfG4yxgz1suxFZjOhlKe1qV+fcoFBekUWlUoUqXapQ1FphJWZMmdFdwOYJwxZkNGLe63jTE/F0NsBaazoZSnhQYH07VBA5J0NbcqrLyKLIVXKt44isjd21CLROQWEdFxClXqxEdEsPfUKbYdP+7rUFQJ5LLIEgCCOX6k2OMpLHcbi78BXwIXReS0iJwRkdNejEspvxEfGQlA0vbAXqGrvMNVkSW5rgOkXcT24ZvY164qEb3WfGdDZYxZ9DDG/FAM8Sjld2pVrEjLmjWZt3Mnj99wg6/DUSWQNarVJTOfTPsu2OdMwzF/JmbnVqy9b0PKV/BRhPlzd8zi3WKIJU8iYhGRMSLyjojc4ctYVOmUEBnJqoMHOXHunK9DUQFCKlTEOnQklh59s5Vw/dXXYeXJ62MWIjJRRI6KyJZc23uIyHYR2SUiT+Rzmj5AbSAdSC5oDEoVVXxkJA5jWLBrl69DUQFERLBe14GgkY9AxYrYp03EPm8GJu2ir0O7REHHLNIKMWYxCeiRfYOIWIFxQE+gKTBYRJqKSJSIzMv1qIGz9vcqY8w/gNFuXlcpj2l11VXULF9eEwsqr5Aa2Uq4rv/RL0u4urWC2xhT6BtpxpjlIlIv1+Y2wC5jzB4AEZkG9DHGvIwzo20OIpIMpGU8zbee4cXUk+xZPT3HtvArG1G1XgwOezr71nx1yTGVazejcp3m2NLOcWB90iWvV7m6BZWuakza+dMkb7x0MU21BrFUvKIhF8/+yaHN313yeo2ItpSvdjXnU45yeOuSS16/otENhFWpReqfhziyfeUlr1/ZtDMh4TU4e3w/R3f+eMnrtaK6UbZ8FU4f2c3xPZemOqkd05MyIRU59ftv/Ln/l0ter9sqgaAyoZw8uIWTyZd2heu16Y/FGsyJfRtJOXzpQG+DdrcBcGz3Ws4c3ZPjNbEGUb/NLQAc2bma1OM5/xFYy4RwdaveAPzx2wrOncyZhC24XAXqtLwZgN9/XcKF0zlrTJQJq0zt6JsASN60kLTUnPPXy1WswVXNOgNw8Of5pF84k+P10MpXUbNxBwD2r5+LPe18jtfDqtXlioh2xEdE0ObcdnatmoYlWye7Qo0GVG/oLPeS+/cO9HdPf/fc+9078Mt87OXPY25sjjl+FFZMoXzVOlwRNxixWNm7ZibGnrMuvLd/97JzN5GgiMjtIvJsxvM6ItLGnWPzUAs4mO15csa2vHwFdBeRd4DlecQ4SkTWici69PQ0V7soVSTxkZHYHA5SLvrfLQIVOKRcCHJVHSQsDLN3p9+UcHU3keD7gAPoYoxpIiKVgYXGGLeq52X0LOYZY5pnPB+Ac4bVPRnPhwHXGWPuL9S7+Os6CUDCNddcM3KnZgpVHpaalkbV//2P0bGxvNmjR/4HKFVEjs0bsH89E4zB2rMf0iIWby53u1wiQbcr5Rlj7gMuABhjTgJlihDTIaBOtue1M7Yp5bfCypShS/36uppbFZucJVyn+bSEq7uNRXrGoLQBEJHqOHsahbUWiBCR+iJSBhgEzC3C+QBN96G8LyEykt0nT7L9hO9vC6jSwV9KuLrbWLwNzAJqiMgYYCXwf+4cKCKJwGqgkYgki8jdxhgbcD/wLbAN+MIYU+QJxppIUHlbr4zV3JpYUBUnfyjh6taYBYCINAbicNaxWGSM2ebNwIpCix8pb4r54APCy5Vj2Z13+joUVQqZ9HQc3yXhWPsDVK9J0C1DkSuu8si5PVH8CGPMb8BvHonIS7INcPs6FBXAEiIjeXnlSv48f54qISG+DkeVMpklXCWiCfa507F9NBZLXC9MWBhm8Tdeq8bn7m2oEkHHLFRxiI+MxG4M3+hqbuVDlogmBN37KBLRBMfCuZjZ0/6qkZFRjc++eb3nruexMylVSrSuVYsaYWG6mlv5XGYJV0JCIfeQgoer8QVUY6ED3Ko4WEToFRHBN7t2kZ5XFTSliomIwPk8Elx6sBpfQDUWehtKFZeEyEhOXbjADwcP5r+zUt6WZzW+PLYXQkA1FtqzUMWlW8OGlLFadQqt8gsuq/EFBzu3e+oaHjuTH9CehSou5cuUoXO9ejpuofyCq2p8loSBHp0N5fbUWaVUTvGRkTywYAE7TpwgsmpVX4ejSjlX1fg8KaB6FkoVpwRdza1KkYBqLHTMQhWnqytVIqpGDW0sVKkQUI2Fjlmo4hYfGcmKAwc4deGCr0NRyqsCqrFQqrglZBRE+lZXc6sAp42FUkXQplYtqoWG6qwoFfC0sVCqCKwWC70iIpi/cyc2R1FKvCjl3wKqsdABbuUL8ZGRnLxwgdW6mlsFsIBqLHSAW/nCTQ0bEmyx6K0oFdACqrFQyhcqli1Lp3r1dAqtCmjaWCjlAfGRkWw7fpzdf/7p61CU8gptLJTygHhdza0CnDYWSnlAg8qVaVq9uo5bqIBVIhoLEekgIh+IyMcissrX8SjlSkJkJMv27+f0xYu+DkUpj/N6YyEiE0XkqIhsybW9h4hsF5FdIvLE5c5hjFlhjLkXmAd86s14lSqseF3NrQJYcfQsJgE9sm8QESswDugJNAUGi0hTEYkSkXm5HjWyHToEmFoMMStVYO1q16ZKSIjeilIByev1LIwxy0WkXq7NbYBdxpg9ACIyDehjjHkZiHd1HhGpC6QYY87k8fooYBRA3bp1PRO8UgVgtVi4OWM1t93hwGopEXd5lXKLr36bawHZl7smZ2y7nLuBT/J60Rgz3hgTa4yJrV69ugdCVKrgEiIjOXH+PD8mJ/s6FKU8qsR89THG/NsYc9nBbU33oXyte8OGBFksOoVWBRxfNRaHgDrZntfO2KZUiRZerhw3Xn21jluogOOrxmItECEi9UWkDDAImFvUk2puKOUPEiIj+fXYMfaePOnrUJTymOKYOpsIrAYaiUiyiNxtjLEB9wPfAtuAL4wxv3rgWnobSvmcruZWgUiMMb6OweNiY2PNunXrfB2GKsWajBtH3fBwvr39dl+HopTbRGS9MSbW1WslZoDbHdqzUP4iPiKCpfv2cUZXc6sAEVCNhY5ZKH+R0KgRaXY7C3fv9nUoSnlEQDUW2rNQ/uL6OnWoXK4c83bu9HUoSnlEQDUW2rNQ/iLIYqFnRARf79iBXWtzqwAQUI2FUv4kITKSY+fOsfb3330dilJFFlCNhd6GUv6ke8OGWEVI2r7d16EoVWQB1VjobSjlTyqHhNDh6qt13EIFhIBqLJTyN/EREWw6coT9p075hNSZEgAACDxJREFUOhSlikQbC6W8KKFRIwC+1t6FKuECqrHQMQvlbyKrViWiShVNLKhKvIBqLHTMQvmjhMhIFu/dy9m0NF+HolShBVRjoZQ/ylzN/f2ePb4ORalC08ZCKS9rX6cO4WXL6hRaVaJpY6GUlwVbrc7V3Dt34gjALM+qdAioxkIHuJW/io+I4EhqKut0NbcqoQKqsdABbuWvekZEYBHRgkiqxAqoxkIpf1UlJISIypV5ZeVKLC+8QL2xY5myebOvw1LKbUG+DkCp0mDK5s3sOXWK9IwMtPtTUhiVlATA0KgoX4amlFu0Z6FUMXh60aKshiLTufR0nl60yEcRKVUw2lgoVQwO5DHpIq/tSvmbEtFYiEhdEZktIhNF5Alfx6NUQdXNY9JFXtuV8jdebywy/sAfFZEtubb3EJHtIrLLjQYgCphhjBkBtPRasEp5yZi4OEKDg3NsCw0OZkxcnI8iUqpgimOAexLwLvBZ5gYRsQLjgG5AMrBWROYCVuDlXMePAH4EZojICGByMcSslEdlDmI/vWgRB1JSqBsezpi4OB3cViWGmGJYUSoi9YB5xpjmGc/bAc8bY7pnPH8SwBiTu6HIPP5RYI0xZrmIzDDGDHCxzyhgFEDdunVb7d+/3xtvRSmlApaIrDfGxLp6zVdjFrWAg9meJ2dsy8s3wIMi8gGwz9UOxpjxwAvAhjJlyngoTKWUUlBC1lkYY7YAl/QmXOyXBCTFxsaO9H5USilVeviqZ3EIqJPtee2MbUWiuaGUUso7fNVYrAUiRKS+iJQBBgFzi3pSzQ2llFLeURxTZxOB1UAjEUkWkbuNMTbgfuBbYBvwhTHmVw9cS3sWSinlBcUyG6q4icgxIHM6VDiQvfXI/dzVtmrAca8FmHccnj7Wnf3y2qcg2/P7jAPl83Rn30D5PF1d1xvHFfbzzOu1wmzz98+zIMcW9fOsZIyp7vJVY0xAP4Dxl3uexz7rijsubxzrzn557VOQ7fl9xoHyebqzb6B8nkX5TIvj83T3s3Nnm79/ngU51tOfZ/ZHiUj3UURJ+TzPa5u3FeWa7h7rzn557VOQ7e58xt5WHJ+nO/sGyudZlOsWx+eZ12tF2eZt/vBvvqCfZ5aAvA1VVCKyzuSxMEUVnH6enqWfp2fp5+me0tCzKIzxvg4gwOjn6Vn6eXqWfp5u0J6FUkqpfGnPQimlVL60sVBKKZUvbSyUUkrlSxsLpZRS+dLGIh8iEiYin4rIRyIy1NfxBAKR/2/v3kKsKsMwjv+fupDKsjK8yAIphbTQsaiLSvGiJAlMsgOBdrACvSgIFCOiNAoRb0IK9cIwUzCQlMAKidIMEqzRNBM6WEQHiA52kk72dLE/87hnuW3P7JnZzw8G9qz9fd96eVkz715rw/fqIknLJa1tdSz9gaQp5fp8UdLEVsfT10kaKWmppLWSZrU6nt6iLYtFg61eb6bW0vV+YHKPB9tHNJJT23tt39uaSPuGBvO5vlyfM4HbWxFvb9dgPvfYngncBlzTinh7o7YsFtRavd5w+IHDWr1OAkYBd0gaRW379IONmg70YIx9zQpOPKdRbQWN5/PR8n4cawUN5FPSZGAD8ErPhtl7tWWxsP0W8MNRh68CPimfev8E1gA3Uevid0EZ05b5OhEN5jQqNJJP1SwEXrXd2dOx9gWNXp+2X7Y9Ccij5yL//A6p1+r1JWCqpCW0bo+evuq4OZU0uLTIHXuw/3qckHrX6APAdcAtkma2IrA+qt71OUHSYknLyJ3Ff/pEW9VWsv0bcE+r4+hPbH9P7fl6NIHtxcDiVsfRX9jeBGxqcRi9Tu4sDumWVq9tLjltruSzuZLPBqRYHNItrV7bXHLaXMlncyWfDWjLYtGTrV7bRXLaXMlncyWf/192nY2IiEpteWcRERGNSbGIiIhKKRYREVEpxSIiIiqlWERERKUUi4iIqJRiEdEkkgZIel3SDkmVW4VLulvSM+X1TEl3dlNcKyR9JWlA+f08SZ93x7mi/8reUBHNMxbAdkejE20vbX44RzgAzACWdPN5op/KnUW0BUnDJO0pHeV2S9oo6bTyXoekrZJ2Slon6ZyKtc6VtL6M3ypptKQhwCrgynJncfFRcx6U9GGZs+Y4a86TNLu8Hl7uUN6X1HlwLUlzJG0ra8wvx86QtKGM/aCLO5qngYck5QNinJQUi2gnI4BnbV8K7AOmluMrgbm2RwO7gMcr1pkPbC/jHwFW2v4WuA/YYrvD9qdHzXkYGFvmVO24u7rEOQa4GvhGtXapI6j1YOgArpA0nlpDn69tj7F9GfBanTW/AN4GplecO+K4UiyinXxme0d5/R4wTNIg4Gzbm8vx54HxFetcC7wAYPsNYLCksyrm7ARWS5oG/F1vkKQzgaG215X1f7e9H5hYfrYDncAl1IrHLuB6SQsljbP9UxcxLADmkL/7OAm5aKKd/HHY6wP07Hd2N1Jr4Xk5sO0kHgcJWFDuWjpsD7e93PZHZc1dwJOSHqu3gO2PgR3UektHNCTFItpa+ST+o6Rx5dB0YHMXUwC2UNptSpoAfGf753qDJZ0CXGj7TWAuMAgYWCeeX4AvJU0pcwdIOp3azqgzJA0sx4dKGiLpfGC/7VXAImqFoytPAbMrxkQcI192RcBdwNLyT3kvpTOipCeAd20f3eNgHvCcpJ3A/jK/K6cCq8ojLwGLbe+TVG/8dGBZOf9fwK22N0oaCbxT5v0KTAOGA4sk/VPGzuoqENu7JXVSXVQijpAtyiMiolIeQ0VERKUUi4iIqJRiERERlVIsIiKiUopFRERUSrGIiIhKKRYREVHpX1/hjeksA3EdAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Looking at the graph below, we can say that ATR uses more number of steps to achieve the required accuracy than RI. Making RI more efficient when greater accuracy is required since it takes less time for the computer to compute."
      ],
      "metadata": {
        "id": "-_cyxW2Bpo4v"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "m2a6sg90R-NX"
      },
      "source": [
        "### PROBLEM 3 (10 pts.)\n",
        "\n",
        "**Integration using Gaussian quadrature**\n",
        "\n",
        "Calculate the integral from Problem 1 using the Gaussian quadrature method: $ 2x^4 + 3x^2 + 4x + 5$ from $x=0$ to $x=5$.\n",
        "\n",
        "1. Define the function $\\tt gaussxw$ from Appendix E of the Newman text (and run the cell). *(2 pts.)*\n",
        "\n",
        "2. Calculate the integral using the $\\tt gaussxw$ function in (1) and number of slices $N = 3$. Describe each step in the calculation clearly and concisely. *(8 pts.; 4 pts. for the code & 4 pts. for the discussion text)*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YO_PDJ4VIVs5"
      },
      "source": [
        "#copying the gaussxw function straight from the book \n",
        "from numpy import ones,copy,cos,tan,pi,linspace\n",
        "\n",
        "def gaussxw(N):\n",
        "\n",
        "    # Initial approximation to roots of the Legendre polynomial\n",
        "    a = linspace(3,4*N-1,N)/(4*N+2)\n",
        "    x = cos(pi*a+1/(8*N*N*tan(a)))\n",
        "\n",
        "    # Find roots using Newton's method\n",
        "    epsilon = 1e-15\n",
        "    delta = 1.0\n",
        "    while delta>epsilon:\n",
        "        p0 = ones(N,float)\n",
        "        p1 = copy(x)\n",
        "        for k in range(1,N):\n",
        "            p0,p1 = p1,((2*k+1)*x*p1-k*p0)/(k+1)\n",
        "        dp = (N+1)*(p0-x*p1)/(1-x*x)\n",
        "        dx = p1/dp\n",
        "        x -= dx\n",
        "        delta = max(abs(dx))\n",
        "\n",
        "    # Calculate the weights\n",
        "    w = 2*(N+1)*(N+1)/(N*N*(1-x*x)*dp*dp)\n",
        "\n",
        "    return x,w\n",
        "\n",
        "\n",
        "def gaussxwab(N,a,b):\n",
        "  x,w = gaussxw(N)\n",
        "  return 0.5*(b-a)*x+0.5*(b+a),0.5*(b-a)*w\n",
        "  "
      ],
      "execution_count": 992,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "N = 3 # number of slices\n",
        "a = 0 #lower limit integration\n",
        "b = 5 #upper limit integration\n",
        "\n",
        "#calculating the sample points and weights \n",
        "x, w = gaussxw(N) #assigning gaussxw(N) values as x and w\n",
        "xp = 0.5 * (b-a) * x + 0.5 * (b+a) \n",
        "wp = 0.5 * (b-a) * w \n",
        "\n",
        "# Performing the integration\n",
        "s = 0.0\n",
        "for k in range(N):\n",
        "  s += wp[k] * f(xp[k]) #using f(x) from the first problem\n",
        "\n",
        "print(\"Using the gaussxw function, the integral is\" , s)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "OSD55s0TXJz_",
        "outputId": "5a2c1f9d-b4d1-4dc9-e7cd-143db5fe7218"
      },
      "execution_count": 993,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Using the gaussxw function, the integral is 1450.000000000002\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "HceOjNu-VfN9"
      },
      "source": [
        "### PROBLEM 4 (30 pts.)\n",
        "\n",
        "**Physics Application: Anharmonic Oscillator**\n",
        "\n",
        "_Refer to Exercise 5.10 of the Newman text._ Consider an anharmonic oscillator with potential $V(x)$ symmetric about $x = 0$. Let us set our anharmonic oscillator going with amplitude $a$. That is, at $t = 0$ we release it from rest at position $x = a$ and it swings back towards the origin. Then at $t = 0$ we have $dx/dt = 0$ and the equation for its total energy becomes $E = V(a)$, i.e., the total energy of the particle is given in terms of the amplitude.\n",
        "\n",
        "When the particle reaches the origin for the first time, it has gone through one quarter of a period of the oscillator. The period $T$ is given by: $$ T = \\sqrt{8m} \\int_0^a \\frac{dx}{\\sqrt{V(a)-V(x)}}. $$\n",
        "\n",
        "1. Suppose the potential is $V(x) = x^4$ and the mass of the particle is $m = 1$. Write a Python function that calculates the period of the oscillator for given amplitude a using Gaussian quadrature with $N = 20$ points. Describe each step clearly and concisely. *(16 pts.; 8 pts. for the code & 8 pts. for the discussion text)*\n",
        "\n",
        "2. Use your function to make a graph of the period for amplitudes ranging from $a = 0$ to $a = 2$. Describe the behavior you find and explain these results clearly and concisely. *(14 pts.; 7 pts. for the plot & 7 pts. for the discussion text)*"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "def gquadrature(function,a, b): # calculates position and weights and does the mapping \n",
        "\n",
        "    xp, wp = gaussxwab(N,a,b)\n",
        "\n",
        "    s = 0.0\n",
        "\n",
        "    for k in range(N):\n",
        "        s += wp[k]*function(xp[k])\n",
        "\n",
        "    return s\n"
      ],
      "metadata": {
        "id": "yUc_yB7nXHdv"
      },
      "execution_count": 994,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "p28JF-OQYock",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 404
        },
        "outputId": "cfb158d5-0574-4366-e184-8ac5328a6897"
      },
      "source": [
        "#4\n",
        "m = 1 #mass of the particle\n",
        "N = 20 #number of points\n",
        "a = 0 # initial position\n",
        "amplitude = linspace(0.1,2,N) # start, stop, number of samples to generate\n",
        "\n",
        "def k(x): #The function for period T V(a) equal to the amplitude raised to 4\n",
        "  return sqrt(8*m) / sqrt(amp**4 - x**4)\n",
        "\n",
        "S = [] #list to store the period of the oscillator given amplitude\n",
        "for amp in amplitude:\n",
        "  S.append(gquadrature(k,a,amp)) #adds the period of the oscillator to the list S\n",
        "\n",
        "#plotting the graph \n",
        "plt.figure(figsize=(8,6))\n",
        "plt.title(\"Period of an Anharmonic Oscillator\")\n",
        "plt.xlabel(\"Amplitude (m)\")\n",
        "plt.ylabel(\"Period (s)\")\n",
        "plt.plot(amplitude,S, 'salmon', linewidth = '4', alpha = 0.75)\n",
        "plt.grid(color = 'teal', linestyle = '--', linewidth = 0.5, alpha = 0.75)\n",
        "plt.show()\n",
        "\n"
      ],
      "execution_count": 995,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 576x432 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe4AAAGDCAYAAADtffPSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9e3hc513g//nO6G5ZtizZkS+KbMWxaZymbuO0SShpWvfiBkzLrSywbMtls7uwLCwsD/zY/tgCC8vy47YsLA+XdtOlLEuhBZq0pC1uWm/rtInTqkmcxIojx5GiTCIpsmVZsizNfH9/nCNpNLqNrDnvmffr9/M882jO7T3fz5yxv/Oe815EVQkEAoFAIOAHmbQDCAQCgUAgUD4hcQcCgUAg4BEhcQcCgUAg4BEhcQcCgUAg4BEhcQcCgUAg4BEhcQcCgUAg4BEhcQdMISL/KCLvv8pjnxeRt1/lsf9ZRIZFJHc1x7tCRD4gIl9OO46rQUROicjdaccxS3E8IvIhEflY/H63iKiI1KQaYMAsIXEHUidOmJMiMi4iL4vIfSLSfDVlqeq7VfWjlY5xJUTkeuDngJtUtcPROfeISEFE/tjF+aoBVT2gql9c63EiUi8i/0VEXoi/Z8+KyM+LiKQRT0lsd4vIwHrKCFx7hMQdqBaOqmoz8AbgEPDBtRwsEWl9n68HRlT1FYfn/BfAKPD9IlLv8LxzeFSj/BvgMHAPsBH4YeBe4L+lGVQl8OgaBCpISNyBqkJVXwT+EbgZQERuF5ETInJeRL5ZfKtURL4oIr8uIl8BJoDueN2Px9szIvJBETknIq+IyP8SkU1Fx/9wvG1ERP7jSnGJyKb4+KH4mA/G5b8d+DywI75jcN8Sx7aKyAPxsaPx+10lHr8mIl8RkYsi8jkRaV8hFiFK3B8EpoGjJdtVRP51XLM8LyJ/VFq7FJHfjmM5KyLvLlr/IyLydBxHn4j8q6Jtd4vIgIj8QvxI4H/Gt4j/RkQ+Fh/zhIjsE5H/J/7M+0XknUVl7BCRT4nIqyJyRkT+ZdG2D4nIx+PP+WJ8K/pQ0fa5RxkikhWRXxKR5+J9HxORziU+q8PAO4HvUdUnVXVGVb8K/HPgJ0Vkb7zfB2Lfi/Fn8kNFZfzLos/kKRF5Q2k8K7HcZyoiG4i+67PfnfH486kXkd8XkcH49fuzP86WugarnT9gj5C4A1VF/J/vPcA3RGQn8GngPwNbgP8AfEJEthYdMlt72gicKynuA/HrrUA30Az8YXyem4A/jo/fAbQBu1ie/w5sist5C1Hi/BFV/Sfg3cCgqjar6geWODZD9B9sF1HtfHI2jiJ+EPgRYBtQF7sux5vjWP8P8HFgqWf63wHcBtwCvA94V9G2NwGngXbgt4APFyX2V+JjW+J4fm82UcV0EF2LLqLPHaIfDn8BtALfAD4bO+8EfhX4k6Lj/w8wQPSZfy/wGyLytqLt3xnvsxn4FIs/p1l+FvgBou9KC/CjRD/eSnkH8DVV7S9eqapfi+M4HCfQPwDeraobgTuBHgAR+T7gQ0TXuyWOb2SZmJZjyc9UVS+x8LvTrKqDwH8EbgcOAq8D3sjCO1BLXYPAtYSqhld4pfoCngfGgfNEyfd/AI3ALwB/UbLvZ4H3x++/CPxqyfYvAj8evz8G/ETRtv1ENdQa4JeB/1O0bQNwBXj7EvFl4203Fa37V8AX4/d3AwNr8D0IjJbE/MGi5Z8AHlzh+D8H/j5+f0fstK1ouwJvLlr+OPCL8fsPAGeKtjXF+3csc66/B366yPMK0FC0/UPA54uWj8bXMhsvb4zL3wx0AnlgY9H+/wW4r6isfyradhMwWfI9eXv8/jTwnjI+6z8vvs4l275KlCQ3xN+97wEal/i+/fQK39u3F8X+sfj97ti5pszPdKBk+3PAPUXL7wKeX+4ahNe19wo17kC18F5V3ayqXar6E6o6SVSj+L74du95ETlPVNvcXnRc/5KlRexgYS38HFHSvi7eNnesRrWf5WpS7UDtEmXtLEdMRJpE5E/iW+xjwHFgs4hki3Yrbo0+QXR3YKmyGoHvA/4yjvth4AWiGnsxK5U3t01VZ2upzXH57xaRr8a3ss8T1WiLb9sPqerlknO9XPR+EhhW1XzR8mz5O4BXVfVi0f6ln2Np3A2y9HPcTqIEtxrDLPy+FLM9jvUS8P3AvwZeEpFPi8i3rPE8y1LGZ1rKUt/bHUXLS12DwDVESNyBaqafqMa9uei1QVV/s2iflaa3GyRK/rNcD8wQJZqXiP5TBqLkSnS7fCmGiWq1pWW9WKbHzxHV9t+kqi3AXbOnLfP4Yr6L6Jbr/xCRXPyccydL3y5fE/Fz1E8Avw1cp6qbgc+UxLme6QQHgS0isrFo3Vo+x2L6gRvK2O+fgDeVPv8WkTcRXf8vAKjqZ1X1HUTJ/Bngz9Z4niUp4zNd6vNc6ns7WLQcpnS8xgmJO1DNfAw4KiLvihsjNcSNc1Z6Fl3MXwH/XqKuU83AbwB/raozwN8C3yEibxaROqJnsUv+e4hrjx8Hfl1ENopIF9Ez1o+VGcdGoprneRHZAvynMo9bivcDHwFeS3TL/SDwrcDrROS16ygXomfr9cAQMBM3WnvnyoeUj0bPmU8A/yW+lrcAP0b5n2Mxfw78mojcKBG3iMiiH14atUE4RtQ24kD8Pbo9Pucfq+qzInKdiLwnftY9RXSrv1B0nv8gIrfG59kbX/9yWe0zfRlok6JGk0Tf2w+KyFaJGin+Mlf3GQWMEhJ3oGqJ/6N/D/BLRP/x9QM/T/nf248QNZo6DpwFLgM/FZd9CvhJ4H8T1b5HiRorLcdPAZeAPuDL8XEfKTOO3yd6Zj9M9Fz1wTKPW0DcWO8w8Puqmit6PRaXua5ad3wL+98R/UgZJbr9/qn1lLkEP0D0DHgQ+DvgP8XJda38LlGcnwPGgA8TfcZL8T3AQ0Sf0ThREvww8XeB6Pv0s3FMrxI1Pvw3AKr6N8CvE13vi0TPp7eUG+Rqn6mqPkOUqPvix0E7iBpjngQeB54Avh6vCwQAENVw1yUQCAQCAV8INe5AIBAIBDwiJO5AIBAIBDwiJO5AIBAIBDwiJO5AIBAIBDwiJO5AIBAIBDzCi5ll2tvbdffu3WmHUTGUqxt5o9qx6BWc/MGil0UnsOlVaafHHntsWFW3LrXNi8S9e/duTp48mXYYFeP4uXPc1bWWMRz8wKJXcPIHi14WncCmV6WdRKR00qQ5wq3yFOgbHU07hESw6BWc/MGil0UnsOnl0smLAVgOHTqklmrcIxMTtDU1pR1GxbHoFZz8waKXRSew6VVpJxF5TFUPLbUt1LhT4P7e3rRDSASLXsHJHyx6WXQCm14unULiToGW+vq0Q0gEi17ByR8sell0ApteLp1C4k6Bgx0daYeQCBa9gpM/WPSy6AQ2vVw6hcSdAsfPLdtY0GssegUnf7DoZdEJbHq5dAqJOwUs/toEm17ByR8sell0ApteocZtnJGJibRDSASLXsHJHyx6WXQCm14unULiToH+sbG0Q0gEi17ByR8sell0ApteLp1CP+4UsNiHEWx6BSd/sOhl0QlseoV+3AmgqhTO9VF48hvkH/4S+c/djxYKqcRisQ8j2PQKTv5g0cuiE9j0cunkxVjllUBEyH/1S3B5cn7lpXHY2OI8lrbGRufndIFFr+DkDxa9LDqBTS+XTtdMjRtAmkuS9MULqcSxv709lfMmjUWv4OQPFr0sOoFNL5dO11TiZuOmBYt6MZ0GEif6+1M5b9JY9ApO/mDRy6IT2PRy6XRNJW4puS2eVuK+bceOVM6bNBa9gpM/WPSy6AQ2vVw6XVOJm5aFNe60bpVb7AoBNr2Ckz9Y9LLoBDa9XDpdU4l7cY07ncSdGx9P5bxJY9ErOPmDRS+LTmDTy6XTNZW4S59xc3GMNPqxH923z/k5XWDRKzj5g0Uvi05g08ul07WVuBsaoaZ2fnlmGqYuOw/DYh9GsOkVnPzBopdFJ7DpFebjTggRgdIuYWPub5d3NDc7P6cLLHoFJ3+w6GXRCWx6uXS6phI3VMdz7s4W94O+uMCiV3DyB4teFp3AppdLp2sucS9uWe6+deOjg4POz+kCi17ByR8sell0ApteLp0SS9wi0iAij4jIN0XklIj8Srz+PhE5KyI98etgUjEsGVcV1Ljv7Ox0fk4XWPQKTv5g0cuiE9j0cumUZI17Cnibqr4OOAgcEZHb420/r6oH41dPgjEsZomW5a45PTzs/JwusOgVnPzBopdFJ7Dp5dIpsUlGNOpnNduxrTZ+pT6HqCwa9tR9jXtkcnL1nTzEoldw8geLXhadwKaXS6dE5+MWkSzwGLAX+CNV/QURuQ+4g6hGfgz4RVWdWuLYe4F7Abbv3HnrbzzwAHd1ddGTyzE2NcXRffu4v7eX7tZWGmtqODU0xOE9e3h4YIDpfJ7D3d08eOYM+9raAOgdGeHI3r184bkzvOXEP7Gpro7Ry5fZUFvLU4ePcubi+FyZLfX1HOzo4Pi5cxzs6GBkYoL+sbG57W2Njexvb+dEfz+37dhB/9gYufH54zuam+lsaeHRwUHu7Ozk9PAwI5OTc9s31dezp7WVnlyuIk7H+vqozWa5Y9cujp09y4GtW5mcmaFvdNSZU2dLC1kRnr9wITgFJ+dObU1NfPmFF7jnxhuDU5U79eRyvHbbNs5duBCcVnC6fvPmZefjTjRxz51EZDPwd8BPASNADqgD/hR4TlV/daXjDx06pCdPnqxYPDOf/NiCW+TZo+9Dtrib2eW+nh4+cNDpo30nWPQKTv5g0cuiE9j0qrSTiCybuJ20KlfV88BDwBFVfUkjpoD/CbzRRQzFSHO6s4RZ7AoBNr2Ckz9Y9LLoBDa9THQHE5GtcU0bEWkE3gE8IyLb43UCvBd4MqkYliXlyUbampqcns8VFr2Ckz9Y9LLoBDa9XDolWePeDjwkIo8DjwKfV9UHgL8UkSeAJ4B24D8nGMOSpD29Z08u5/R8rrDoFZz8waKXRSew6eXSKclW5Y8Dr19i/duSOmfZlCRu1zXuu7q6nJ7PFRa9gpM/WPSy6AQ2vVw6XXsjp7FUl7BQ464EFr2Ckz9Y9LLoBDa9XDpdk4l7UY370kW0kHd2+rGpRb3fTGDRKzj5g0Uvi05g08ul0zWZuKWmFhqLGhKogsNJ0C3ORQs2vYKTP1j0sugENr3CfNwOkOb0xiy3OBct2PQKTv5g0cuiE9j0CvNxuyDFMcu7W1udncslFr2Ckz9Y9LLoBDa9XDpds4lbWtIbs7yxJrHG/Kli0Ss4+YNFL4tOYNPLpdM1m7gXdwlzV+M+NTTk7FwusegVnPzBopdFJ7Dp5dLpmk3cac7LfXjPHmfncolFr+DkDxa9LDqBTS+XTtds4l7qGbeLCVcAHh4YcHIe11j0Ck7+YNHLohPY9HLpdO0m7voGqK2bX87PwOSEk1NP5931GXeJRa/g5A8WvSw6gU0vl07XbOIWEWhO5zn34e5uJ+dxjUWv4OQPFr0sOoFNL5dO12zihvSecz945oyT87jGoldw8geLXhadwKaXS6drOnGn1Zd7X1ubk/O4xqJXcPIHi14WncCml0unazpxp9mXOxAIBAKBq+HaTtwp9eXuHRlxch7XWPQKTv5g0cuiE9j0cul0TSfu0kFYXNW4j+zd6+Q8rrHoFZz8waKXRSew6eXS6dpO3BuaIVP0EUxdRqevJH7aY319iZ8jDSx6BSd/sOhl0Qlserl0uqYTt0hmcZewseRr3bXZbOLnSAOLXsHJHyx6WXQCm14una7pxA1LdQlL/jn3Hbt2JX6ONLDoFZz8waKXRSew6eXS6ZpP3Iu7hCVf4z529mzi50gDi17ByR8sell0ApteLp2u+cSdxiAsB7ZuTfwcaWDRKzj5g0Uvi05g08ul0zWfuNOY3nNyZibxc6SBRa/g5A8WvSw6gU0vl07XfOKWje4HYekbHU38HGlg0Ss4+YNFL4tOYNPLpdM1n7gXtSq/dAlNeJaXo/v2JVp+Wlj0Ck7+YNHLohPY9HLpdM0nbqmpgaYNRWsUxpO9XX5/b2+i5aeFRa/g5A8WvSw6gU0vl07XfOIGkObS2+XJJu6W+vpEy08Li17ByR8sell0ApteLp1C4gZoKW2gluxz7oMdHYmWnxYWvYKTP1j0sugENr1cOoXEzVIN1JKtcR8/dy7R8tPColdw8geLXhadwKaXS6eQuMF5lzCLvzbBpldw8geLXhadwKZXqHE7xvUgLCMTE4mWnxYWvYKTP1j0sugENr1cOoXEDUsMezqGqiZ2uv4xN/N+u8aiV3DyB4teFp3AppdLp5C4AalvgLqiFoGFPExcSux8Fvswgk2v4OQPFr0sOoFNr9CPOwWkdCCWBJ9zW+zDCDa9gpM/WPSy6AQ2vUz04xaRBhF5RES+KSKnRORX4vV7RORrInJGRP5aROqSimFNOHzO3dbYmFjZaWLRKzj5g0Uvi05g08ulU5I17ingbar6OuAgcEREbgf+K/B7qroXGAV+LMEYykZa3E3vub+9PbGy08SiV3DyB4teFp3AppdLp8QSt0aMx4u18UuBtwF/G6//KPDepGJYEw77cp/o70+s7DSx6BWc/MGil0UnsOnl0inRZ9wikhWRHuAV4PPAc8B5VZ2d/2wA2JlkDOVS2iUsyWfct+3YkVjZaWLRKzj5g0Uvi05g08ulU02ShatqHjgoIpuBvwO+pdxjReRe4F6A7Tt3cl9PD3d1ddGTyzE2NcXRffu4v7eX7tZWGmtqODU0xOE9e3h4YIDpfJ7D3d08eOYM+9raAOgdGeHI3r0c6+ujNpvljl27OHb2LAe2bmVyZoYXc4O8N19geHKCGsmwQYb5ZE8PBzs6GJmYoH9sbO6cbY2N7G9v50R/P7ft2EH/2Bi58fG57R3NzXS2tPDo4CB3dnZyeniYkcnJue0vjY8zXSjQk8sl6tQ3OjpXZkt9PQc7Ojh+7lwiTp0tLTwzPMyjg4PBKTg5d2prauIvvvlNfvKNbwxOVe7Uk8uhqpweGQlOKzitmB+T7K+84EQivwxMAr8AdKjqjIjcAXxIVd+10rGHDh3SkydPJhqfqpL/2J9GXcFisv/sR6OuYhXmvp4ePnDwYMXLTRuLXsHJHyx6WXQCm16VdhKRx1T10FLbkmxVvjWuaSMijcA7gKeBh4DvjXd7P/APScWwFkTE2dCnFvswgk2v4OQPFr0sOoFNLyv9uLcDD4nI48CjwOdV9QGiGvfPisgZoA34cIIxrInSvtxJNVCz2IcRbHoFJ3+w6GXRCWx6uXRK7Bm3qj4OvH6J9X3AG5M677pYVONOpktYR3NzIuWmjUWv4OQPFr0sOoFNL5dOYeS0Ikr7cic1CEtn6fzfRrDoFZz8waKXRSew6eXSKSTuYpaYbCQJHh0cTKTctLHoFZz8waKXRSew6eXSKSTuIhZP75lM4r6zszORctPGoldw8geLXhadwKaXS6eQuItpbgFkfnniEpqfWXb3q+X08HDFy6wGLHoFJ3+w6GXRCWx6uXQKibsIyWZhw4aiNZrI7fKRycmKl1kNWPQKTv5g0cuiE9j0cukUEncJ4mDMcot9GMGmV3DyB4teFp3AppeVftx+4mAQFot9GMGmV3DyB4teFp3AppeJ+bh9ZXGNu/Jdwix2hQCbXsHJHyx6WXQCm16hO1iaOBiEpa2pqeJlVgMWvYKTP1j0sugENr1cOoXEXYKLZ9w9uVzFy6wGLHoFJ3+w6GXRCWx6uXQKibuU0hr3+BiqhYqe4q6uroqWVy1Y9ApO/mDRy6IT2PRy6RQSdwlSVw/FU3kWCnDpUkXPYfHXJtj0Ck7+YNHLohPY9Ao17pQpnSWs0s+5x6amKlpetWDRKzj5g0Uvi05g08ulU0jcS7FospHKPue22IcRbHoFJ3+w6GXRCWx6hX7cKVM6Znmla9wW+zCCTa/g5A8WvSw6gU2v0I87bRJuWd7d2lrR8qoFi17ByR8sell0ApteLp1C4l6CpGcJa6ypqWh51YJFr+DkDxa9LDqBTS+XTiFxL8WiebkvoKoVK/7U0FDFyqomLHoFJ3+w6GXRCWx6uXQKiXspGpsgW/TrafoKTF2uWPGH9+ypWFnVhEWv4OQPFr0sOoFNL5dOIXEvgYgkOtnIwwMDFSurmrDoFZz8waKXRSew6eXSKSTuZSjty13JyUam8/mKlVVNWPQKTv5g0cuiE9j0cukUEvdytJQ+565cjftwd3fFyqomLHoFJ3+w6GXRCWx6uXQKiXsZFrcsr1yN+8EzZypWVjVh0Ss4+YNFL4tOYNPLpVNI3MuxqGV55Wrc+9raKlZWNWHRKzj5g0Uvi05g08ulU0jcy5B0X+5AIBAIBK6GkLiXo3kjIPPLk5fQmZmKFN07MlKRcqoNi17ByR8sell0ApteLp1C4l4GyWShuXnhygo95z6yd29Fyqk2LHoFJ3+w6GXRCWx6uXQKiXsFpDmZMcuP9fVVpJxqw6JXcPIHi14WncCml0unkLhXIqFBWGqz2YqUU21Y9ApO/mDRy6IT2PRy6RQS9wrIonm5K3Or/I5duypSTrVh0Ss4+YNFL4tOYNPLpVNI3CuR0Lzcx86erUg51YZFr+DkDxa9LDqBTS+XTiFxr4AkNC/3ga1bK1JOtWHRKzj5g0Uvi05g08ulU0jcK1Fa4x6/iBYK6y52skLdyqoNi17ByR8sell0ApteLp1C4l4Bqa2Dhsb5FVqAS+PrLrdvdHTdZVQjFr2Ckz9Y9LLoBDa9XDollrhFpFNEHhKRp0TklIj8dLz+QyLyooj0xK97koqhEpTeLq/Ec+6j+/atu4xqxKJXcPIHi14WncCml0unJGvcM8DPqepNwO3AT4rITfG231PVg/HrMwnGsH4SGPr0/t7edZdRjVj0Ck7+YNHLohPY9HLpVJNUwar6EvBS/P6iiDwN7EzqfEkhG1vQouVKJO6W+vp1l1GNWPQKTv5g0cuiE9j0cunk5Bm3iOwGXg98LV71b0XkcRH5iIi0uojhqkngVvnBjo51l1GNWPQKTv5g0cuiE9j0cumUWI17FhFpBj4B/IyqjonIHwO/Bmj893eAH13iuHuBewG279zJfT093NXVRU8ux9jUFEf37eP+3l66W1tprKnh1NAQh/fs4eGBAabzeQ53d/PgmTNzU631joxwZO9ejvX1UZvNcseuXRw7e5YDW7cyOTND3+joXJkt9fUc7Ojg+LlzvKmuhvbJy1zOz9De2MTzz/dx5vRp9re3c6K/n9t27KB/bIzc+Pjc8R3NzXS2tPDo4CB3dnZyeniYkcnJue1PDQ3xg699LT25XCpOBzs6GJmYoH9sbG57W2Pjupw6W1r4Sn8/12/aFJyCk3OntqYm/vCRR/ilb/u24FTlTj25HK9cusTeLVuC0wpOK+ZVVV1xh/UgIrXAA8BnVfV3l9i+G3hAVW9eqZxDhw7pyZMnE4lxNXRygvzH75tfUVNL9gd/HBFZ9pjV6MnlTP7itOgVnPzBopdFJ7DpVWknEXlMVQ8ttS3JVuUCfBh4ujhpi8j2ot2+C3gyqRgqQkMj1NTOL89Mw+XJdRU5MjGxzqCqE4tewckfLHpZdAKbXi6dkrxV/q3ADwNPiEhPvO6XgB8QkYNEt8qfB/5VgjGsGxGB5hY4XzTX6sUxaGy66jL7xyozAlu1YdErOPmDRS+LTmDTy6VTorfKK0Wat8oB8g89iL4wP2Vb5s2Hydyw/6rLG5mYoK3p6hN/tWLRKzj5g0Uvi05g06vSTqncKjdFhScbsdiHEWx6BSd/sOhl0Qlserl0Com7DKTCg7C0NTauvpOHWPQKTv5g0cuiE9j0cukUEnc5LOrLvb7Evb+9fV3HVysWvYKTP1j0sugENr1cOoXEXQaLp/dc363yE/396zq+WrHoFZz8waKXRSew6eXSKSTuctjQDFL0UV2eRKevXHVxt+3YUYGgqg+LXsHJHyx6WXQCm14unULiLgPJZKC5eeHKddwut9gVAmx6BSd/sOhl0Qlserl0Com7TKS59Hb51V+k3Pj65/SuRix6BSd/sOhl0Qlserl0Com7XFoqN9mIxblowaZXcPIHi14WncCml5X5uE1RyS5hFvswgk2v4OQPFr0sOoFNr9CPuxqp4CAsHaXPy41g0Ss4+YNFL4tOYNPLpVNI3GWyuEvY1de4O1taVt/JQyx6BSd/sOhl0Qlserl0Com7XEpr3JcuooX8VRX16OBgBQKqPix6BSd/sOhl0Qlserl0Com7TKSmduGMYKpwla0I7+zsrFBU1YVFr+DkDxa9LDqBTS+XTiFxrwFpLm2gdnXPuU8PD1cinKrDoldw8geLXhadwKaXS6eQuNdChcYsH5mcrEAw1YdFr+DkDxa9LDqBTS+XTiFxrwFpqcyY5Rb7MIJNr+DkDxa9LDqBTa/Qj7taWdQl7Opq3Bb7MIJNr+DkDxa9LDqBTa/Qj7tKqdQsYRa7QoBNr+DkDxa9LDqBTa/QHaxaWaLGraprLqatqWn1nTzEoldw8geLXhadwKaXS6eQuNdCfQPU1s0v52dgcmLNxfTkchUMqnqw6BWc/MGil0UnsOnl0ikk7jUgIhVpWX5XV1eFIqouLHoFJ3+w6GXRCWx6uXQKiXuNLJ5sZO3PuS3+2gSbXsHJHyx6WXQCm16hxl3FlCbuq6lxj01NVSia6sKiV3DyB4teFp3AppdLp5pydhKRbcC3AjuASeBJ4KSqFhKMrTqpQMtyi30YwaZXcPIHi14WncCmV9X04xaRt4rIZ4FPA+8GtgM3AR8EnhCRXxERe+36V2BRjXts7YnbYh9GsOkVnPzBopdFJ7Dp5dJptRr3PcC/VNUXSjeISA3wHcA7gE8kEFt1UlrjHl/7rfLu1tZKRVNVWPQKTv5g0cuiE9j0cum0YuJW1Z9fYdsM8PcVj6ja2bABMhkoxE8Jpi6jV6aQuvqyi2isKesJhXdY9ApO/mDRy6IT2PRy6VRW4zQR+WkRaZGID4vI10XknUkHV42IZKB5fQ3UTg0NVTCi6sGiV3DyB4teFp3AppdLp3Jblf+oqo4B7wRagR8GfjOxqKqc9XYJO7xnTyXDqUibb/QAACAASURBVBosegUnf7DoZdEJbHq5dCo3cUv89x7gL1T1VNG6a491DsLy8MBABYOpHix6BSd/sOhl0Qlserl0KjdxPyYinyNK3J8VkY3AtdcVLGa9Ne7pfL6S4VQNFr2Ckz9Y9LLoBDa9XDqV+zT9x4CDQJ+qTohIG/AjyYVV5axzEJbD3d0VDKZ6sOgVnPzBopdFJ7Dp5dJptX7cuwFUtaCqX1fV8/HyiKo+HjdW25V8mNXFeqf3fPDMmUqGUzVY9ApO/mDRy6IT2PRy6bRajfv/E5EM8A/AY8AQ0ADsBd4KHAb+E2DvgcVKlLYqv3QJzeeRbLasw/e1tSUQVPpY9ApO/mDRy6IT2PRy6bRaP+7vE5GbgB8CfpRo5LQJ4GngM8Cvq+rlpY4VkU7gfwHXAQr8qar+NxHZAvw1sBt4Hnifqo5WxMYRUlMDTRtg4lK8RmF8DDbZG1QgEAgEAtXFqo3TVPUpVf2Pqnq3qu5X1der6g+q6seWS9oxM8DPqepNwO3AT8Y/An4ROKaqNwLH4mXvkObS2+XlP+fuHRmpdDhVgUWv4OQPFr0sOoFNL5dOic0OpqovqerX4/cXiWrpO4H3AB+Nd/so8N6kYkiUltIGauU/5z6yd2+Fg6kOLHoFJ3+w6GXRCWx6uXRyMkZb3Mjt9cDXgOtU9aV4U47oVvpSx9wL3AuwfedO7uvp4a6uLnpyOcampji6bx/39/bS3dpKY00Np4aGOLxnDw8PDDCdz3O4u5sHz5yZe+7QOzLCkb17OdbXR202yx27dnHs7FkObN3K5MwMfaOjc2W21NdzsKOD4+fOcbCjg5GJCfrHxua2tzU28gbJUhgfp6Wunqn8DGeefor9XXu5v7eXjuZmOltaeHRwkDs7Ozk9PMzI5OTc8aeHh/n+m2+mJ5erKqf97e2c6O/nth076B8bIzc+Prd9NafOlha+NjDAjpaW4BScnDu1NTXxx48+yi+8+c3BqcqdenI5Ricn6dq8OTit4LRiTlXVteTgNSMizcCXiJ6Hf1JEzqvq5qLto6q64sPhQ4cO6cmTJxONc60Uzj5L4fjn55Zl126yh+8p69iPnzrF+w4cSCq01LDoFZz8waKXRSew6VVpJxF5TFUPLbVtxRq3iLxhpe2zt8JXOL6WaOawv1TVT8arXxaR7ar6kohsB15ZqYxqZT2DsKz2a8pXLHoFJ3+w6GXRCWx6uXRa7Rn378SvPyK6zf2nwJ/F7/9opQNFRIAPA0+r6u8WbfoU8P74/fuJupr5xxLDnpZ79+LY2bMJBJQ+Fr2Ckz9Y9LLoBDa9XDqtmLhV9a2q+lbgJeANqnpIVW8lel794iplfyvRZCRvE5Ge+HUP0eQk7xCRZ4G34+lkJVLfAMVTeRbyRd3DVubA1q0JRZUuFr2Ckz9Y9LLoBDa9XDqV2zhtv6o+Mbugqk+KyGtWOkBVv8zyE5EcLvO8VY00t6CvFk3ldnEMNjSvetzkzEyCUaWHRa/g5A8WvSw6gU0vl07ldgd7XET+XETujl9/BjyeZGBecJXPuftGvRpvpmwsegUnf7DoZdEJbHq5dCq3xv0jwL8BfjpePg78cSIReYS0bGLBU+0yE/fRffsSiSdtLHoFJ3+w6GXRCWx6uXQqq8Ydj5D2R8AvA/8v8IerjJp2bbBospHyRk+7v7c3iWhSx6JXcPIHi14WncCml0unsmrcInI30ShnzxM9t+4Ukfer6vHkQqt+SmcJK3d6z5b6+tV38hCLXsHJHyx6WXQCm14uncq9Vf47wDtV9TSAiOwD/gq4NanAvOAqn3Ef7OhIIprUsegVnPzBopdFJ7Dp5dKp3MZptbNJG0BVe4HaZELyiKYNkCmayvPKFDq1+hOE4+fOJRhUelj0Ck7+YNHLohPY9HLpVG6N+6SI/DnwsXj5h4DqGoM0BUQkqnVfKGpNeHEM6htWPM7ir02w6RWc/MGil0UnsOlVjTXufwM8Bfy7+PVUvO6ap/Q5dzkN1EYmJpIKJ1UsegUnf7DoZdEJbHq5dCqrxq2qU8Dvxq9AMRvXPr1n/1j5c3f7hEWv4OQPFr0sOoFNL5dOq00y8nFVfZ+IPAEsGohbVW9JLDJPkI0L+3Lri/1wy8pt9iz2YQSbXsHJHyx6WXQCm17V1I97dsCV7wCOLvG65pHtOxcs6yuD6PlXVzzGYh9GsOkVnPzBopdFJ7Dp5dJptUlGXhKRLHCfqp4rfTmKsaqRzVuQbTsWrCucPrXiMW2NjUmGlBoWvYKTP1j0sugENr1cOq3aOE1V80BBRDattu+1inzLwsnT9bnT6PT0svvvb29POqRUsOgVnPzBopdFJ7Dp5dKp3Fbl48ATIvJhEfmD2VeSgfmEdHVDQ9Gvrekr6Nlnl93/RH+/g6jcY9ErOPmDRS+LTmDTy6VTuf24Pxm/AksgmSyZG19D4Ymvz60rnD6F3PiaqK93Cbft2LFonQUsegUnf7DoZdEJbHq5dCp3kpGPAh8HvqqqH519JRuaX8i+AyyYfvzVIRh+Zcl9LXaFAJtewckfLHpZdAKbXi6dykrcInIU6AEejJcPisinkgzMN6R5I7Kra8G6wuknl9w3Nz7uIiTnWPQKTv5g0cuiE9j0culU7jPuDwFvBM4DqGoP0J1QTN4i+0saqT1/Zsmxyy32YQSbXsHJHyx6WXQCm17V1I97lmlVLR0SrFDpYHxHdnZC88b5Ffk8euaZRftZ7MMINr2Ckz9Y9LLoBDa9qqYfdxGnROQHgayI3Cgi/x04kWBcXiKSIbNvYa27cPoUqgsHnetobnYZljMsegUnf7DoZdEJbHq5dCo3cf8UcACYAv43cAH4maSC8hm58TWQKfpYL15AXxpYsE9nS8n45kaw6BWc/MGil0UnsOnl0mnFxC0iDSLyM8BvAS8Ad6jqbar6QVVdfeLpaxBpaES6bliwTktGUnt0cNBlSM6w6BWc/MGil0UnsOnl0mm1GvdHgUPAE8C7gd9OPCIDZEobqfWfRSfmWxze2dnpOiQnWPQKTv5g0cuiE9j0cum0WuK+SVX/uar+CfC9wF0OYvKfbdth85b5ZVW096m5xdPDwykElTwWvYKTP1j0sugENr1cOq2WuOcG3FbVmYRjMYOIkNl/84J1hd6n0UIegJHJyTTCShyLXsHJHyx6WXQCm14unaS0xfOCjSJ54NLsItAITMTvVVWdPI0/dOiQnjx50sWpKoZOXyH/8Y/CzPxkI5m7j5Dp6mZkYoK2pqb0gksIi17ByR8sell0AptelXYSkcdU9dBS21ab1jOrqi3xa6Oq1hS9t9cssIJIbR3SvbBDvj4TjaRmsQ8j2PQKTv5g0cuiE9j0qsZ+3IGrYFEjtdwAOnbeZFcICF08fMGiE9j0sugENr2qpjtYYH3IlnZka8eCdYXTp8zdIprFoldw8geLXhadwKaXS6eQuBNm0fjlZ57h8RdfTCmaZOnJ5dIOoeIEJ3+w6GXRCWx6uXQKiTthZPcNUN8wv+LKFG/PLN8g0Gfu6upafSfPCE7+YNHLohPY9HLpFBJ3wki2Btn7LQvWXTr1jZSiSZbwK9oPLDqBTS+LTmDTK9S4jVHaSK3hwig6MpRSNMkxNjWVdggVJzj5g0Uvi05g08ulU0jcDpCNm5Ad88PhtTc2UTj9ZIoRJUOYY9cPLDqBTS+LTmDTqxrn414zIvIREXlFRJ4sWvchEXlRRHri1z1Jnb/akKKR1IYnJ9C+Z9Ertn51hr6ZfmDRCWx6WXQCm15W+nHfBxxZYv3vqerB+PWZBM9fVciuLmiK5mttrKmB/Az63OmUo6os3a2taYdQcYKTP1j0sugENr1cOiWWuFX1OPBqUuX7hmQyZPbdBEBGBIj6dK805KxvNNbUpB1CxQlO/mDRy6IT2PRy6ZTGp/dvReRfACeBn1PV0aV2EpF7gXsBtu/cyX09PdzV1UVPLsfY1BRH9+3j/t5eultbaayp4dTQEIf37OHhgQGm83kOd3fz4Jkz7GtrA6B3ZIQje/dyrK+P2myWO3bt4tjZsxzYupXJmRn6Rkfnymypr+dgRwfHz53jYEcHIxMT9I+NzW1va2xkf3s7J/r7uW3HDvrHxsiNj89t72huprOlhUcHB7mzs5PTw8OMTE7ynZ27Gf/yQ0xeuUJGhIsv9sOzz/BoHm+dZrd3trTwlf5+Tg0NeX+dgpN/Tm1NTfzJY4/R1tQUnKrcqSeX45VLl+gfGwtOKzitxIqTjKwXEdkNPKCqN8fL1wHDgAK/BmxX1R9drRwfJxlZjvwXP8vkc6dpiH+dSdcNZO9+V8pRVYb+Cxfo3LQp7TAqSnDyB4teFp3Aplelna56kpFKo6ovq2peVQvAnwFvdHn+aiCz/8CCbgP6wll0ciLFiCrHwwMDaYdQcYKTP1j0sugENr1cOjlN3CKyvWjxuwB7faJWo2MnlxqKxrTVAvrs0+nFU0Gm8/m0Q6g4wckfLHpZdAKbXi6dkuwO9lfAw8B+ERkQkR8DfktEnhCRx4G3Av8+qfNXKyJC28HbFqwr9J4iugnhN4e7u9MOoeIEJ3+w6GXRCWx6uXRKslX5D6jqdlWtVdVdqvphVf1hVX2tqt6iqt+pqi8ldf5q5jMFgWxRu8BL4+jAufQCqhAPnjmTdggVJzj5g0Uvi05g08ulUxg5LQW6r+tA9uxdsE6fOZVSNJVjtoWkJYKTP1j0sugENr1cOoXEnRKZopHUAHSwH704llI0gUAgEPCFkLhToHdkBGnfhrRtK1qrFE77XevuHRlJO4SKE5z8waKXRSew6eXSKSTuFDiyN7pNLiWzhumZp1GPW1vOelkiOPmDRS+LTmDTy6VTSNwpcKyvDwDZcyPU1c9vmLqMnnsupajWz6yXJYKTP1j0sugENr1cOoXEnQK12SwAUlOD3LB/wTZ9xt+u7bNelghO/mDRy6IT2PRy6RQSdwoUj0ObKb1dPpRDR/18/rPa+Lo+Epz8waKXRSew6eXSKSTuFDh29uzce9nUinTsXLDd10ZqxV5WCE7+YNHLohPY9HLpFBJ3ChzYunXBspR2DXvuNDp9xWVIFaHUywLByR8sell0ApteLp1C4k6ByZmZBcty/W5o3DC/YmYafa7XbVAVoNTLAsHJHyx6WXQCm14unULiToG+0YVTkEsmS+bG1yxYVzh9iiSnXE2CUi8LBCd/sOhl0Qlserl0Cok7BY7u27donex7DSDzK86PwFDOXVAVYCkv3wlO/mDRy6IT2PRy6RQSdwrc37v4Nrhs2Ih07l6wruDZ+OVLeflOcPIHi14WncCml0unkLhToKW+fsn1i0ZSO3cGvTzpIqSKsJyXzwQnf7DoZdEJbHq5dAqJOwUOdnQsuV52dMLGlvkVhQJ65hlHUa2f5bx8Jjj5g0Uvi05g08ulU0jcKXD83NJzb4sImX0La90+NVJbzstngpM/WPSy6AQ2vVw6hcSdAiv9MpMbXwOZoqHzxsfQwX4HUa2f8CvaDyw6gU0vi05g0yvUuI0zMjGx7Dapb0B2L5xlRp/8hhe17pW8fCU4+YNFL4tOYNPLpVNI3CnQPza24vZF45fnXkR7HkkypIqwmpePBCd/sOhl0Qlserl0Eh9qcocOHdKTJ0+mHUbFGJmYoK2padntqkrhwb9HX3lpwfrMXe8ks6d657FdzctHgpM/WPSy6AQ2vSrtJCKPqeqhpbaFGncKrNbfT0TI3PV2aGhcsL7wlS+gI0NJhrYuQt9MP7DoBDa9LDqBTa/Qj9s4bY2Nq+4jGzaSfesRyBRdovwM+S/8IzpZnc+HyvHyjeDkDxa9LDqBTS+XTiFxp8D+9vay9pNt28nc/paFKyfGKTz0IJrPJxDZ+ijXyyeCkz9Y9LLoBDa9XDqFxJ0CJ/rL796VufE1yGtuWbBOh3IUHv5S1bU0X4uXLwQnf7DoZdEJbHq5dAqJOwVu27FjTftnDt2JbO9csE6fewZ9+vFKhrVu1urlA8HJHyx6WXQCm14unULiToG1dhuQTIbMW94BGzctWF949ASFF1+oZGjrInTx8AOLTmDTy6IT2PRy6RQSdwrkxsfXfIzUN5A9fA/U1hWtVQpf+hw6dr5ywa2Dq/GqdoKTP1j0sugENr1cOoV+3Cmwnv5+hYFzFI59Bii6bi2byX779yB16c64E/pm+oFFJ7DpZdEJbHqFftzGWU9/v8yuLjK33r5w5dh5Cl/6PKqFdUa2PkLfTD+w6AQ2vSw6gU2v0I/bOB3Nzes6Xg4cRLr3L1ingy9QeOyr6yp3vazXqxoJTv5g0cuiE9j0cukUEncKdLa0rL7TCogImTvfgrRft2C9nuqhkOL83ev1qkaCkz9Y9LLoBDa9XDqFxJ0Cjw4OrrsMydaQeesRaNywYH3h4S+iQ7l1l381VMKr2ghO/mDRy6IT2PRy6RQSdwrc2dm5+k5lIE0byL7tCGSL5u8uFMg/9CB6yX2rzUp5VRPByR8sell0ApteLp1C4k6B08PDFStL2q8jc+dbF66cnKDwhX9EZ2Yqdp5yqKRXtRCc/MGil0UnsOnl0imxxC0iHxGRV0TkyaJ1W0Tk8yLybPy3NanzVzMjk5MVLS/TvY/MzW9YsE5fHaLwlYecDotaaa9qIDj5g0Uvi05g08ulU5I17vuAIyXrfhE4pqo3Asfi5WuOo/v2VbxMecObkF27F6zT559Fn/h6xc+1HEl4pU1w8geLXhadwKaXS6fEEreqHgdeLVn9HuCj8fuPAu9N6vzVTBL9/ebm8N608CZG4RuPUOg/W/HzLUXom+kHFp3AppdFJ7Dp5dKpxtmZIq5T1Zfi9znguuV2FJF7gXsBtu/cyX09PdzV1UVPLsfY1BRH9+3j/t5eultbaayp4dTQEIf37OHhgQGm83kOd3fz4Jkz7GtrA6B3ZIQje/dyrK+P2myWO3bt4tjZsxzYupXJmRn6Rkfnymypr+dgRwfHz53jYEcHIxMT9I+NzW1va2xkf3s7J/r7uW3HDvrHxsiNj89t72huprOlhUcHB7mzs5PTw8OMTE7ObR+6dImeXI6eXK7iTieu6+Lw6KtwZYrJmRnaG5t4+YFP8vRt38b+G25MzKmzpYWpmRlT1yk4+ePU1tTEU0ND9I2OBqcqd+rJ5ciK8Mmnnw5OKzitRKJDnorIbuABVb05Xj6vqpuLto+q6qrPua0NedqTy3GwoyOx8vWlAfKfvx+Kr21zC9nv+F6kviGx8ybtlQbByR8sell0AptelXaqpiFPXxaR7QDx31ccn78q6Mkl289atu8i88Y3L1w5Pkbhi59FC/nEzpu0VxoEJ3+w6GXRCWx6uXRynbg/Bbw/fv9+4B8cn78quKurK/FzyP6bkRtvWrBOcy9SeOQriZ3ThZdrgpM/WPSy6AQ2vVw6Jdkd7K+Ah4H9IjIgIj8G/CbwDhF5Fnh7vHzN4eKXmYiQuf3bkG0LJ3fX009SOH0qkXOGX9F+YNEJbHpZdAKbXi6dEmucpqo/sMymw0md0xfGpqacnEcyWTJvfRf5B/4WLl2cW1/42v9FNrci1+1Y4ei148rLJcHJHyx6WXQCm14uncJ83Cngei5aHR0h/5lPwsz0/MraOjJ33k1m996KnSfMsesHFp3AppdFJ7DpFebjNo7rPozS2kbmzSU3OqavUPjS58if+GLFhkYNfTP9wKIT2PSy6AQ2vcJ83MbpbnU/0mumq5vMwTcuWq/PPkX+gb9BR0fWfY40vJImOPmDRS+LTmDTy6VTSNwp0FjjetybCLnlVjK3vRkyJZf9wij5T/8thWeeXNfY5ml5JUlw8geLXhadwKaXS6eQuFPg1NBQKucVETI33UL2nu+Bls0LN+bzFL52POrrPXX5qspPyytJgpM/WPSy6AQ2vVw6hcSdAof37En1/NK2NRpF7Yb9i7bpC33k7/84+vJLSxy5Mml7JUFw8geLXhadwKaXS6eQuFPg4YGBtENAauvIvvlw1Gitpnbhxkvj5B/8ewrfPIlqoewyq8Gr0gQnf7DoZdEJbHq5dAqJOwWm88kNO7pWMjfsJ3v0fciWrSVblELPIxQ+ez86MV5WWdXkVSmCkz9Y9LLoBDa9XDqFftwpkBsfp6O5Oe0wFqCFPIXHvoY+1bN4Y30DmW99G5nO3SuWUY1e6yU4+YNFL4tOYNOr0k6hH3eV8eCZM2mHsAjJZMnedieZw98ODY0LN05dpvCFz5D/2v9FV/hVWY1e6yU4+YNFL4tOYNPLpVNI3CkwOwdrNZLZ1RXdOu9YPB+sPvME+U9/Ar0wuuSx1ex1tQQnf7DoZdEJbHq5dAqJO7AIadpA5p1Hybz+TSCycOPoMPkH/pbCmWfW1ec7EAgEAldHSNwp0Duy/lHKkkZEyNxyK9kj3wUbNi7cODNN4StfoPB/j6HTV+ZW++C1VoKTP1j0sugENr1cOoXGaSngW8MMvTJF4cQX0XPPLd64sYXsXe9E2rd551UOwckfLHpZdAKbXqFxmnGO9fWlHcKakLp6Mm95J5k77oZsybB+F8fI/+MnKZzq4dhzSyR2z/HtWpWDRSew6WXRCWx6uXSyN2CsB9Rms2mHsGZEBNl3E7Ktg/yXPgfnX53fWChQOHmC19c1op07kS3t6QVaYXy8Vqth0Qlsell0ApteLp1CjTsF7ti1uMW2L8jmLWS//XuR/Tcv2nb9xEXy93+c/LHPoMMvpxBd5fH5Wi2HRSew6WXRCWx6uXQKiTsFjp09m3YI60JqasjefheZu49AXf3c+tHL0eQkOvA8+U9/gvzn70dfWfuY59WE79dqKSw6gU0vi05g08ulU7hVngIHtpYOL+onma5upH0rhePH0FcG2VC7cMxzHewnP9iPXLeTzOtuhY6dSGn3sirHyrUqxqIT2PSy6AQ2vVw6hcSdApMzM2mHUDFkw0YyR96DDpzj0le+SMvUxKJ99OUXyX/uRWRrB3LLrcjO671J4Jau1SwWncCml0UnsOnl0incKk+BvtGlRx7zFREh07mbY/tfR+YdR5Ft25fcT4dyFI59msKnP0HhhbNeDOBi7VqBTSew6WXRCWx6uXQK/bhTYGRigramprTDqDjFXpobjKYFza0w1d3mNjKvuxXp6kakOn9DWrxWFp3AppdFJ7DpVWmn0I+7yri/tzftEBKh2Es6dpB913eSffd3IzuvX/qA8yMUvvQ58v/w1xSeO40Wyp/72xUWr5VFJ7DpZdEJbHq5dArPuFOgpb5+9Z08ZCkv2dZB9u3fgY4MUXj8MfSFJQYpuDBK4cvH4JuPknntrcgN+5BMdfTztHitLDqBTS+LTmDTy6VTSNwpcLCjI+0QEmElL2nbSvatR9DRkSiBP/8cUPKY5uIYhRMPRQn85jcgN34LUjpSm2MsXiuLTmDTy6IT2PRy6RRulafA8XPn0g4hEcrxktY2sm95J9n3/jOke//i2ccALo1T+Npx8p/4SwqnvolOXU4g2vKweK0sOoFNL4tOYNPLpVOocaeAxV+bsDYv2dRK9tsOowdvo/DE19HnnoHSZ9yTlyic/Ao89jCyoxPpvhHp3IOU9BdPEovXyqIT2PSy6AQ2vVw6hcSdAiMTi/s6W+BqvGRjC9k770ZfdyuFJ76BPvs0FPILd9IC+uI59MVzkK1BdnVFSXxnF5Lw+MAWr5VFJ7DpZdEJbHq5dAqJOwX6x8bSDiER1uMlGzaSvf0u9JZbKZzqQU+fgvwSAxrkZ9Bzz0VTjNbWIV3dZPbcCNt3JtKlzOK1sugENr0sOoFNL5dOoR93CljswwiV9dLLk+jpUxTOPgsXyhjYoKEJ2X1DlMS3XlexkdksXiuLTmDTy6IT2PQK/biNY7EPI1TWSxoaybzuENn3/DOyR9+HHHg9bFhhkvrLE+gzT5D/x0+S/+THKDz2VXR0ZN2js1m8VhadwKaXRSew6RX6cRunrbEx7RASIQkvEYEt7WS3tKO33g5DOQp9z0a3yi9PLn3Q+EUKT34dnvw6bGols+fG6Jn4xk1rPr/Fa2XRCWx6WXQCm14unULiToH97e1ph5AISXuJCGzbTnbbdvSNb0ZfGkDPPou+cBamryx90IVRCj2PQM8jSNu2KIHv3os0bSjrnBavlUUnsOll0Qlserl0SiVxi8jzwEUgD8wsdx/fKif6+9nX1pZ2GBXHpZdkMtFQqjuvR/Mz6MALURIfeB7y+SWP0ZFX0JFX4NETSNtWZPsuZPtO2LYdqVn6n4LFa2XRCWx6WXQCm14undKscb9VVYdTPH9q3LZjR9ohJEJaXpKtQbq6oasbnb6CvnA2SuKD/bDkM26dT+JPfh0yGWTrdmT7TmT7LmjfhmSi5h8Wr5VFJ7DpZdEJbHq5dAq3ylOgf2yMA9u2pR1GxakGL6mtQ27YDzfsj1qmn3sO7TuDvjK4/EGFAvryi+jLL0LPI1E3s+t2IB07GbqSR7du9Wb+8HKohuuUBBa9LDqBTS+XTmklbgU+JyIK/Imq/mlKcaRCbnw87RASodq8pKER2X8z7L8ZvTSOPn8G7XsWfXVo5QOnr6ADz6MDz3PD+Dj53m8iHTvnb603t3idyKvtOlUKi14WncCml0unVPpxi8hOVX1RRLYBnwd+SlWPl+xzL3AvwPadO2/9jQce4K6uLnpyOcampji6bx/39/bS3dpKY00Np4aGOLxnDw8PDDCdz3O4u5sHz5yZe+bQOzLCkb17OdbXR202yx27dnHs7FkObN3K5MwMfaOjc2W21NdzsKOD4+fOcbCjg5GJCfrHxua2tzU2sr+9nRP9/dy2Ywf9Y2Pkxsfntnc0N9PZ0sKjg4Pc2dnJ6eFhRiYn57Zvqq9nT2srPbmcGafOlhayIjx/4UL1O13fySOPPULX5DjXXbzAxPlX2VRfz8T0NNOFAu2NTQxPTtCQrUGAyfwMm+sbGL9yhRkt0Na+lSczddTt6qKwrYPHL46n72TxOq3B0tTP4QAAFD9JREFUqa2piS+/8AL33HhjcKpyp55cjtdu28a5CxeC0wpO12/evGw/7tQHYBGRDwHjqvrby+1jbQCW+3p6+MDBg2mHUXF89dKLY1EL9ZdeRHMDC7qZvTQ+zvbmFfqPA2xqjWrj1+1A2rfBhuaqrpH7ep1Ww6KXRSew6VVpp5UGYHF+q1xENgAZVb0Yv38n8Kuu40iTjtUSgaf46iUbW5CNN8G+m6IBW86/GiXxlwao6Xt29QIujKIXRtFnnoiW6xuirmftW2H2b+OGqknmvl6n1bDoZdEJbHq5dErjGfd1wN/F/4nVAP9bVR9MIY7U6GxpSTuERLDgJSLQ2oa0tsFNt/DqgRwdWeYSub6SWzwJSilTl9HBF9DBF+bXNTZFybxta9RqvW0r0pjOkI8WrtNSWPSy6AQ2vVw6OU/cqtoHvM71eauJRwcHzbWoBJtej+ZyHDh4ENnaAbfcis7MwFAuSuKDA+jIEFFby1WYnJhr8DZH04bFybwh+dGXLF4nsOll0Qlserl0Ct3BUuDOzs60Q0gEi16lTlJTA9t3Rf293wB6ZQp9eRDNDcLwUNRifWa6vMInLqETZ9H+s/PrNmyMEnicyGnbitQ3VNDI5nUCm14WncCml0unkLhT4PTwsLlRg8Cm12pOUlePdO6Bzj0AqBbgwnl0ZAgdfgVGhtBXh5eeonQpLl1EL11EX+ibX9fYhGxqjRrBbd6CbNoMm7dAQ+NVPTe3eJ3AppdFJ7Dp5dIpJO4UGJlcZnIMz7HotVYnkQxs3oJs3gI37AdAC4WoAdvwK9Gt9eFX0NFhKBTKK3RyAp2cgNyLC2/K19VHCX3zlvhvK7K5FZpWbtVu8TqBTS+LTmDTy6VT6t3BysFadzCLc9GCTa+knLSQh9FX52rmOjIEoyOgZSbzlaipLaqhxwl90xbYuBGRjMnrBOH75xMWvVzOxx1q3Clwf2+vuT6MYNMrKSfJZKPn121bYd9NAGh+BkZH0OGh+VvtY6Pl18xnmZmOxmEfeWVhDT2ThZbNnL08Reu+/dEIcBs3wcaNUUM5yVTMLw3C988fLHq5dAqJOwUsdoUAm14unSRbA+3XIe3Xza3TQgHGx9Dzr8L5uL/4+VG4MFr+c/NZCnk4P0LX5GUKj5cMz5jJRIm8uQU2tiAbW+LEHi/X1lXAMFnC988fLHqZ7g4WwNwtolkseqXtJJkMtGxGWjbD9fPrVTVqyBYncT3/atQo7vyry89NHlObXaJmXSjA2Hl07HxUfun2+oY4qW9CNm6M/jZHf9lQHbX1tK9VElh0ApteLp1C4k6BnlyOgx0daYdRcSx6VauTiMzXkHd1za1XVZicmEvmeuF8NBLchdG5oVwvXrlCc90aa9BTl9Gpy4tvv0NUW9+wEdmwMRrudcMG2BDfft/QDE0booZ0CY8cV63Xaj1YdAKbXi6dQuJOgbu6ulbfyUMsevnmJCJRomzaEPU1L0KnLsP5UbIv9pOhgF68CBcvoONjMHX56k9aKETlXLwQnWepfbI1UVJvao5q6E3NsKE5Xrcher/O5O7btSoHi05g08ulU0jcKdCTy9Hd2pp2GBXHopclJ6lvgOu289VXz/Pdr3nNgm06fQUujqEXx6Jn6hfH4qR+EcbH1t5ArpT8zMq34qEouW+IurRtiP7S2BQNDxu/pGbp/7YsXatZLDqBTS+XTiFxp8DY1FTaISSCRa9rxUlq62BLO7KlfdE21QJMXJpP7BfHolr67N/LFeq/Wk5yB6itixJ4QxM0Ns4l9uaXchQ2FiX5hsaojYDHWPz+gU0vl06hH3cKWOzDCDa9gtPq6PQVGL+IXhqHS+MwMR6/v4ROxOvW2gL+KpjOF0oa3gk0NEBDSY29sSle3xjdhWhojN4vU5NPE4vfP7DpFfpxG8diH0aw6RWcVkdq6+ZnVFsCVYUrU9HY7JfGo9r7pfG5pK7xctljvC/D8OREydzpGt0NuDyJnh8pXrs02Zooodc3InFip74hmvhl7v18onfR4M7i9w9seoV+3Max9mxnFotewWn9iAjUN0SJb6XkPn1lYSK/NI5OXoLJyflhXycnlh1drnG9Neb8zPx5i2Nb9oC4Rl/fgNQ3Ru/r6uPl6C91dVGtPl5PfX00sl2ZCd/i9w9serl0Cok7Bdb9H0yVYtErOLlBRKLkVle/bHKHotp7nMS16O+ll3O01tbMJ/j1tJQvi6IaPaOlW5ZHMlECr6svSupFy3NJv57NFy9EXflq66L9stnEa/kuqMbv4Hpx6WTv0/OAU0ND3LZzZ9phVByLXsGpuiiuvbN5C8Up7KGengW3KrWQjxLr5OR8Mp9N9pcvRwl3ajJ6PzW5/pbz5aKFooRfsqlkuXV8nHzx7f9MZv4HTm0d1NXNJ/W6eqQuXldXH9X26+qLtkf7VkPi9/k7uBwunULiToHDe/akHUIiWPQKTv5Q6iWZbNSdrKmZ1VKVqkbP2GcT+uXJqMZ+eRKdTeylyX6VEeoqQWtDyVzsheWTPqxS0wdAoLZ2PonX1MaJvS5aX1sX/SCI30f71iNzy0X7ZbNX7WXxO+jSKSTuFHh4YIDOTZvSDqPiWPQKTv6wHi8RmU9MG1tWTfQwW6O/PJ/gpy7D1BRcmUKnpqL1V6biUeem5t6vpYX92NQUDRW9BRu3JYh/dFxd8o/JZBYk87mEX1M791dKlmf3PXXuHLu+5TXz22prIVtTFXcDrhaX/65C4k6B6Xw+7RASwaJXcPIH115RjT4apQ4oK9lDPAtcUSLXK1PR8lTp8mXGXhpkW8vGaN8rU+5u55dDoRD9EInbEqzlR8Brxy+Rf/aJkrUCNTXztfqa2uiOQE1N/Jr9AVATdd2b2zb/d8H+sz8Wamogk3zbAJffv5C4U+Bwd3faISSCRa/g5A++eEm2Bppqykr4m8fHqYmfcasq5PNRbTlO5HrlCly5ssS6qWj9lamon/3sPg5u75dDa2PDEmvjxxUz0zA5t2ZJ1j76iMwlfWpqIFsbJ/+aqBtg8Y+DbDbali3+UTC/nyxzzOE9u9cc1dUSEncKPHjmjLk+jGDTKzj5g0WvYicRmU8WjdFAH2utQ2qhECfw6blErleuRMnyyuz6KXR2+5X5fbXoGKavwDoG7xqZnCzpc580uuiHy0rRX43Z9Pg4My0tZL/7n0eT6yRISNwpsK9t+e4uPmPRKzj5g0WvSjtJJjPfKn923VWUs6D2P/dDYBqdmZ5fnpme+6vT0wvWyfnRqNvbTNF+DkbXS5Kmmtro8cE6Gu2VS0jcgUAgEFgTS9X+ofwfAS/097Ors3PBOi0UFibymZmopj8zs3Bdfib+kRCvn/sbvdeSZWam3bYNqKlN/hSJnyGwiN6REe4s+dJawKJXcPIHi14WnWBpLynqo75gfQXOp4V8UYKfmX/lZ6K7BDP5aFt+4V8t2m/hMbPrZvedYWJmnE3Uhxq3VY7s3Zt2CIlg0Ss4+YNFL4tO4N5LMlmoyy76UQCV+WEAUHvxItmmRidd2vye885TjvX1pR1CIlj0Ck7+YNHLohPY9Dp29mzUY8ABIXGnQK2DWylpYNErOPmDRS+LTmDTy6VTSNwpcMeuXWmHkAgWvYKTP1j0sugENr1cOoXEnQLHzp5NO4REsOgVnPzBopdFJ7Dp5dIpJO4UOLB1a9ohJIJFr+DkDxa9LDqBTS+XTiFxp8DkjN8DDSyHRa/g5A8WvSw6gU0vl04hcadA3+ho2iEkgkWv4OQPFr0sOoFNL5dOousYb9YVhw4d0pMnT6YdRsUYmZigralp9R09w6JXcPIHi14WncCmV6WdROQxVT201LZQ406B+3t70w4hESx6BSd/sOhl0Qlserl0Cok7BVrqF4/eYwGLXsHJHyx6WXQCm14unULiToGDHR1ph5AIFr2Ckz9Y9LLoBDa9XDqFxJ0Cx8+dSzuERLDoFZz8waKXRSew6eXSyYvGaSIyBFi60u3AcNpBJIBFr+DkDxa9LDqBTa9KO3Wp6pKdw71I3NYQkZPLtRb0GYtewckfLHpZdAKbXi6dwq3yQCAQCAQ8IiTuQCAQCAQ8IiTudPjTtANICItewckfLHpZdAKbXs6cwjPuQCAQCAQ8ItS4A4FAIBDwiJC4K4yIHBGR0yJyRkR+cYntPysiT4nI4yJyTES6irblRaQnfn3KbeTLU4bTB0RkqCj2Hy/a9n4ReTZ+vd9t5CtThtfvFTn1isj5om3Veq0+IiKviMiTy2wXEfmD2PlxEXlD0baqvFZlOP1Q7PKEiJwQkdcVbXs+Xt8jIlUz4UEZTneLyIWi79gvF21b8XubJmV4/XyR05Pxv6Mt8bZqvVadIvJQ/P/2KRH56SX2cfvvSlXDq0IvIAs8B3QDdfD/t3e3MVYddRzHvz9TlBaIhW5VbLHYB9IHrS00DQIimEZrH1hMMQG1CcZEIbHGaHijRiO+kKTGVNMoJqbBGmW1VCpiSyHBii1s6ZbaUkUprARLSNpATUEqZuXni5kbDte7u/cuy95z1v8nIXvOnDn3zn/nzs49cw4zPA9cW5dnPnBB3l4O/KJw7Hi7YxhiTEuB+xucOwnozT8n5u2J7Y6p2bjq8t8DPFDmusrlmgtMB17s5/htwGOAgJnA0xWoq8FimlUrK/DRWkx5/wDQ0e4YhhDTPGBjg/SWPrdli6su753A1grU1WRget6eAOxt8DdwRNtVXHEPr5uBfbZ7bf8b6AI6ixls/872ibzbDVw6wmVs1aAxDeAjwBbbR22/BmwBbj1H5WxVq3EtAdaOSMnOgu1twNEBsnQCDzrpBi6UNJkS19VgMdnenssM1WhTzdRTf86mPZ5zLcZVlTZ12PauvH0M2ANcUpdtRNtVdNzD6xLg74X9l/nfCi76DOlbWs1YST2SuiUtPBcFHIJmY7orDxGtkzSlxXPboemy5dsZ7wa2FpLLWFfN6C/uMtdVK+rblIHNkp6V9Nk2lWmo3i/peUmPSboup42KepJ0AakDe7iQXPq6kjQVuBF4uu7QiLar8872BcLQSPoUcBPwwULyZbYPSboc2Cppt+397SlhS34DrLV9UtLngJ8AH2pzmYbTYmCd7f8U0qpaV6OWpPmkjntOIXlOrqe3AVsk/SVfFZbdLtJn7Lik24BHgKvaXKbhdCfwlO3i1Xmp60rSeNIXjS/afr2dZYkr7uF1CJhS2L80p51B0i3AV4EFtk/W0m0fyj97gSdI3+zabdCYbB8pxPFjYEaz57ZRK2VbTN2QXknrqhn9xV3muhqUpOtJn71O20dq6YV6egVYTxpqLj3br9s+nrcfBcZI6qDi9VQwUJsqXV1JGkPqtH9m+1cNsoxsu2r3jf/R9I80gtFLGlatPThyXV2eG0kPl1xVlz4ReEve7gBeogQPnTQZ0+TC9seA7rw9Cfhbjm1i3p7U7piajSvnu5r00IzKXleF8k2l/4eebufMh2h2lr2umojpXcA+YFZd+jhgQmF7O3Bru2NpMqZ31D5zpA7sYK6zpj63ZY0rH38r6T74uCrUVf69PwjcN0CeEW1XMVQ+jGz3Sfo88Djp6c8HbP9J0kqgx/YG4F5gPPCQJICDthcA1wA/knSKNBKyyvaf2xJIQZMxfUHSAqCP1CCX5nOPSvoW8Ex+uZU+c2isbZqMC9KVQZdzK8xKWVcAktaSnkjukPQy8A1gDIDt1cCjpCdg9wEngE/nY6WtqyZi+jpwEfCD3Kb6nBZ7eDuwPqedB/zc9qYRD6CBJmJaBCyX1Ae8ASzOn8GGn9s2hNBQE3FB+nK/2fY/C6eWtq6A2cDdwG5Jf8xpXyF9YWxLu4qZ00IIIYQKiXvcIYQQQoVExx1CCCFUSHTcIYQQQoVExx1CCCFUSHTcIYQQQoVExx1CCUlaKMmSrh7G15wnaWPeXlBbWSq/17VDeL0nJN3U4jn3SZrbQv6LJZXlvwWFUArRcYdQTkuAJ/PPYWd7g+1VeXch0HLH3SpJFwEz3cI0lrZfBQ5Lmn3uShZCtUTHHULJ5DmR55Dm3V5cSJ8n6feSfi2pV9IqpbWod+Z1jK/I+dZIWp0XQdkr6Y4G77FU0v2SZgELgHvzOshXFK+kJXVIOpC3z5fUJWmPpPXA+YXX+7CkHZJ2SXoox1DvLmBT4ZwDkr6d37dH0nRJj0vaL2lZ4bxHgE8O+RcawigTHXcI5dMJbLK9FzgiaUbh2PuAZaTZ2+4Gptm+mTRP9z2FfFNJU2XeDqyWNLbRG9neDmwAVti+wQMvlLIcOGH7GtKMWDMgde7A14BbbE8HeoAvNTh/NvBsXdpB2zcAfwDWkGYMmwl8s5CnB/jAAOUK4f9KTHkaQvksAb6Xt7vyfq3De8b2YQBJ+4HNOX03ML/wGr+0fQp4SVIvac71szUX+D6A7RckvZDTZ5KG2p/KU1a+GdjR4PzJwKt1abWpZXcD453WOz4m6aSkC23/A3gFeOcwlD+EUSE67hBKRNIk0pKo75Vk0nzUlrQiZzlZyH6qsH+KM9tz/VzGrcxt3Mfp0biGV+r1xQa22B7sfvwbDV6vWP762GrxjM3nhhCIofIQymYR8FPbl9meansKaUWhVoeKPy7pTfm+9+XAXwfIewyYUNg/wOmlWRcV0rcBnwCQ9B7g+pzeDcyWdGU+Nk7StAbvswe4ssU4AKYBLw7hvBBGpei4QyiXJaS1iIsepvWnyw8CO0lLDS6z/a8B8nYBKyQ9lzv675BWpnqOtGxpzQ+B8ZL2ACvJw/f5ye+lwNo8fL6DxkPzvyWtHNWq+fncEAKxOlgIo46kNcBG2+vaXZZ6kp4E7sj3rps9ZxvQafu1c1eyEKojrrhDCCPpy+R1jJsh6WLgu9Fph3BaXHGHEEIIFRJX3CGEEEKFRMcdQgghVEh03CGEEEKFRMcdQgghVEh03CGEEEKFRMcdQgghVMh/AfWqeif7vy3KAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "From the graph, we can see that as the amplitude increases from 0 to 2, the period exponentially decreases. "
      ],
      "metadata": {
        "id": "oHm-mgeivMgN"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BY5x-mZsYoqG"
      },
      "source": [
        "### PROBLEM 5 (20 pts.)\n",
        "\n",
        "**Differentiation using central difference**\n",
        "\n",
        "1. Calculate the derivative of the function $1 + \\frac{1}{2} \\tanh 2x$ in the range $-2 \\le x \\le 2$ using central difference. Describe each step of the calculation clearly and concisely. *(10 pts.; 5 pts. for the code & 5 pts. for the discussion text)*\n",
        "\n",
        "2. Calculate an analytic formula for the derivative and make a graph with your numerical result and the analytic answer on the same plot. Plot the exact answer as lines and the numerical one as dots. Describe how the two results compare. *(10 pts.; 5 pts. for the code & 5 pts. for the discussion text)*"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "The **Central Difference** is given by the equation \n",
        "\\begin{align}\n",
        "  \\frac{df}{dx} \\simeq \\frac{f(x+h/2)- f(x-h/2)}{h}\n",
        "  \\end{align}"
      ],
      "metadata": {
        "id": "CVOcpoMB44U_"
      }
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "ZWXWJyYdLS5G"
      },
      "source": [
        "\n",
        "h = 0.000000000001 # approximation that approaches 0\n",
        "def func(x):\n",
        "  '''the function given in the problem'''\n",
        "  return 1 + 0.5 * tanh(2*x) \n",
        "\n",
        "def CentralDif(f,x,h):\n",
        "  '''the function for calculating the derivative using the central difference'''\n",
        "  derivative = (f(x+h) - f(x-h))/ (2 * h) #function for the central difference as shown above\n",
        "  return derivative\n",
        "  \n",
        "calculated_answers = [] #stores the calculated derivative using central difference in the range -2 to 2 with step 1\n",
        "for i in range(-2,3,1): #appends the derivative of the function in the list \n",
        "  calculated_answers.append(CentralDif(func,i,h))"
      ],
      "execution_count": 996,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "def deriv_m(x): \n",
        "  ''' produces the actual value of the derivative'''\n",
        "  return ((1/ cosh(2*x)) **  2 )\n",
        "\n",
        "exact_answers = [] # stores the exact derivative in the range -2 to 2 with step 1\n",
        "for i in range(-2,3,1):\n",
        "  exact_answers.append(deriv_m(i))"
      ],
      "metadata": {
        "id": "UVTVOiu6Xsmv"
      },
      "execution_count": 997,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "errors = [] # stores the percent error between the actual and calculated derivatives\n",
        "for i in range(0,5):\n",
        "  error = abs((exact_answers[i] - calculated_answers[i]) / exact_answers[i]) * 100\n",
        "  errors.append(error) "
      ],
      "metadata": {
        "id": "N5XZpLAe7kJz"
      },
      "execution_count": 998,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "df = pd.DataFrame(\n",
        "    {'Actual value': exact_answers,\n",
        "     'Calculated value' : calculated_answers,\n",
        "     'Error': errors\n",
        "    })\n",
        "df.index += 1\n",
        "df #displays the integral given the number of slices and its error"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 204
        },
        "id": "Af21FhGF7K-c",
        "outputId": "f79e5bae-27a1-4093-9937-a7152031808d"
      },
      "execution_count": 999,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "   Actual value  Calculated value     Error\n",
              "1      0.001341          0.001332  0.647530\n",
              "2      0.070651          0.070610  0.057523\n",
              "3      1.000000          1.000033  0.003339\n",
              "4      0.070651          0.070610  0.057523\n",
              "5      0.001341          0.001332  0.647530"
            ],
            "text/html": [
              "\n",
              "  <div id=\"df-75cc6637-cc30-40a5-bc0f-3296fbb90eb9\">\n",
              "    <div class=\"colab-df-container\">\n",
              "      <div>\n",
              "<style scoped>\n",
              "    .dataframe tbody tr th:only-of-type {\n",
              "        vertical-align: middle;\n",
              "    }\n",
              "\n",
              "    .dataframe tbody tr th {\n",
              "        vertical-align: top;\n",
              "    }\n",
              "\n",
              "    .dataframe thead th {\n",
              "        text-align: right;\n",
              "    }\n",
              "</style>\n",
              "<table border=\"1\" class=\"dataframe\">\n",
              "  <thead>\n",
              "    <tr style=\"text-align: right;\">\n",
              "      <th></th>\n",
              "      <th>Actual value</th>\n",
              "      <th>Calculated value</th>\n",
              "      <th>Error</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>1</th>\n",
              "      <td>0.001341</td>\n",
              "      <td>0.001332</td>\n",
              "      <td>0.647530</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>2</th>\n",
              "      <td>0.070651</td>\n",
              "      <td>0.070610</td>\n",
              "      <td>0.057523</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>3</th>\n",
              "      <td>1.000000</td>\n",
              "      <td>1.000033</td>\n",
              "      <td>0.003339</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>4</th>\n",
              "      <td>0.070651</td>\n",
              "      <td>0.070610</td>\n",
              "      <td>0.057523</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>5</th>\n",
              "      <td>0.001341</td>\n",
              "      <td>0.001332</td>\n",
              "      <td>0.647530</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>\n",
              "      <button class=\"colab-df-convert\" onclick=\"convertToInteractive('df-75cc6637-cc30-40a5-bc0f-3296fbb90eb9')\"\n",
              "              title=\"Convert this dataframe to an interactive table.\"\n",
              "              style=\"display:none;\">\n",
              "        \n",
              "  <svg xmlns=\"http://www.w3.org/2000/svg\" height=\"24px\"viewBox=\"0 0 24 24\"\n",
              "       width=\"24px\">\n",
              "    <path d=\"M0 0h24v24H0V0z\" fill=\"none\"/>\n",
              "    <path d=\"M18.56 5.44l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94zm-11 1L8.5 8.5l.94-2.06 2.06-.94-2.06-.94L8.5 2.5l-.94 2.06-2.06.94zm10 10l.94 2.06.94-2.06 2.06-.94-2.06-.94-.94-2.06-.94 2.06-2.06.94z\"/><path d=\"M17.41 7.96l-1.37-1.37c-.4-.4-.92-.59-1.43-.59-.52 0-1.04.2-1.43.59L10.3 9.45l-7.72 7.72c-.78.78-.78 2.05 0 2.83L4 21.41c.39.39.9.59 1.41.59.51 0 1.02-.2 1.41-.59l7.78-7.78 2.81-2.81c.8-.78.8-2.07 0-2.86zM5.41 20L4 18.59l7.72-7.72 1.47 1.35L5.41 20z\"/>\n",
              "  </svg>\n",
              "      </button>\n",
              "      \n",
              "  <style>\n",
              "    .colab-df-container {\n",
              "      display:flex;\n",
              "      flex-wrap:wrap;\n",
              "      gap: 12px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert {\n",
              "      background-color: #E8F0FE;\n",
              "      border: none;\n",
              "      border-radius: 50%;\n",
              "      cursor: pointer;\n",
              "      display: none;\n",
              "      fill: #1967D2;\n",
              "      height: 32px;\n",
              "      padding: 0 0 0 0;\n",
              "      width: 32px;\n",
              "    }\n",
              "\n",
              "    .colab-df-convert:hover {\n",
              "      background-color: #E2EBFA;\n",
              "      box-shadow: 0px 1px 2px rgba(60, 64, 67, 0.3), 0px 1px 3px 1px rgba(60, 64, 67, 0.15);\n",
              "      fill: #174EA6;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert {\n",
              "      background-color: #3B4455;\n",
              "      fill: #D2E3FC;\n",
              "    }\n",
              "\n",
              "    [theme=dark] .colab-df-convert:hover {\n",
              "      background-color: #434B5C;\n",
              "      box-shadow: 0px 1px 3px 1px rgba(0, 0, 0, 0.15);\n",
              "      filter: drop-shadow(0px 1px 2px rgba(0, 0, 0, 0.3));\n",
              "      fill: #FFFFFF;\n",
              "    }\n",
              "  </style>\n",
              "\n",
              "      <script>\n",
              "        const buttonEl =\n",
              "          document.querySelector('#df-75cc6637-cc30-40a5-bc0f-3296fbb90eb9 button.colab-df-convert');\n",
              "        buttonEl.style.display =\n",
              "          google.colab.kernel.accessAllowed ? 'block' : 'none';\n",
              "\n",
              "        async function convertToInteractive(key) {\n",
              "          const element = document.querySelector('#df-75cc6637-cc30-40a5-bc0f-3296fbb90eb9');\n",
              "          const dataTable =\n",
              "            await google.colab.kernel.invokeFunction('convertToInteractive',\n",
              "                                                     [key], {});\n",
              "          if (!dataTable) return;\n",
              "\n",
              "          const docLinkHtml = 'Like what you see? Visit the ' +\n",
              "            '<a target=\"_blank\" href=https://colab.research.google.com/notebooks/data_table.ipynb>data table notebook</a>'\n",
              "            + ' to learn more about interactive tables.';\n",
              "          element.innerHTML = '';\n",
              "          dataTable['output_type'] = 'display_data';\n",
              "          await google.colab.output.renderOutput(dataTable, element);\n",
              "          const docLink = document.createElement('div');\n",
              "          docLink.innerHTML = docLinkHtml;\n",
              "          element.appendChild(docLink);\n",
              "        }\n",
              "      </script>\n",
              "    </div>\n",
              "  </div>\n",
              "  "
            ]
          },
          "metadata": {},
          "execution_count": 999
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "x = [-2, -1, 0, 1, 2]\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "ax.plot(x,exact_answers, color = 'salmon', marker = 'o' ,linewidth = 9, alpha = 0.5, label = 'Actual value')\n",
        "ax.plot(x,calculated_answers, '--', color = 'slategray', marker ='v', alpha = 1, label = 'Calculated value')\n",
        "ax.legend()\n",
        "ax.set_title(\"Actual vs Calculated Value of Derivatives using Central Difference\")\n",
        "ax.set_xlabel('range')\n",
        "ax.set_ylabel('derivatives')"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 312
        },
        "id": "0DhjKPbzXxha",
        "outputId": "1c1ef999-3dad-438b-d5cc-f0a6ee8b5155"
      },
      "execution_count": 1000,
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "Text(0, 0.5, 'derivatives')"
            ]
          },
          "metadata": {},
          "execution_count": 1000
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEWCAYAAADPZygPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydeZwdVZX4v+ctvaa39L5kgRCWBBJCQiABFRQFFWFQGeKOGzrjNoOMo/4UGUZHR53RUVzAGUVHCSAKBMURGUEYkwCJxAiEJUCW3pf0vvd75/dHVd6rV/26+6XTr9fz/Xz6069u3bp1qurWPXXvPfccUVUMwzAMYyYIzLQAhmEYxsLFlJBhGIYxY5gSMgzDMGYMU0KGYRjGjGFKyDAMw5gxTAkZhmEYM8acVkIicoOI/HSm5RgPEVEROek4y5i26xSRC0SkNg3lLnfvRWiqy57gvNkicp+IdIrIz6fxvE+LyAVpKHepiPSISHCqy55qXDlPnGk5ZgIROSAiFx3H8e8QkQc82+eJyAvuPf0rESkXkUdEpFtE/m1qpJ4ZjksJicjDItIuIpkp5r9aRP7veM45E4jI20Vkl1sBGkTkNyJy/kzLlQwRuVVEvpimsp8VkfclSf+EiOxKxzmngLcC5UCxql7p3+kq+GH3Ze4WkedF5CYRqTyek6rqalV9+HjKcOVLaMxU9ZCqLlLVyPGWnW5cOV9KR9kicrKI/FxEWt0PjL0icu3xKufp+Fhy39EhT517SkS+LCIFR/Oo6s9U9XWew24EbnLv6T3ANUArkK+qn0yXrNPBpJWQiCwHXgEocNkUyTPrEJFrgW8C/4LTmC0FvgtcPpNyzRA/Bt6dJP1d7r7ZyDLgeVUdGSfPHaqaBywGrgAqgN2TUUTT3dNbiIjICuAx4DBwhqoWAFcCG4C8aTj/VDzjr7p1rhR4L3Au8EcRyR0j/zLgad/2MzoJbwOzro6q6qT+gOuBPwL/DvzKt28J8EugBWgDbgJOAwaACNADdLh5HwY+4Dn2auD/PNv/gVPZuoDdwCs8+24AfjqGfPuASz3bIVees4As4KeubB3AE0B5kjIKXFmvHOc+bAR2uOU0uNea4dmvwEnu72zg34CDQCfwf27aBUCtr9wDwEXJrhP4OdDolvEIsNpNvwYYBoZcue9z06uAX7jX/zLwcU9Z2cCtQDvwDPAPflk8eWuAEWCZJ22Ve74S4I3Ak+6zOgzc4Mm33L0XIf/1jXGN5wLb3fv6Z+CCcZ7BaW496sB5US9z0//JlW3YvR/vT3LsqDoEBN1zft2Tdimwxz3HdmCN71n9I7AXGMSpaweAi9x73w8s9uRfh/MVGwZWAL/HqYutwM+AQjfffwNR9/ge4FPe+whcBezyyf73wDb3dybwdeAQ0AR8H8h295UAv3Kv5wjwKBBIcn8Snpv/nQVOAv6AUxdbcRR6srp/K/Ad4NdAN44SWeHJ+zrgObec77plfsAvj5v3p8CvJ2ifxqw/rvz/jNN+dQMPACXuvkOu3D3u3yacNumPwDfc5/TF8Z5bsvrtk+1W4Iu+tDyc9uOj/nYQeNFXD7aS+J5fhNOh+LSbtw24E7fOeZ7h+93re8RNfx9OO9kO/JbE91qBDwMvuPfwO4B49n/QPbYbp904a6K2ZsxnNVGGcR7yfuBvgfXuDSn3vcDfAHJxGvzz/Tc2WYVOlgd4J1CM89J9EqfxzRqrAfEcdz3wM8/2G4F97u8PAfcBOa6863G6tf4yLsFpdEPj3If1OBU+5D7sfcDfjfEifse93mr3vJtxGooLODYl9D6cSpuJ00vbM1YFx6mcu937kQGcCLwEXOzu/wpOA7QY5+PhKb8sPrl+B3zOs/1l4B739wXAGe451+A0fH+VrDFjHCXk3p824A1uWa91t0uTyBPGqYufda/v1TgvxikT1ZHx9uMMfzzm/l4HNAPnuM/tPa78mZ5r2ePev2z/9eE0Vh/0lP014Pvu75Pc68vE+Sp+BPjmWI0ZiUoox73WlZ79TwBb3N/fALa5zzYPp85/2fPcvu/evzDOqIYkuQ8Jz83/zuI0iP/PfU6xdz1J3b/VfYYbXdl/Btzu7ivB+XB5s7vvEzhtylhKqBF47zjPdNz648r/InAyzkfYw8BXxrneq3HagY+58mUf63PzyXcrPiXkpv8EV4kzuh3014OEMtx7thPnQzETuBnY6rumn+C0ydk4Izn7cT7gQsDngO2+Z/croBBn9KcFuMTddyVQB5wNiHsvljFBWzPW36SG49z5kGXAnaq6G+eBvt3dvRFHG/6Dqvaq6oCqTnoeSFV/qqptqjqiqv+Gc4NPSeHQ24DLRCTH3X47zgsDTgUvxnlBIqq6W1W7kpRRDLTqOEM57rE7XfkO4Dz8V/nziUgAR3l8QlXr3PNuV9XBFK7Ff84fqmq3e+wNwFrveLKPs3FevhtVdUidMfofAFvc/X8NfElVj6jqYeBbE5z+xzjDb0ev6R1uGqr6sKr+RVWjqroX536Puhcp8E7gflW93y3rd8AunEbFz7nAIpxGZEhVf4/z8rxtEuf1Uo/TeIPTw7xZVR9zn9uPcXo853ryf0tVD6tqf5Kybjsqj4gIzr2/DUBV96vq71R1UFVbcEYWUrpnqtoH3OspeyVwKrDNPc81wN+7z7YbZ0j56HMfBipxvn6HVfVRdVuYY2QYpy2oSuFdv1tVH3ffp58BZ7rpbwCeVtVfuvu+haNoxqIYp9cwFqnUnx+p6vPu87rTI8tY1Kvqt933vP94ntt45yBe546VDwP/T1VrPe3CW31Dbze4bXK/m//LqrrPvef/ApwpIss8+b+iqh2qegh4iPg9+gDOcOIT6rBfVQ8ycVuTlMnOCb0HeEBVW93t29w0cL4GD47XcB8LInKdiOxzJx87cIbISiY6TlX34/RK3uQqostcOcEZ5vgtcLuI1IvIV0UknKSYNqBkvDFUd4L0VyLSKCJdOA8zmXwlOF+KL04k+3iISFBEviIiL7rnO+ApPxnLgCoR6Tj6h9NrKHf3V+EMnR3l4AQi/BKoFJFzcXo+OThDLIjIOSLykIi0iEgnTkWf8FmNIfOVPpnPx2k0/VQBh1U16ruG6kmc10s1zjDVUXk+6ZNniXvuoxz2F+DhF8Amd47plThDK48CuFZOt4tInfs8f8qx3bOYgsP50LrHVU6lOM9mt0fm/3HTwemN7QceEJGXROTTx3BOL5/C+Rp+3LUIHGW44sGrWPpwPh7AVwddZTiehWYbyevCUVKpP2PJMhYJz3cKnlsyvHXuWFkG3O253n04Ux/lnjyHffn/w5P/CM5z9L43Y92jJSRvxyZqa5JyzEpIRLJxvp5f5Ta8jTjj0GtFZC3OhS4do+FO9qXVi/OyHKXCc65X4FTyvwaKVLUQZ8xYUhR3K84LejnOJN5+APfL759UdRXOkNilJJ9w34HzxftX45zje8CzOEMi+Tg3PZl8rThzYiuS7Eu4B66FT2mSfOA0NJfjjAMX4HS18ZzTf48PAy+raqHnL09Vj34VNuBUqqMsHeO8TuFOA3cXzv16F86QypC7+zac4Z8l6kwWf5+xn9WYz92V+b99Mueq6leSlFMPLHF7Zd5rqBvvOsbDLetNuIrCledLPnlyVHWr57AxexGq2o4z73AVzvO73dPr+Bf32DPc+vNOEu/ZRL2T3wGlInImTl0/+qHVijOHsNojc4GqLnJl6lbVT6rqiTgfaNeKyGuSlN/r/k/6rFS1UVU/qKpVOMPc35VjX5LQgDOMBMR6izVjZ+dB4C3j7D+W+uNnrPvtT5/ouR0TIrII551+dKK8Y3AYeL3vmrNU1fseqC//h3z5s1V1e4rnStaOTdTWJGUyPaG/wtGwq3C6Z2fijCs+itMwPY5Tqb4iIrkikiUi57nHNgE1IpLhKW8P8GYRyXEr7/s9+/JwxmJbgJCIXA/kH4Ost+NMeP4N8ZcTEblQRM5wG/sunCGFqP9gVe3EGd/8jji2+TkiEhaR14vIVz0ydgE9InKqe65RuF/qPwT+XUSq3B7NJnHM258HskTkjW6P7HM4w47JyMNRjG04DcO/+PY34YzFHuVxoFtE/lGcNTNBETldRM52998JfEZEikSkBmfceyJ+jNOgvoVEq7g84IiqDojIRuJDtMnYA2xx7+cGHFPqo/wUpwd7sStvljjrl5I1TI/hfKV9yi3rAhwFcnsK15GAiIRE5DScj5cKnCEWcIYUPuz29MSt128UkWOxxLoN5/14K566iHPPeoBOEanGMQzx4n+eCajqMI6hytdwhnJ+56ZHXbm/ISJl7vVVi8jF7u9LReQkt8HvxHmnk70DLTgK/Z3us3gfngZIRK70PJd2nIZuVDkT8GvgDPcdCwEfIfGjxM8XgM0i8jURqXDlOElEfioihRxb/fHT4so/0fqmiZ5bSohIpoisB+7BuX8/mkw5OB98XxJ3OE1ESkVkPAve7+O896vd/AUiMmoJwxj8J3CdiKx334eT3PNO1NYkZTJK6D0446mH3K+gRlVtxLEKewfO18CbcCarDuF0q69yj/09jvVSo4gcHcr7Bo6VRxNOg/Yzz7l+izOE8DzOEMsA4w97JKCqDTi9mc3AHZ5dFThf81043dY/4AzRJSvj34BrcRRDi3v+j+JUGoDrcBrbbpyX/o4kxeDJ+xecyeMjwL/iWCR14hh5/CfOC9/L2MMRP8G5F3U4Vik7ffv/C1jldofvUWc9yaU4Hwsv43wh/ydOLwocC7KD7r4HxroPPh7BabhqVfUJT/rfAjeKSDeO8r5znDI+j9OYtbsyxBpmdeamLsfpVR695/9Akvrq9sLeBLzevbbvAu9W1WdTuI6jXCUiPe41bcNR8OtVtd49xy4ca6CbXHn340wcHwvbgJVAo6r+2ZP+TzgWm504jfEvfcd9Gfic+zyvG6Ps23C+on/uGwb/R1fWne6Q0YPE51NXuts9OO/Id1X1oTHK/yDO/W8DVuNYnR3lbOAx9/5tw5nzPKa1Qe6w/pXAV91zrMKZw0k6X6qqL+JYrS0HnhZn6PcX7jHdx1J/kpTdB3wJx1y6Q5xh52RM9Nwm4lPue9KG807vBjarau/4h43Jf+Dc/wfccnfiGNIkRVXvxml/bnfrxlM479CEqOrPce7RbTjt3j04lngTtTVJEZ3UXKRhGEZ6cIdDa4F3jKMYjXnCnHbbYxjG/MAdOit0h6ePzqv6e/nGPMSUkGEYs4FNOBZXrTjDq3+lyc3djXmGDccZhmEYM4b1hAzDMIwZY3Y5skuBkpISXb58+UyLYRiGMafYvXt3q6qOtf5wxphzSmj58uXs2jVbowYYhmHMTkRkIm8oM4INxxmGYRgzhikhwzAMY8YwJWQYhmHMGHNuTsgwZiM6Mgz9fZCdi4Rm/2s1PDxMbW0tAwMDMy2KMcVkZWVRU1NDOJwsMMDsY/a/LYYxi9GBfvS5p9H2tliaFJciJ69CMrNmULLxqa2tJS8vj+XLl+P4MDXmA6pKW1sbtbW1nHDCCTMtTkqkTQmJyA9xnNk1q+rpSfYLjtO9N+B4Qb5aVf+ULnkMY6rRoUG+/p2fUN/j7008S1XeE1z38fchodn5NTowMGAKaB4iIhQXF9PS0jLToqRMOueEbsUJjz0Wr8fx5LsSJwLk99Ioi2FMOfricyzNChD0hZoJoizLDKAvvTBDkqWGKaD5yVx7rmlTQqr6CONHCbwc+IkbHnYnUChO5EnDmPXo4ADa3MhFmYOjIpkFgIuyBtHGOnR4eCbEM4w5w0zOCVWTGBuo1k0bFTteRK7B6S2xdOm4gT8NY1rQhlpQJX9kkA2LQjze40RCCwIbFgl5w4MQjaKNdciS5TMs7WgWP7WbSNOhaT1n8IKLR6Xdc889XHHFFezbt49TTz113OO/+c1vcs0115CTkzNuvrG49dZb2bVrFzfddNOkjp/qcgyHOWGiraq3qOoGVd1QWjrrvE4YCwyNRtF6N+ZgVhavztFYKFEBXpOjkOUYJWj9YeaDk+BoTxfRA/uJPvuU87+na0rK3bp1K+effz5bt26dMO83v/lN+vr6puS8xuxhJpVQHbDEs13jphnG7KatGYbcoJ8lZTQPRWK7RoD24QiUlDkJ/X3gsZybi0R7uqD2IAwPQ2am87/24HErop6eHv7v//6P//qv/+L22+PR2CORCNdddx2nn346a9as4dvf/jbf+ta3qK+v58ILL+TCCy8EYNGiRbFj7rrrLq6++moA7rvvPs455xzWrVvHRRddRFNT09jXFo2yfPlyOjo6YmkrV66kqakppXKuvvpq7rrrrti2V6avfe1rnH322axZs4YvfOELx36DFggzqYS2Ae92Y5SfC3S64bgNY1ajdfFR5MCifHboIrIFloWjZAA7JI/Aovx4/vqUI9LPTlqbIRiCcBhEnP/BkJN+HNx7771ccsklnHzyyRQXF7N7924AbrnlFg4cOMCePXvYu3cv73jHO/j4xz9OVVUVDz30EA89NH6w1fPPP5+dO3fy5JNPsmXLFr761a+OmTcQCHD55Zdz9913A/DYY4+xbNkyysvLj6kcPw888AAvvPACjz/+OHv27GH37t088sgjKR+/kEinifZW4AKgRERqgS8AYQBV/T5wP4559n4cE+33pksWw5gqtLcH7Yjb27QPDPNM5wAXLinmDSeW8ssXGnmsoZPLh0fIDTuvl7Y2owP9SFb2TIl9fAwMOD0gL6GQk34cbN26lU984hMAbNmyha1bt7J+/XoefPBBPvzhDxNyF/0uXrz4mMqtra3lqquuoqGhgaGhoQnXy1x11VXceOONvPe97+X222/nqquumlQ5Xh544AEeeOAB1q1bBzi9vhdeeIFXvvKVx3QtC4G0KSFVfdsE+xX4SLrObxjpwN+rqesZICMY4NyqQgA2VxWxvb6Dxxs6uXBpcfy4hlrkhJXTKuuUkZXlDMF5V+CPjMTmvSbDkSNH+P3vf89f/vIXRIRIJIKI8LWvfS3lMrymyF7PDx/72Me49tprueyyy3j44Ye54YYbxi1n06ZN7N+/n5aWFu655x4+97nPpVxOKBQiGnVmBKPRKENDQ4CzaPQzn/kMH/rQh1K+noXKnDBMMIzZgEZG0Mb6hLTTS/L4wqaTWJzlNNAVuZmcWJDNn5oT50u0oRaNRpmTlJRBZMRRRKrO/8hIfN5rEtx11128613v4uDBgxw4cIDDhw9zwgkn8Oijj/La176Wm2++mZGREcBRWAB5eXl0d3fHyigvL2ffvn1Eo9HYcBpAZ2cn1dXVAPz4xz+eUBYR4YorruDaa6/ltNNOo7i4OOVyli9fHhtG3LZtG8OuSf7FF1/MD3/4Q3p6egCoq6ujufn4hi/nK6aEDCNVmhqcxtdlYCSCqpIRTHyNrjy5ko+c6VtKMDQErWNPkM9mAovyoWaZ0xMaHHT+1yxLmPc6VrZu3coVV1yRkPaWt7yFrVu38oEPfIClS5eyZs0a1q5dy2233QbANddcwyWXXBIzTPjKV77CpZdeyubNm6msjC8xvOGGG7jyyitZv349JSUlKclz1VVX8dOf/jQ2FJdqOR/84Af5wx/+wNq1a9mxYwe5ubkAvO51r+Ptb387mzZt4owzzuCtb31rggI14shcMx/dsGGDWlA7Y7pRVaK7d0BPvCH57p6DLAqHePfq6pTKkIIiAus2pkvEY6Lp5z+lZJqXOyRbJ2Skh3379nHaaaclpInIblXdMEMijYk5MDWMVOjqSFBADb2DvNTZzxtPTN6QH+7u567nG3n3qmqKszMA0M52tLcbyc2bFpHH48jp6yn3NVKGMRPYcJxhpIDfIGFHfTshETZWFCTNnxcOUd8zyGMNHQnpXvNuwzBMCRnGhOjQINocn88ZGImwu6mLtWV5MTNsP4VZYVYXL+Kxxk5GPAYJ2lSPjowkPcYwFiKmhAxjArSxDjSuSJ5s7mIwEmVzVVFCPqlMnBvaXFVE73CEvS2eCelIBG1KtLAzjIWMKSHDGAdVj584l/XlBbxrVRVL8zzrZLJzkJWrICO+qPOkohxKssNsr/cNyc0Tf3KGMRWYEjKM8TjSBgP9CUkZwQBrS/MTFktKVQ0SCCBVNbG0gAiXLC/lnMrCRKXT2wOd7WkX3TDmAqaEDGMctC4x3MH9LzWPMjYgEEAqnKE4qaxx/Ku5nFmWz9kVBaMCjc15f3JTQGNjI1u2bGHFihWsX7+eN7zhDTz//PPjHuN1EHosPPzww1x66aXj5tmzZw/333//MZd9wQUXMBXLRqaqnLmGmWgbxhhofx96pDW23T00wh9qj4yeCyqrQMKOGbZkZiHFpajHuWffcIQnGjvZWFlAdijolN3SjA4NIhk+n2yzkK/ffBt1jaPDRVdXlHLdh94+qTJVlSuuuIL3vOc9MQ/af/7zn2lqauLkk08+Lnkny549e9i1axdveMMbZuT8CxXrCRnGGGhD4lzQ4w0dRBQ2uX7ijiJVid4RpDpx+8jAMPe91MzuJo8rH42iDXMjcsnymgqCPq8QwWCA5TWTD4T80EMPEQ6H+fCHPxxLW7t2La94xSvo6enhNa95DWeddRZnnHEG9957b9Iy/vVf/5UzzjiDtWvX8ulPfxpI7E20trayfPnyUcc9/vjjbNq0iXXr1rF582aee+45hoaGuP7667njjjs488wzueOOO+jt7eV973sfGzduZN26dTE5+vv72bJlC6eddhpXXHEF/f39o87xP//zP1x55ZWxbW9P7G/+5m/YsGEDq1evHjPEw1hhKlpaWnjLW97C2Wefzdlnn80f//jHsW7xnMF6QoaRBI1GEpRQVJUdDR2sLMyhLCfee5G8fMjzua8pXAw5udDXC0BNXhZL8rLYXt/OeVWFsaE5rT+MLl2OyMx/C950612j0s5cvZLzz17LBZvPYvvuvyTsi0ajlBY7yrinr59b7/x1wv6PXv3Wcc/31FNPsX79+qT7srKyuPvuu8nPz6e1tZVzzz2Xyy67LGFI8ze/+Q333nsvjz32GDk5OTH/cqlw6qmn8uijjxIKhXjwwQf57Gc/yy9+8QtuvPHGhIipn/3sZ3n1q1/ND3/4Qzo6Oti4cSMXXXQRN998Mzk5Oezbt4+9e/dy1llnjTrHRRddxDXXXENvby+5ubnccccdbNmyBYAvfelLLF68mEgkwmte8xr27t3LmjVrUpL9E5/4BH//93/P+eefz6FDh7j44ovZt29fytc+GzElZBjJaGlyHHW67GvroWNwhMtXlCdkk6olo+Z7RASpWoLufzaWtrmqiDuea+DFzj5OKnT8izE4AG2tx+UIdDrIX5RLXm4uXT29sbS83Fyys9IzlKiqfPazn+WRRx4hEAhQV1dHU1MTFRUVsTwPPvgg733ve2Ohvo8l3ENnZyfvec97eOGFFxCRmNNRPw888ADbtm3j61//OuB46j506BCPPPIIH//4xwFYs2ZNUgUSCoW45JJLuO+++3jrW9/Kr3/961g8ojvvvJNbbrmFkZERGhoaeOaZZ1JWQg8++CDPPPNMbLurq4uenp5Jz5XNBkwJGUYS/J4NskJB1pTmsarE87KHQlCWfEhKKqrQl16AqBN19czSPO57sYnt9R1xJYTTG5JZoITG67lkhMNce83b+OK3fsTISIRwKMQnP/Q28hc517EoJ3vCno+f1atXJ0Qk9fKzn/2MlpYWdu/eTTgcZvny5QmhGsbDG1phrGM+//nPc+GFF3L33Xdz4MABLrjggqT5VJVf/OIXnHLKKSmd28+WLVu46aabWLx4MRs2bCAvL4+XX36Zr3/96zzxxBMUFRVx9dVXJ5VzrDAV0WiUnTt3knUcYTRmGzM/DmAYswzt7kK7Ei3gVhTm8O5V1QS9ZtkV1UgwmLQMCYWR8viXezgY4OyKQgZHokQ95tp6pBXt75viK5h6CvJyOefMVYjAxjNXxRTQZHn1q1/N4OAgt9xySyxt7969PProo3R2dlJWVkY4HOahhx7i4MGDo45/7Wtfy49+9CP6+px7d3Q4zhtaYSwl5w3RcOutt8bS/aEiLr74Yr797W/HzOuffPJJAF75ylfGPHs/9dRT7N27N+l5XvWqV/GnP/2JH/zgB7GhuK6uLnJzcykoKKCpqYnf/OY3SY8dK0zF6173Or797W/Htvfs2ZP0+LmEKSHD8OE3n37+SC9dg6Nd7UjVknHL8RssvPHEUj64ZgmBOWqu/dpXnsOJS6t53auO3xO4iHD33Xfz4IMPsmLFClavXs1nPvMZKioqeMc73sGuXbs444wz+MlPfsKpp5466vhLLrmEyy67jA0bNnDmmWfGhsyuu+46vve977Fu3TpaW1tHHQfwqU99is985jOsW7cuFrMI4MILL+SZZ56JGSZ8/vOfZ3h4mDVr1rB69Wo+//nPA45hQU9PD6eddhrXX3/9mHNbwWCQSy+9lN/85jcxo4S1a9eybt06Tj31VN7+9rdz3nnnJT12rDAV3/rWt9i1axdr1qxh1apVfP/730/hbs9uLJSDYXjQkWGi2/8QG0YbjkS5ced+Ti7K5V2r4m55pKiYwNqJveJH/7QT7epMSOscHGZRRijeqwqFCWx61Zi9qnSQzNW/MX+YS6EcrCdkGB60sT6mgAD2tHTRPxIdbZZdPX4vKJbP1xs61NXPl3a+yDNtPfHEkWFoaZy80IYxhzElZBguqjpqaGx7fQdlORmsKMiJJ2ZmQXGKAeHKyp1IpC41eVnkZ4bYXpfotsdCPBgLFVNChnGUjiOxtT3gBKY73D3A5qqiRD9xlTUpr+2RQDDm0gccf3KbKgt5oaOP5r7BWLp2d6LdncmKSBtzbSjeSI259lxNCRmGi78XdLCrn6xggPXlnsWoIqNCNkyE34BhY2UhQYEdfu/a09gbysrKoq2tbc41WMb4qCptbW1zyoTb1gkZBqCDAwn+3gDOr17M+vK4vzcAKSlHMo/tBZfsHGRxScwPXV5GiDWl+exu6uLSE8sIBlwPCs0N6IpTEM/wXbqoqamhtraWlpbRPuGMuU1WVhY1NTUTZ5wlmBIyDFw/cZ5ewXAkSjgYSFBAkLpBgh+pXprgDPWS5SW88cTSmAICIBpFG+uQJfA7YgcAACAASURBVMsndY5jIRwOc8IJJ6T9PIYxETYcZyx4NJoYuC6qyr/vPsCvXkrsGZGTCwVFTIrFxZCVHdsszs6gMHN0j8cC3hkLDVNChtHWDENxI4H97X209A9RlZvoG02ql47yE5cqIokB7wDa+oe4Ze9hDnZ5vDD390F726TOYRhzEVNCxoLHbxCwvb6d3LDjKy5GMIiUVx3XeaSiGjxWdbnhIAe7+tlR7zPXniMeFAxjKjAlZCxotLcH7YiHAegYGObpth7OqSggFIi/HlJehYSObwpVMjKRsrgX7qxQkPXl+exp7qZ3OL5AVlub0YHRMWoMYz5iSshY0Ph7HTvd0N3njgpcNzmDBD/+cjZVFTGiyhONPnNtX0A9w5ivmBIyFiwaGXHc9Hg4r7qId5xWxeKsjFiaFBQii/L8h0+O/ELIjZdVmZvJiQXZ7KjvSPSu3VCLuiEJDGM+Y0rIWLg0N0Ik0Tt2XkaIM8sSI6X6/b8dDyIyysz7wiXFnFtVmKCEGBoC37olw5iPpFUJicglIvKciOwXkU8n2b9URB4SkSdFZK+IvCGd8hjGUVSVaN2hhLR79zclOhYFCGdA6dQGnZPySgjG55dOK17EhUuKE+agALT+kP9Qw5h3pE0JiUgQ+A7wemAV8DYRWeXL9jngTlVdB2wBvpsueQwjge5O6IkHMGvsHeTRunaaPP7cAKSyGglMbYgFCYaQikRLu+FolN1NnXQMxENNa0c72tvtP9ww5hXp7AltBPar6kuqOgTcDlzuy6PA0bGPAqAew5gG/GbZO+rbCYmwsaIgIX2qDBL8+MvtGYpw+7MN7GjwGSjUm4GCMb9JpxKqBrxveq2b5uUG4J0iUgvcD3wsWUEico2I7BKRXebryjhedHgIbY7H7xmMRNnV1MXasjxyw/FhMikuRTxeDqYSyV2EFMa9LxRlhVlVvIjHGjoY8RgkaGM9OjI6qqthzBdm2jDhbcCtqloDvAH4b0niI19Vb1HVDaq6obQ0xTguhjEG2lAHGm/o/9TUyWAkyuaqRJc8Uj11BgnJ8Bs8bK4qomc4wl9aPfNSkRG0uSGtchjGTJJOJVQHeMccatw0L+8H7gRQ1R1AFlCSRpmMBY5qdNTaoJxwkDNL81ia5/GOnZUNRcXpFaakDDLipuAri3IozgonDXhn/uSM+Uo6ldATwEoROUFEMnAMD7b58hwCXgMgIqfhKCEbbzPSx5E28HkjWFuazztXVScGrqteMmk/cakigQBSGf9OC4iwqaqQgUiUgZG4BwV6u6GrI0kJhjH3SZsSUtUR4KPAb4F9OFZwT4vIjSJymZvtk8AHReTPwFbgarVPPiON+HtBz7f3MhjxLQoNBBKioaYTqaoBj7J7Rc1irl2/nCxfCAkL/23MV9IaT0hV78cxOPCmXe/5/QxwXjplMIyjaH8f2hbvaHcPjfBffznMedVFXLYi7tNNyiqQcEayIqYcycxCiktjAfWCrkIaGIkQECEj6HwnaksTOjSIZGSOWZZhzEVm2jDBMKYNvz+2xxs7iSicW5keP3Gp4jeA6BgY5p93vsgTjZ3xRI06BhWGMc8wJWQsCDQaSVBCUVV21rdzUmEOZTme3sWifMgrSFJCGilcDNk58c2sMKU5GWyvb08wSHAC3pk/OWN+YUrIWBi0NMFw3BvBs0d6aR8cYbPPW3ZgGgwS/Dj+5Pzm2oU09Q3xUqfHiGJwANpaMYz5hCkhY0Hgn9jf39FLfkaI1cUe79ihEJRVTrNkDlJeBR73QGeW5pMdCrDdAt4Z8xxTQsa8R3u6UJ+J82Uryrl2/XKCAY9ZdkU1EpxaP3GpIuEwUl4R284IBji7ooC/tHbTOxz3mKBHWtH+vpkQ0TDSgikhY97j7wVFos48y6KMRONQqaqZNpmS4feg8MqaxfzdWcsTXAmB9YaM+YUpIWNeoyPDaFPc7c1wJMq/PPYif/R5JZCiYiRn0XSLlyhDXj6SHzeKKMwMU7Uoa1Q+bahDI5FR6YYxFzElZMxrtLEeovEG+88t3XQOjVCWk7gOaLrNssfCL0ffcITb9tXzdKsnpMPIMLQ0YhjzAVNCxrxFVUcNXW2vb6c0O4OTCuMm0WRmQckscYxbWgGhcGwzMxTgxc4+/ug3UDAPCsY8wZSQMX/pOAJ9vbHN2u4BDnUPsLmqMNFPXGUNSZy3zwgSDCKVcZdBQRE2VRbyfHsfLX1DsXTt7kS7O5MVYRhzitnx5hlGGvAHhNtR3044IGzwBq4TSWj0ZwP+IbmNlYUExJHfixkoGPMBU0LGvEQHB9DWpoS0C5YUc9UplWR7nINKSTmSOXryfyaR7BxkcTyiSX5GiDUleTzR1MmQx9mqNjWingW4hjEXSasDU8OYKZzAdYkO2UtzMij1GyRUzw6DBD9StQQ9EveOcH51EQWZYUaiSsZRHRqNoE31SM2ymRHSMKYA6wkZ8w6NJgauU1V++UIjh7oS4wiRkwsFRcxKikscgwmX5QU5vGlFGTlhf4iHQxbwzpjTmBIy5h9tLTA0GNvc39HH9voOWvqHErJJ9dJp9xOXKiKBUXNDUVWeO9JDU2/82ujvcwwwDGOOYkrImHckM8vOCQVZU+rxExcMIuUz4ycuVaSyGjxWe0ORKD95pp6HDicqHa07NN2iGcaUYUrImFdoXw/a3hbb7hgc5unWHjZWFhAOxKu7lFchnvU4sxHJyERK48H2skJBzirPZ09LF73D8QW42tqMDg7MhIiGcdyYEjLmFX6z7McaOlBg06jAdTPrJy5V/IYTmysLGYkqTzQmOmT1X7dhzBVMCRnzBo2MoI2J0UdzQkE2VBRQnB23ipOCQmRR/nSLNznyCyE3PoxYuSiLE/Kz2VHfQdQb8K7hMBq1gHfG3MOUkDF/aG6EkZGEpFfULOaqUxLnfmaLn7hUcALe+XpD1UWMRJX2Ac8aoaEhaG2eZukM4/gxJWTMC1SVqM+f2osdfUT85svhMHjmWeYCUlYJwfiSvjWleXz23BUJvTsArTcDBWPuYUrImB90d0JPV2yzsXeQ7/350OiQDZU1SGBmAtdNFgmFkIqq2HZQhKAIkagy6PWg0NGO9vbMhIiGMWlMCRnzAr9X6Z0NHQRFOKssce5nLg3FefEbUgxFonz58Rf5/aG2hHTzJ2fMNUwJGXMeHR5CPfF1BiNRdjV2srY0LyF6qhSXIlnZMyHicSO5eUhh3LtDRjBA9aIsHmvoYMRjkKCN9WhkJFkRhjErMSVkzHm0oQ48DfGTTV0MRKJsrvKbZc/NXtBR/OG/N1cV0jMc4S+tniG4yEhCJFnDmO2YEjLmNMkC1+070kNVbibL8j29nqxs8HimnpOUlEFG3BhhZVEuxVlhticJeGf+5Iy5gikhY25zpBUGEh2Tvmd1Ne8/oyYxcF3VklnrJy5VJBBAKuNzQwERNlUV8nJnP41ef3K93dDVkaQEw5h9WCgHY07j7wVFVQmIUJDpcckTCMy6wHWTRSpr0IMvxbY3VhRStSiLMl+ICq07jMxWD+GG4cF6QsacRQf60baW2Hb30Ahf3PkiT7d2J+ST0goknOE/fE4iWdlISVlsOycc5OSiXAK+Xp62NKEeT+KGMVtJSQmJyFdFJF9EwiLyvyLSIiLvTLdwhjEe/l7QE42ddA2NUJI9NwLXTRa/gcVIVLnvxWaeaOyMJ2rUMdgwjFlOqj2h16lqF3ApcAA4CfiHiQ4SkUtE5DkR2S8inx4jz1+LyDMi8rSI3Jaq4MbCRqORhEY2qsqO+nZWFOZQnpsZz7goH/IKZkDCNFJUDNk5sc2gON4hHj7clmCQoA21qJo/OWN2k6oSOjp39Ebg56raOV5mABEJAt8BXg+sAt4mIqt8eVYCnwHOU9XVwN+lKrixwGlphuF4kLpnj/TSPjgyyiw7UD33DRL8iEhCb0hE2FxVSFPfEC91eow0BvrhSFuSEgxj9pCqEvqViDwLrAf+V0RKgYkCmGwE9qvqS6o6BNwOXO7L80HgO6raDqCq5oHRSAm/n7Tt9e3kZQQ5vdgTuC4UgrKKaZZsepCKavC4HzqzLJ/sUIAdo8y1zZ+cMbtJSQmp6qeBzcAGVR0G+hitUPxUA95B+1o3zcvJwMki8kcR2SkilyQrSESuEZFdIrKrpaUlWRZjAaE9XWhnognyRUuLueKkcoIBj1l2RTUSnJ8GoBIOIx4FmxEMcHZ5AXtbu+kaintM0COtaH/fTIhoGCmRqmFCDvC3wPfcpCpgwxScPwSsBC4A3gb8QEQK/ZlU9RZV3aCqG0pLS6fgtMZcJlkAt+UFOawp9fuJmxuB6yaL3+BiU1UR68ryGYkmLlQ1f3LGbCbV4bgfAUM4vSGAOuCLExxTB3jfkho3zUstsE1Vh1X1ZeB5HKVkGEnRkWG0qT62PRyNsm1/E639Qwn5pKgYyVk03eJNK5JXgHiMLkpzMnjbqVUszkoMW64NdWgk4j/cMGYFqSqhFar6VWAYQFX7gIlme58AVorICSKSAWwBtvny3IPTC0JESnCG517CMMZAmxrA06Dubenmkbr2xABvzH0/camSzPy8rmeA2m7PlO3IMLQ0TaNUhpE6qSqhIRHJBhRARFYA466EU9UR4KPAb4F9wJ2q+rSI3Cgil7nZfgu0icgzwEPAP6iqmfMYSVHVURPt2+vaKc3O4KTCuMkyGZlQskCGbUsrIBTv+URV+eFTtdz/cuLcqQW8M2Yrqc7a3gD8D7BERH4GnAdcPdFBqno/cL8v7XrPbwWudf8MY3w626GvN7ZZ1z3Awe4BLl9RlsRP3MJwBiLBIFJZjR4+ADj+5M6tLOS3B1pp6Rui1HXno12daHdnwvCdYcwGUrWOewB4M47i2YpjJfdw+sQyjNH4A9dtb2gnHBA2VHgaVpF54ycuVfxDj+dUFBIQ2NHgM9c2AwVjFpKqddx9wOuAh1X1V6raml6xDCMRHRxAWxPnNTKDAc6tLCQ7FF8vIyVlSGbWdIs3o0h2DuIJU5GfGeKMkjyeaOxkyBv+u6kRHRlOVoRhzBipjll8HXgF8IyI3CUibxWRhfWmGzOKNtSBL0bOZSvKufyk8oS0hWKQ4Md/3Zuriogq1Pd4pm6jEbSxHsOYTaQ6HPcHVf1b4ETgZuCvAfNuYEwLqlG0odazrRzq6h8duC0nFwoXT7N0s4TiEvD0AE8syOb6TSexvCAxnLnWHbKAd8asIuXZW9c67i3Ah4GzgR+nSyjDSKC1BQbjJsf7O/r41pMHeaqtJyHbfAhcN1lEAqP8yWUGA6gqg54hOfr7oOPIDEhoGMlJdU7oThwz61cDN+GsG/pYOgUzjKP4J9S313eQEwpy6uLceGIgiFRUTbNkswuprAaPVaCqctOeQ9yzP3EuzfzJGbOJVHtC/4WjeD6sqg+p+Yc3pgnt60Hb40vHOgeHebq1m40VBYQD8eor5ZVIKJysiAWDZGQipfGAdyJCVW4mTzZ30TccX+CrbS3o4ET+hw1jehhXCYnIq92fucDlIvJm71/6xTMWOn4/cTsbOlDgXF/IhvkWuG6y+A0UNlUVMhJVX8A7Tep/zzBmgokWq74K+D3wpiT7FPjllEtkGC4aGUEb4+4GVZWnWns4ZXFuQvRUyS9EFuUnK2LhUVAEuYug15kvq1qUxfL8bHY0tPOKmqJYGHBtOIwuOxEJLIxFvcbsZVwlpKpfcH/e6DoYjSEiJ6RNKsMAaG6EkXhYAhHh4+uW0Tuc6IzTekFxjga80xf2xdI2VxVy27MNvNDexylH59GGhqC1ed7GWzLmDql+Bv0iSdpdUymIYXhRVaI+DwmqSjgYoNDrJTochtJyjDhSXgXB+ALeNaV5vOO0KlYU+sy1zZ+cMQsYtyckIqcCq4EC3xxQPmCLVY300d0FPV2xzabeQX70dC1vP7WKpfnxxlQqaxBPhFEDJBRCyqtiVoWhQIB1ZaOHK7WjHe3tQXLnd8gLY3Yz0ZzQKcClQCGJ80LdOKG5DSMt+L/SdzR00D4wQpEvVo5Uzu/AdZNFqpckmLarKn+oPUIoIJxfHV/Qq/WHkZWnzYSIhgFMPCd0L3CviGxS1R3TJJOxwNHhIbS5MbY9GImyq7GTNaV55GXEq6wUlyLZOcmKWPBIbh5SUIR2Ok5MRYSXOvs53NXPuZVFhNww6NpYj564ct6GQTdmP6nOCT0pIh8Rke+KyA+P/qVVMmPBoo11EI0vRXuyuYuBSJTNfrPsBeonLlX8BhubqwrpHo7wVGt3PDEyAk0N0yyZYcRJVQn9N1ABXAz8ASdUd/e4RxjGJHAC1yUOI22vb6cyN5PlnrkgsrJhcfEMSDiHKCmHcNyU/eSiXIqzwmyvTwzxEK07bP7kjBkjVSV0kqp+HuhV1R8DbwTOSZ9YxoKlvQ0G+hOSXrushNefULpgA9dNFgkEkKr4nFlAhHOrCnmps5/GXo937d5u6OqYAQkNI3UldDQISYeInA4UAGXj5DeMSeH3ayYinFGSx6pijwVXILDgAtdNFr/hxsaKAk5bnEvE1/OxgHfGTJGqErpFRIqAzwPbgGeAf02bVMaCRAf60baW2Hbv8Ai/PdBC99BIQj4prUA8w0zG2EhWNlIS/17MDYd4/xlLqF6UuMJCm5vQoUH/4YaRdlJVQj9S1XY3rtCJqlqmqjenVTJjweH/Gn+8oZPfHWwb7SGhysyyj4VkBhwdA8Mc6vIMe2o0wUWSYUwXqSqhl0XkFhF5jSzUgC1GWtFoxIme6hJVZUdDBysKcqjIzYxnXJQH+YVJSjDGpKgYfKbstz5Tx8+fb0wwSND6WsxBvjHdpKqETgUeBD4CHBCRm0Tk/PSJZSw4WppheCi2+dyRXo4MDLO5OlHhBBZw4LrJctSfnJdNlYU09A5ywNsbGuiHI20YxnSSanjvPlW9U1XfDJyJ47bnD2mVzFhQ+D0kbK9vJy8jyOnFefHEYAjKK6dZsvmBVFSDx2P2mWX5ZAUDbK9PtIqzgHfGdHMs4b1fJSLfBXbj+I3767RJZSwotKcb7Yw3hlFVwoEAm6uKCAY8ZtkVVbayf5JIOIyUxRV4ZjDA2RUF7G3pSjD80COtaH/fTIhoLFBSDe99APg74FHgDFX9a1VN5lnbMI4Zv0FCQIR3r67moqWJi1EtZMPxkSzgHcDLnYlKx8y1jekk1c/KNaraNXE2wzg2dGQEbaqPbY9Eo7QPjFCak5G4OLVwMZJj3p6PB8kvQPLy0W7nVS7LyeT6TSvJDSd6IdfGOvSEk8w7uTEtTBTK4VOq+lXgSyIyyq+Hqn48bZIZCwJtqodI3AR7b0s3tz3bwCfOWsaSPE/IBusFTQlStQR97unY9lEFNBKNEjo6ZzQ8DM1NUFE1EyIaC4yJekJHwzPuSrcgxsLD8RPnN0jooDQ7I3ExZUYmFJuDjimhrBJefB5GhmNJt+2rp3c4wgfXxBW91h9CTAkZ08BEoRzuE5EgzjzQddMkk7FQ6GyHvt7YZl3PAAe6+rlsRRmBBD9xNUjA/MRNBRIMIhVVaO3BWFpJdgZ/am6ltX+IkmzHE4V2daLdXUje6GB4hjGVTPhmq2oEOG8aZDEWGOoL372jvoNwQNhQXhBPFLHAdVOM30DhnMoCAjj334uF/zamg1Q/L/eIyDYReZeIvPnoX1olM+Y1OjiAtjbFtkeiyt6WbtaV5ZPjmSiXkjIk0yLJTyWSk4sUxS0PCzLDnF6SxxONHQxH4h4TtKkR9QzbGUY6SFUJZQFtwKtxwny/CSfs97iIyCUi8pyI7BeRT4+T7y0ioiKyIUV5jDmONtaBx2VMKCD848YTuXh5SUI+C1yXHpIFvOsbibKnxWMEG42gjfUYRjpJyURbVd97rAW7c0nfAV4L1AJPiMg2VX3Gly8P+ATw2LGew5ibqEbR+tpR6Y6llscsOCcXChdPn2ALieJSyMyCwQEAVhTmcOXJFaz2eqjAWTOk1UvNVZKRNlJdrHqyiPyviDzlbq8Rkc9NcNhGYL+qvqSqQ8DtwOVJ8v0zTliIgWOQ25jLtLXEGj+A/R29fPtPB2jtH0rIJuYnLm2IJAa8ExHOqSxMGAoFHMORjiPTLJ2xkEh1OO4HwGdwg9up6l5gywTHVAPemedaNy2GiJwFLFHVX49XkIhcIyK7RGRXS0vLeFmNOYDfIGF7fQfN/UMUZHg65oGgmQinGamsAZ+Sf7K5i98fSnRiah4UjHSSqhLKUdXHfWkjSXOmiDixmf8d+OREeVX1FlXdoKobSktLj+e0xgyjfb1oe7yR6xwc5qnWbjZWFBIOxqujlFciofBMiLhgkIxMpLQ8Ie3Fjj5+d7CVPk8MJ21tRgdtoMJID6kqoVYRWQEogIi8FWiY4Jg6wDv7WeOmHSUPOB142PVNdy6wzYwT5jf+r+rHGjqJatyP2VHMQ8L0kMyf3HBU2dXUGU9URRtGz+EZxlSQqhL6CHAzcKqI1OE4M/3wBMc8AawUkRNEJANn+G7b0Z2q2qmqJaq6XFWXAzuBy1TVvDPMUzQSSYjeGYkqOxs6OKUoN7ZIEkDyC5FFtkhyWigogty4T77qRVksy89me307UX/Au6gFvDOmnol8x13r2bwfeAhHcfUCb8EZTkuKqo6IyEeB3+KYPP1QVZ8WkRuBXaq6baxjjXlKcwOMeMIGAK9dVkxpTkZCNusFTR9HA97pC/tiaZurCtn6bAP7O/o4uSjXSRwahLZmKK2YIUmN+cpEJtpH7TVPAc4G7gUEeBfgnyMaharej6O8vGnXj5H3gonKM+YuqkrUZ5AQCgibqooSM4bD4JunMNKLlFehLz0fcyS7pjSPXY2d7uB7HK07jJgSMqaYiXzH/ROAiDwCnKWq3e72DcC4Fm2GkUB3F/TEF0K29A3xfHsvGyoKyPQaJFRUWwiBaUZCIUcRufN14UCAD61dOiqfdhxBe3uQXAupYUwdqc4JlQPeRRxDbpphpESy8N3bXmxiKJI4z2AeEmaGZPe9fyTCwa7+hDQz1zammlSV0E+Ax0XkBrcX9Bhwa7qEMuYXOjyENjfGtgcjUZ5o7GRNaT55nrVBsrgEyc6ZCREXPLIoDylItFC887kGbn26lpGox0ChsR6NHNfqDMNIICUlpKpfAt4LtLt/71XVL6dTMGP+oI114LGs2tPcxUAkyuZRZtmjh4CM6UOqEu//xopCuociPNXaHU+MjEDTRKszDCN1Ug3vjar+CfhTGmUx5iGqmuAnTlXZXt9OZW4my/PjkVPJyobFxUlKMKaN0jLYnwHDzsj7KYtzWZwVZnt9O2eWxU3mo/WHCVTWmEslY0qwSGFGemlvg/6+2OZAJEo4EGBTVWFCIyZVNThONIyZQgJBpDLuWSsgwqbKQl7q7KexdzCesacbujqSlGAYx4699UZa8Yfvzg4F+ei6ZWyq9AzFSQCpqMaYefwGChsrCwiJ8Hx7b0K6GSgYU4UpISNt6EA/2hZ3ODswEqFnyJnUTugFlZUjGZnTLp8xGsnKRorj/hlzwyE+e84KXlmTGFJDm5vQoUH/4YZxzJgSMtKG/2t5R0MH/7zzRToHE6N1mln27MJvIJKf6UwdRzxufNBoggsmw5gspoSMtKDRCNoQb6Siquyo72BZfhYFmR7v2IvyIL8wSQnGjFFUDD5T+f95uYVvP3kQ9fuTU/MnZxwfpoSM9NDSHLOyAnjuSC9HBobZ7HPTE7DAdbMOx59cTUJaYVaY2u4BDngXrw70w5E2DON4MCVkpAX/UNz2+nbywkFOL/GEjw6GoLxymiUzUkEqqiEQbx7WleWTFQywvT7RKs4MFIzjxZSQMeVoTzfa2R7b7h4a4bkjvZxTWUgo4DFIqKhCgikvVTOmEQlnIGVxZ6WZwQAbKgrY29JF95DHE3pbC+oxwTeMY8WUkDHl+L+O8zJCfGrjiZxXnTgUZwYJsxu/B4XNVYVEFB5v7ExIt4B3xvFgSsiYUnRkBG2qH5Vekp2R6CeucLF5Y57t5OUjeXFPCWU5mbxlZTnryvISsmlDLRqN+I82jJQwJWRMKdpUH4tLA/Bkcxc/eqqWvuHERsoC181+jga887KpqojFWYlBCBkehpamaZTMmE+YEjKmDMdPXOJQ3B/r2mnqGyQr5KlqGZlQXDbN0hmToqwSQonzdvvbe3ngQGtCmtaZgYIxOUwJGVNHZzv09sQ263sck95NlUUE/H7iAlb15gISDI5yqfRSZz+/O9hKa3/cBF+7OtDuLv/hhjEh1hIYU8Zos+wOQgHh7IqCeKIIUlmDMXfwD8mdU1mAADvNXNuYAkwJGVOCDg2innmB/pEIf2rqZF1ZPjnheLhuKS5FMrNmQkRjkkhOLlIUD7NRkBlmdUkejzd2MOyJjKtNDejIcLIiDGNMTAkZU4I21ILXtxhw4dJizvd5SLDAdXMTvyHJ5qpC+kai/LnFE/AuGkEbR1tGGsZ4mBIyjhvVaELgOnBCNrx2WQnVeZ5eT04uFC7GmIMUl4KnB3tSYQ4nF+Xi97ik9YcT/MsZxkSYEjKOn7YWGByIbdZ2D7C3pSvR6zLO3IL5iZubiAQS5vJEhGvWLGF9eUFixr5e6DgyzdIZcxlTQsZx4zfPfehwGz9/vpFo1KOEAkGkomqaJTOmEqmsxt/1GYkqh7xOTTEDBePYMCVkHBfa14u2xz0pdw2O8JfWbjZWFBIOxquXlFcgoXCyIow5gmRmISXlCWm/fqmZ7/35UMJiZG1tRj09Y8MYD1NCxnHh/+p9rLGDqMK5lYkxgvx+yIy5id9AYUNFAcNRZVeTx5+cqvmTM1LGlJAxaTSSaA0VUWVnfQcnF+VSmhN37SL5BQk+yIw5TEGRY2DiUr0oi2X52eyo7xgd8C5qAe+MiTElZEye5gbwrAvpGBgmHBQ2V1kvaL4iIqPM7DdXFdLSP8T+Dk9Ih6FBaGueZumMuYgpIWPS+IfiirMz+NTZJ7Kq2OMdOxyGtur0ngAAGcxJREFUsnKM+YOUV0EwvgB5TWkeueEgT7V2J+Qzf3JGKlhEMWNSaFdngq+wvuEIoYCQEUz8rpGKaiQQ9B9uzGEkFELKq2IfIeFAgI+vW8birETDE+04gvb1IDkWssMYm7T2hETkEhF5TkT2i8ink+y/VkSeEZG9IvK/IrIsnfIYU4e/F/S7g618+bEXGfbNA1jguvmJ/7kWZ2cgIkR9a8P8i5gNw0/alJCIBIHvAK8HVgFvE5FVvmxPAhtUdQ1wF/DVdMljTB06PIw2N8S2ByNRnmjsZEVhDmGPd2xZXIJk58yEiEaakUV5SEHi3N9jDR18fdfLRDzrw7SxDo2M+A83jBjp7AltBPar6kuqOgTcDlzuzaCqD6nq0dnMnYC5V54DaGMdeHo8e5q7GIhE2Wx+4hYUfoOTvIwQzX1DPNXmmRsaGYHmxmmWzJhLpFMJVQPeMZtaN20s3g/8Jo3yGFOAP3CdqrK9vp2K3ExOKMiOZ8zKhsXFSUow5g2lZRCOm+KfujiXoqww2+sSQzxE6w6ZPzljTGaFdZyIvBPYAHxtjP3XiMguEdnV0tIyvcIZibS3QX/cFLe+d5C6nkE2VxUm+IWTqhpEZkX1MtKEBIKOKx+XgAibKgt5sbOPxt7BeMaebujuTFKCYaRXCdUB3tnLGjctARG5CPh/wGWqOujfD6Cqt6jqBlXdUFpamhZhjdTwGyRU5WbysXXLOKvMsxhVAqOicRrzE7+BwsaKAoIi7KhvT0g3c21jLNKphJ4AVorICSKSAWwBtnkziMg64GYcBWQr22Y5OtCPtiY+JhFhWX42WSFP4LqyciQjc7rFM2YAycpGiuMfhosyQlyxspyzKxKNFrS5ER0e8h9uGOlTQqo6AnwU+C2wD7hTVZ8WkRtF5DI329eARcDPRWSPiGwbozhjFuD3B/bHunZ+/nxDgjUUmFn2QsNvgHJuZSE1eb7ouRpFG0YNhBhGeherqur9wP2+tOs9vy9K5/mNqUOjiYHroqo8UnuE/MwQwYDHvX9uHuQXJinBmLcUFTuGKAPxkA4NPQP8uaWbi5eXxOYKtf4wumSZzRUaCVhtMFKjtQk8wynPt/fSNjCcxCzbAtctNERkVO/3YNcADx5q42CXJ6TDQD8cacMwvJgSMlLCP7G8vb6DReEgZ5TkxRODIaS8cpolM2YDUlkNnoXK68rzyQoG2O43ULCAd4YPU0LGhGhvN9oZb0yODAyxr62HcyoLCXmG4qSiCgmaO8KFiIQzkNKK2HZmMMCGigL+3NJNz1DcY4K2taAD/cmKMBYopoSMCfH3ggIibKoqTBK4zgwSFjL+gHebKguJqPJ4Y+IaIesNGV5MCRnjoiMjaFN9QlphZpg3r6ygyOM1WQqLkFzzlrygySuARfH1YuW5mZxevAi/rwRtqEWjEQwDTAkZE6BN9RCJNxgvd/bxUkffKDcsFrjOEBECvt7Q1afX8JqlPvdNw8PQ0jSNkhmzGVNCxpj4/cQB3P9yC3c815D4dZuRASVl0yqbMUspq4RQ4rygqlLfM5CYZh4UDBdTQsbYdLZDb09ss6FngJc7+9lUVUjA6yeucgkSsKpkgASDo1w2PVrXzjd2H6CtP27ir10d/P/2zjxIjru6458397n3arW7um3ZlmxjxyzCGJMQrFDGRWES7kAwhRMCwUVSQBJS5IKqVEKoooqEI3EMxFCUITgcTgrCTaiAD8mXZNmyLdsSOlYrraQ9ZnZ3do6XP7p3pmd2RprV7uzM7r5P1dT29vT0vO7p7vd7v9/7fZ+mJio/bqxB7Mlh1KQyCvrliTECPimXZBFBBqwCh1Gi8np4UW8SAR4YLlfXtmjIAHNCRg10NoOeLunEzeTyPDwyzrW9bcSDHp247l4kHKm2C2ONIrEE0lkaB+oIB7myJ8FDJ8fLKu/qyDCayzbDRKOFMCdkVEWHj4GWHhgn07ME/T5uGKhIy7bCdUYVKtP1bxjoJJ3Ns++0p+BdIY+ePIGxtjEnZMxDtVwnDmBLe5S/uv5SNrV5CtdFY9DRtczWGSuCnl7wKKlf2hGjNxrisVPl40B64qgVvFvj2PR2Yz5nRiFTymaayuaJBHxl6gjgREGmE2dUQ8SHDGxEDx9y/xfefdWGsrllAEylnQQYa8ysWSwSMuZRmZDwzWdP8ulHDpe3WH1+pG9gmS0zVhLSPwieRkpvLETAJ/MiH0tQWNuYEzLK0Okp9Oxo8f+J2Rz7Rie5pCNWXr67bz0SDFbbhWEAIOEI0tNXtu6pMyk+ufcFpnOlCdA6OoJmZio/bqwRzAkZZVRGQQ8Nj1FQeFl/RckGU0gw6qBSTy4ZCnBqapa9Xj05VSt4t4axMaE1jKo6NV5Sk2h60vnrqfeSV+X+4TEu64zTGwsV10tbO5Jsq7ZLwyinvRNicWfsB9iQjLApGeGXJ8a4cbCzVPDu6GEK6UmIJ5FEwimOGInamOMawJzQGkFzWdfZpCA1CWl3OV8uJPmpvS9wIp0pWzeeyfGpvS/wwaGtgEVBRv2ICDK4CX32qeK6idkcY5kcf/rzp8u2HYiH+eDQ1pIklD/giOImkpBIOsvxJBKwx9Zqwn7NVYZqAaanIJUqRTfpVFnp5fOxKepnJK3kKbVA/Sibo+4E1UAQ1vXV+LRhzEf6+tHnnyk2eC5LhHgok4Va19gc+Rw6MQYTjtJC0TlFokgiWR41RaNWNnyFYk5oBaPZ2erRjWdW+kJ5RWiGPRXrfMDuiBMdSf8g4vPP+5xh1EICQWRdvzMBGnh1JMNewHuVCqVr7ILMTDuF8UZPlRyTz+9GTW60NOekLHmm5TEntALQQgGm0+XRTWoSZuu8aatQUEUBvwhHJqb56dEzDKcynJkpUN5ChaGEkMy6TsgK1xkXgQxuLDqhtlyGXYkAe1Iw1xmcQ/j8cI7+meMMJMJc1hlns3di9IUo5NHJcZh0Eh6KzikccR1SouiYiMUsamohzAm1GDqb8SQKuE4nnS6T0Fko2UKBoxMzDKczDKcznEjNcDKd4W1XDHB1b5J8QRlJzzKYjDAUztMmeb41JuQAn8DuuEIkgnT1INHY0h2ssWaQRBvS1uF0r0Ui7CbL3rSAOg2dl8eVs+pjOD3DE6OTZPMFNrdFyeQL/Nu+o/THw/QnwgzEw6yPh4kE6ozGMzNO+veZ056oyQexRLErr+ikPAoPxvJhTqhJaCHvZAylUo6kfTrlRDfZ2Qt/uAZ5Vc5Mz3Ii5TibzW1RdnYnGM/k+NzjvwIgGvAxEI+wq7+jOHt9W0eMP9+1DYBCagKOHeFoLMCDU8JQVEmSg55BZOv2xR+4sWaRbdvRx/ZAzzrajh1hKOpcY7tiymsTOdiwGV+ijUy+QL7guIypbB4RePTUBPcPlxpib75sPbv6O0jN5nhufIqBeITuaLCsxEhNCgVITRRLSRSdUyhcJWqKW5mSBmNOqMGoqtNtVhndTKVhEZpZ6WyOmVyB7miIvCr//OgRTqYz5Nyb1we8alM3O7sTdEWC3H7VBgYSYdpCgfOmvfoSbRQ2bGb3yRFGsll2dwWhbyOBoestLdtYFNLRhWzfge/QwSrX2CA+tzR42O9zwiOgMxLk/dduRlUZy+ScaD41U9QwfGF8mq886YigBn3C+rgTLf3mpm56oiFUtf4079kMejYDZ0dLjkl8jiOqiJoIhS19fImQlSYeODQ0pHv37m22GVXRfB6mUqUxm/QkmkrBEsjV7z89yZHJaYbdKGdiNseOrji3X+2M0dxz8ASJYKDYbdEXCxFYTAsuEHRuvPYuZP2AdcMZS4amU07Z+PExp0GWy130vrKFAiPpWU6kZ4r3xolUhg9ct5meaIj/O36Onx0949wX8Qj9iTD98TDrYqH6oqZaBINIoq08aorHWzppR0QeVtWhZttRiUVCF4GqOgKf1aKbRexzMpt3byRn/CZfUN6x06lS+fPjZ/nVxAx98RDbO2P0x8NlA7dvu+IiddxEnJZePAmJhPs3aS09o2FIPIFsuwzw3Etz3dFptwFX570U9PnYkIywIVmqaeVtWPdGg2xrjzGczvD0uTO4HQX83Y2XEfYL+05PMJbJuU4qTCJU5yMxm0XPnYFzZzxRk0A05kkfd6OmcMTupfNgTugCaC5XSn32Rjf5xbXeTqVnGZnKcF1fOwD3PnOSBz1SJu2hABuTkWJ3wjt3DhIL+PH7FtN6C3n6vNucdNZYa7fejNWNiEAk6sz96e4tri/rVUinSg2+7IV7FbwP/Mu7ElzelQAgV1BOTWUYnZ51uvyA/aMpHvWUl0iG/Gxtj/FOt/E3lsmSCAbmKchXRRWm0uhUGjhZck6BQLFxVxY1+e3xC+aEilSVsEmnnImfi9knzk3x7Lk0Dw6PMZzOcHpqtjhHYntnnGQowFU9Sda7XWn98UhZ9VJwNLfqRnzORV4R3Vj2j7FSEL8fku2QbC9OGCiOr5ZFTSnHWdUxrBDwCQOJCAOJUtT09h0DvO6SdZxMl7ry/B5/86UnjjGcztAXCxejpS3tUba2L6B7OpdDx8/B+DnKrCyLmtauVNGadEL1StgshGyhUOyTnkuDHk5neO81mxhMRBjL5DgyMU1/IsLVPcli3/Scs9nRnbi4L57L6Ekki06HqGX0GKsPEYFwxOne6uoprtdCwc009TYgJ2G2vkzTZChAMhRge2d83nu7N/VwLDXDidQMz49P8cipCa7tTRad0N0HjtERDhbv5/WxMEF/nffe9BQ6PQWnR6pLFXnGm1azVNHqPTJqSNikJssKti18n8rZmWzR2ezoSrAhGeGFsWnu3O8oUId8Qn8iwjW9bQTdMH6or42XrG+/+IPx+YraWUUdrUQSCYYu/FnDWMWIz1fSl/Os12pRUzq1oDl3V/cmubo3Wfx/Kptn1lUkyeQLjGdyHDybJusONglw85YebtrcQ65Q4OmzafoTETrD589KLbIGpYpWrRMq7H/EGThchITNTC5PtqAkQwEmZ3PcfeA4w+kMmbyzTwFiAT8bkhE2tkW47cpB+uNhuiLz5yssKMSOREtOphjd2Cxvw1gIEgo7PQWd3cV1jvrIVCkBYq5HpM6GaSzoJ+bmj4f9Pj5w3RYKcw3TVIYT6ZliwtBIepYvHXBKVET8vmLm6vX9HWVdgnVRS6qopxffzmsWtq8Wo6FOSERuBj6Nk/V/l6r+Q8X7YeDLwIuBM8BbVPXwkny5KhQKVVWhoaTYW9pc2T86WZzoOZzOcHYmy42Dnbz+0j4nKUCEF/e1MeCmeq6Ph4sDnNGAn6t7kvO+57z4/dWjm4DpXRlGIxCfz8lYiyeQdf3F9ZrNlhxTMeM1BYULd9H7ROiJhuiJhsqipnXxEHdcu5nh9Awn0hmGUxkeHpngqm5nm6fOpPj2oZFiV97cc6WyEXu+Z9iHXtM7b/1Ko2FOSET8wGeB3wKOAXtE5D5VfdKz2e3AOVW9VETeCnwCeMuSGBBPwNnRqqrQPpQQBb7x9DDxoJ9btq1DRPjOc6eYyOTojYXYmIzw0vXtXNLh9P36fcL7rl1ECYNorCK6WZuDkIbRikgwCB1dSEdXcZ3TnT/t6cpznVSdivRBn48t7VG2tJemUqir2QgQCfgYTEYYTmU4MJoqrv/w0FbWx8M8PzbFcDpDZwBGqKFsn7jIseQWopGR0C7gkKo+DyAiXwNuBbxO6Fbgb93le4HPiIjoEsyglUQSBXaHHcVeLwWEw+kspzIpdnaXBiPfd80m2kOB+gcWq+FNx/QOMFo6pmGsKMRVSyAWRzwBh+aynrEmT/p4HYlNIlJ0JVvbY8UEh9l8oZih1xt1xnkPnEnxv8fOzn2ybD9zyvYSX2DvSwvSyCfjIOCtFX0MeGmtbVQ1JyLjQDcw6t1IRN4DvAdg06Y6oxH3x2nLZRjyKPYKcHkE3pDI0bHz8rJIpCe6wEH+apM8bWKaYaxqJBCE9k6kvVTyvjjFozJqqnOKR8jvY1NbtChHBPDabb28YkMnxw8+w4+m/ByfdcpflCnbJ8wJLQuqeidwJziyPXV9KBZ35stUKPYGBN7UriQjC3AWgaAnDdqNbmIJZy6DYRhrHnHVEojGkJ5S0UfN5xwV/LJEiPqkikSEjnCQtvYwg+Es/3BKKKhH2T4ac1LWVziNdELHAW/xmQ3uumrbHBORANCOk6CwaMTnhtIVir1eVej5HzIJG8Mwlg7xB6CtHdoqJt0uRKqoxjNMNu9YFc+lRjqhPcB2EdmK42zeCvxuxTb3AbcB9wNvBH6yFONBc0giiS9dTRV6EF9Xj2fsJmkSNoZhLAs1pYoK+fKoKe2UeamubD+Ib3ARiVItRMOckDvGcwfwfZxuzC+q6gER+TiwV1XvA74AfEVEDgFncRzV0pFIwCkfvvWDdF5yBXckrICVYRitifj8kGyDZNs8qSJfOkVnapI7ilJFaWdoYBWwqks5aD4PIiZhYxjGqkILBVBd0Li0lXJoApY4YBjGamQ1NaxXz5EYhmEYKw5zQoZhGEbTMCdkGIZhNI0Vl5ggIqeBIxf58R4q1BhaBLNrYZhdC6dVbTO7FsZi7Nqsqi2neLrinNBiEJG9rZgdYnYtDLNr4bSqbWbXwmhVuxaDdccZhmEYTcOckGEYhtE01poTurPZBtTA7FoYZtfCaVXbzK6F0ap2XTRrakzIMAzDaC3WWiRkGIZhtBDmhAzDMIymsaqdkIh8UkQOisg+EfmWiHTU2O5mEXlaRA6JyEeWwa43icgBESmISM10SxE5LCL7ReQxEalPtXV57Fru89UlIj8UkWfdv501tsu75+oxEbmvgfac9/hFJCwiX3fff1BEtjTKlgXa9S4ROe05R7+/THZ9UUROicgTNd4XEfkn1+59InJdi9j1ShEZ95yvv14GmzaKyE9F5En3XvzjKts05Xw1DFVdtS/g1UDAXf4E8Ikq2/iB54BtQAh4HNjZYLt2AJcDPwOGzrPdYaBnGc/XBe1q0vn6R+Aj7vJHqv2O7nupZThHFzx+4I+Af3GX3wp8vUXsehfwmeW6njzf++vAdcATNd6/BfgeIMD1wIMtYtcrgf9e5nPVD1znLieBZ6r8jk05X416repISFV/oKpzdXQfwKnuWsku4JCqPq+qs8DXgFsbbNdTqvp0I7/jYqjTrmU/X+7+73aX7wZe3+DvOx/1HL/X3nuBm6TxJTCb8bvUhar+HKdeWC1uBb6sDg8AHSLS3wJ2LTuqOqyqj7jLk8BTQGUZ6Kacr0axqp1QBe/GaT1UMggc9fx/jPk/erNQ4Aci8rCIvKfZxrg043z1qeqwu3wS6KuxXURE9orIAyLSKEdVz/EXt3EbQeNAd4PsWYhdAG9wu3DuFZGNDbapXlr5HnyZiDwuIt8TkSuX84vdbtxfAx6seKuVz9eCWfH1hETkR8D6Km99VFW/427zUSAHfLWV7KqDG1X1uIisA34oIgfd1luz7VpyzmeX9x9VVRGpNa9gs3u+tgE/EZH9qvrcUtu6gvkv4B5VzYjIH+JEa69qsk2tzCM411RKRG4Bvg1sX44vFpEE8J/An6jqxHJ8Z7NY8U5IVXef730ReRfwWuAmdTtUKzgOeFuEG9x1DbWrzn0cd/+eEpFv4XS5LMoJLYFdy36+RGRERPpVddjtdjhVYx9z5+t5EfkZTityqZ1QPcc/t80xEQkA7cCZJbZjwXapqteGu3DG2lqBhlxTi8X78FfV74rI50SkR1UbKmwqIkEcB/RVVf1mlU1a8nxdLKu6O05Ebgb+DHidqk7V2GwPsF1EtopICGcguWGZVfUiInERSc4t4yRZVM3iWWaacb7uA25zl28D5kVsItIpImF3uQd4OfBkA2yp5/i99r4R+EmNBtCy2lUxbvA6nPGGVuA+4J1u1tf1wLin+7VpiMj6ubE8EdmF87xsaGPC/b4vAE+p6qdqbNaS5+uiaXZmRCNfwCGcvtPH3NdcxtIA8F3PdrfgZKE8h9Mt1Wi7fhunHzcDjADfr7QLJ8vpcfd1oFXsatL56gZ+DDwL/AjoctcPAXe5yzcA+93ztR+4vYH2zDt+4OM4jR2ACPAN9/p7CNjW6HNUp11/715LjwM/Ba5YJrvuAYaBrHt93Q68F3iv+74An3Xt3s95MkaX2a47POfrAeCGZbDpRpyx4H2e59YtrXC+GvUy2R7DMAyjaazq7jjDMAyjtTEnZBiGYTQNc0KGYRhG0zAnZBiGYTQNc0KGYRhG0zAnZBiGYTQNc0KGcQHcSYF2rxhGA7AbyzCqICJb3No8X8ZRqviCK456QEQ+5tnusIh8TEQeEaf20xXu+l5xah8dEJG7ROSIq+SAiLxDRB5ya9T8q4j4m3OUhtF8zAkZRm22A59T1SuBD6nqEPAi4DdE5EWe7UZV9Trg88CH3XV/gyPXcyVOOYdNACKyA3gL8HJVvRbIA29flqMxjBZkxQuYGkYDOaJOvRaAN7vlNAI4hcd24kirAMyJTD4M/I67fCOODBKq+j8ics5dfxPwYmCPK0sWpYYgq2GsBcwJGUZt0gAishUnwnmJqp4TkX/H0YebI+P+zXPhe0qAu1X1L5bYVsNYkVh3nGFcmDYchzQuIn3Aa+r4zC+ANwOIyKuBTnf9j4E3ujWiEJEuEdm89CYbxsrAIiHDuACq+riIPAocxFFl/0UdH/sYcI+I/B5wP05F2ElVHRWRv8SpmOvDUXB+P3CkMdYbRmtjKtqG0QDc2kZ5Vc2JyMuAz7uJCIZheLBIyDAawybgP9xoZxb4gybbYxgtiUVChmEYRtOwxATDMAyjaZgTMgzDMJqGOSHDMAyjaZgTMgzDMJqGOSHDMAyjafw/VVdDvdYilyoAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Comparing the two results from the graph, we can see that the two are almost the same as seen also from their percent error. "
      ],
      "metadata": {
        "id": "8Cm9cUHD_6VX"
      }
    }
  ]
}