Suico_JohnSamuel_module[2].ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "G0BcK0H_G1Iu"
   },
   "source": [
    "# **AP155 Lab Assignment** \n",
    "## Module 2: 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 October 8 (Thu)**\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PeHZ9kP8JaVR"
   },
   "source": [
    "### Student Information \n",
    "\n",
    "_Full Name (Last Name, First Name)_: Suico, John Samuel O. <br>\n",
    "_Student No._: 2018-05897 <br>\n",
    "_Section_: WFU-2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "Hc183co-J4eI"
   },
   "source": [
    "### Submission Information\n",
    "\n",
    "_Date and Time Submitted (most recent upload)_: October 9, 2020 9:00 PM\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:** John Samuel O. Suico"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "X1UGsRpW-eQs"
   },
   "source": [
    "### Grading Information (c/o Lab Instructor)\n",
    "\n",
    "TOTAL SCORE: **XX**/100\n",
    "\n",
    "Score breakdown:\n",
    "* Problem 1 - XX/20\n",
    "* Problem 2 - XX/20\n",
    "* Problem 3 - XX/10\n",
    "* Problem 4 - XX/30\n",
    "* Problem 5 - XX/20\n",
    "\n",
    "_Date and Time Scored (MM/DD/YYYY HH:MM AM/PM):_ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "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": "markdown",
   "metadata": {},
   "source": [
    "1. Here I've defined a function `intTR` which calculates the integral of a mathematical function using the trapezoidal rule. The trapezoidal rule calculates the integral by approximating the area under the curve using trapezoids.\n",
    "   \n",
    "   I first create the function `function1` to represent the integrand. This returns the value of the mathematical function above for an argument x.\n",
    "   \n",
    "   I then created the function `intTR` with the parameters `lowerB`, `upperB`, `N`, `function`, and `adaptive` (i.e. with a default value of `False`). The `lowerB` and `upperB` parameters stand in for the lower and upper bounds of integration, respectively. The `N` parameter represents the number of slices for a given range `lowerB` to `upperB`. The `function` parameter stands in for the integrand. The `adaptive` parameter (i.e. which is not used in this example) is a boolean parameter used to increase the reusability of the `intTR` function for the adaptive trapezoidal rule.\n",
    "   \n",
    "   The integral is calculated using the trapezoidal rule via the equation\n",
    "   $$I(a,b) = h \\left[ \\frac{1}{2}f(a) + \\frac{1}{2}f(b) + \\sum_{k=\\text{odd} \\\\ 1...N-2} f(a+kh) + \\sum_{k=\\text{even} \\\\ 2...N-2} f(a+kh)\\right]$$\n",
    "    \n",
    "    , where $h = \\frac{b-a}{N}$. In the function `intTR`, $a$ and $b$ are replaced by `lowerB` and `upperB` (respectively), and $h$ is replaced by `width`. The summation terms with odd and even $k$ are calculated using for loops and assigned to the variables `oddSum` and `evenSum`, respectively. <sub>[1]</sub>\n",
    "    \n",
    "    The terms are then added up to the variable `riemannSum` to calculate the full value of the integral. The variable `riemannSum` is then returned.\n",
    "    \n",
    "    I then proceed to calculate the integral using the function `intTR` for 10<sub>[2]</sub>, 100 and 1,000 slices, and print out my results.\n",
    "    \n",
    "    <sub> [1] It is useful to calculate the trapezoidal rule this way so that the `oddSum` can be returned when calculating using the adaptive trapezoidal rule.</sub>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "hMqWOsbIMoUh"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral using Trapezoidal Rule for 10 slices: 1471.4375\n",
      "Integral using Trapezoidal Rule for 100 slices: 1450.2145812500003\n",
      "Integral using Trapezoidal Rule for 1000 slices: 1450.0021458331257\n"
     ]
    }
   ],
   "source": [
    "from numpy import abs\n",
    "\n",
    "def function1(x):\n",
    "    y = 2*x**4 + 3*x**2 + 4*x + 5\n",
    "    return y\n",
    "\n",
    "def intTR(lowerB,upperB,N,function,adaptive = False):\n",
    "    width = (upperB - lowerB)/N\n",
    "    riemannSum, oddSum, evenSum = 0,0,0\n",
    "    \n",
    "    for k in range(1,N,2):\n",
    "        oddSum += function(lowerB + k*width)\n",
    "    if(adaptive == True):\n",
    "        return width*oddSum\n",
    "    else:\n",
    "        for k in range(2,N,2):\n",
    "            evenSum += function(lowerB + k*width)\n",
    "        riemannSum = width*(0.5*(function(lowerB) + function(upperB)) + oddSum + evenSum)\n",
    "        return riemannSum\n",
    "\n",
    "    \n",
    "print(\"Integral using Trapezoidal Rule for 10 slices:\", intTR(0,5,10,function1))\n",
    "print(\"Integral using Trapezoidal Rule for 100 slices:\", intTR(0,5,100,function1))\n",
    "print(\"Integral using Trapezoidal Rule for 1000 slices:\", intTR(0,5,1000,function1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Here I've defined the function `percentError` which returns the percentage error given the arguments `Exp` for the experimental value and `Theo` for the actual value. \n",
    "\n",
    "    I first calculate the value of the integral analyticaly by hand to be $1450$. I then check the percentage error of the integral calculated using the trapezoidal rule via my `percentError` function for 10, 100 and 1,000 slices and print out my results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage error of Trapezoidal Rule for 10 slices: 1.478448275862069\n",
      "Percentage error of Trapezoidal Rule for 100 slices: 0.014798706896569988\n",
      "Percentage error of Trapezoidal Rule for 1000 slices: 0.0001479884914292127\n"
     ]
    }
   ],
   "source": [
    "def percentError(Exp,Theo):\n",
    "    return abs((Exp-Theo)/Theo)*100\n",
    "\n",
    "#calculate theoretical value by hand\n",
    "print(\"Percentage error of Trapezoidal Rule for 10 slices:\", percentError(intTR(0,5,10,function1), 1450))\n",
    "print(\"Percentage error of Trapezoidal Rule for 100 slices:\", percentError(intTR(0,5,100,function1), 1450))\n",
    "print(\"Percentage error of Trapezoidal Rule for 1000 slices:\", percentError(intTR(0,5,1000,function1), 1450))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Here I've defined a function `intSR` which calculates the integral of a mathematical function using Simpson's rule. Simpson's rule calculates the integral by approximating the curve using quadratics and then  calculating the area under the quadratics.\n",
    "   \n",
    "   The `intSR` function has the parameters `lowerB`, `upperB`, `N`, and `function`. The `lowerB` and `upperB` parameters stand in for the lower and upper bounds of integration, respectively. The `N` parameter represents the number of slices for a given range `lowerB` to `upperB`. The `function` parameter stands in for the integrand.\n",
    "   \n",
    "   The integral is calculated using Simpson's rule via the equation\n",
    "   $$I(a,b) \\approx \\frac{1}{3}h \\left[ f(a) + f(b) + 4\\sum_{k=odd \\\\ 1...N-2} f(a+kh) + 2\\sum_{k=even \\\\ 2...N-2} f(a+kh)\\right]$$\n",
    "    \n",
    "    , where $h = \\frac{b-a}{N}$. In the function `intSR`, $a$ and $b$ are replaced by `lowerB` and `upperB` (respectively), and $h$ is replaced by `width`. The summation terms with odd and even $k$ are calculated using for loops and assigned to the variables `oddSum` and `evenSum`, respectively.\n",
    "    \n",
    "    The terms are then added up to the variable `riemannSum` to calculate the full value of the integral using the trapezoidal rule. The variable `riemannSum` is then returned.\n",
    "    \n",
    "    I then proceed to calculate the integral using the function `intSR` for 10, 100 and 1,000 slices, and print out my results.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral using Simpson's Rule for 10 slices: 1450.0833333333333\n",
      "Integral using Simpson's Rule for 100 slices: 1450.0000083333337\n",
      "Integral using Simpson's Rule for 1000 slices: 1450.0000000008342\n"
     ]
    }
   ],
   "source": [
    "def intSR(lowerB,upperB,N,function):\n",
    "    width = (upperB-lowerB)/N\n",
    "    riemannSum, oddSum, evenSum = 0,0,0\n",
    "    for k in range(1,N,2):\n",
    "        oddSum += function(lowerB + k*width)\n",
    "    for k in range(2,N-1,2):\n",
    "        evenSum += function(lowerB + k*width)\n",
    "    riemannSum = (3**-1)*width*( function(lowerB) + function(upperB) + 4*oddSum + 2*evenSum )\n",
    "    return riemannSum\n",
    "    \n",
    "print(\"Integral using Simpson's Rule for 10 slices:\", intSR(0,5,10,function1))\n",
    "print(\"Integral using Simpson's Rule for 100 slices:\", intSR(0,5,100,function1))\n",
    "print(\"Integral using Simpson's Rule for 1000 slices:\", intSR(0,5,1000,function1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. Here, I check the percentage error of the integral calculated using Simpson's rule via my `percentError` function for 10, 100 and 1,000 slices and print out my results. (Actual value = $1450$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Percentage error of Simpson's Rule for 10 slices: 0.005747126436776382\n",
      "Percentage error of Simpson's Rule for 100 slices: 5.747126657485657e-07\n",
      "Percentage error of Simpson's Rule for 1000 slices: 5.753338035870472e-11\n"
     ]
    }
   ],
   "source": [
    "def percentError(Exp,Theo):\n",
    "    return abs((Exp-Theo)/Theo)*100\n",
    "\n",
    "#calculate theoretical value by hand\n",
    "print(\"Percentage error of Simpson's Rule for 10 slices:\", percentError(intSR(0,5,10,function1), 1450))\n",
    "print(\"Percentage error of Simpson's Rule for 100 slices:\", percentError(intSR(0,5,100,function1), 1450))\n",
    "print(\"Percentage error of Simpson's Rule for 1000 slices:\", percentError(intSR(0,5,1000,function1), 1450))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "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": "markdown",
   "metadata": {},
   "source": [
    "1. Here, I've defined a class called `aintTR` with class methods used to calculate the integral using the adaptive trapezoidal rule. The class has the instance variables `lowerB`, `upperB`, `function`, `accuracy`, `I_list`, and `error_list`. The variables `lowerB` and `upperB` represent the lower and upper bounds of integration, respectively. The `function` variable stands in for the integrand. The `accuracy` variable stands in for the accuracy we are aiming for when calculating using the adaptive trapezoidal rule. The `I_list` and `error_list` variables keep track of the integral estimates and error values, which will be useful for graphing Error vs N slices.\n",
    "\n",
    "    The adaptive trapezoidal rule calculates the $i$th estimate of the integral via the equation $$ I_i = \\frac{1}{2}I_{i-1} + h_i \\sum_{k=\\text{odd} \\\\ 1...N-2} f(a+kh_i) $$\n",
    "    , with the first integral estimate $I_1$ being calculated by the normal trapezoidal rule (which has already been coded in Problem #1.1). Each subsequent calculation $I_i$ doubles the number of slices `N` and makes use of the fact that the previous integral estimates are already embedded in the integral estimate of the doubled `N`, and we only have to account for the summation term with odd $k$. This is done by calling the same function `intTR` in Problem #1.1 but with the `adaptive` argument set to`True`, which then returns the second term in the RHS of the equation above (of course, we are also doubling the value of `N` with each subsequent calculation of $I_i$). Each estimate of $I_i$ also has an associated error $\\epsilon_i$ inherent within the trapezoidal rule, and this error is calculated using a practical equation with respect to the previous estimate $I_{i-1}$ $$ \\epsilon_i = \\frac{1}{3}(I_i - I_{i-1})$$\n",
    "    As such, we will only be calculating errors associated with the 2nd and above integral estimates (i.e. $i$>2).\n",
    "    \n",
    "    I've also defined a new function `function2` to return the value of the new integrand given above per argument `x`. This is the new function we pass on to the `aintTR`'s `function` parameter.\n",
    "    \n",
    "    <h3>aintTR Class Methods</h3>\n",
    "    \n",
    "    I'll first be explaining the `approximate` function which calculates the integral to a desired accuracy set by the `accuracy` instance variable. Firstly, we reset the instance variables `I_list` and `error_list` to make sure that we are getting freshly calculated integral estimates and error values, respectively. Afterwards, we calculate the first integral estimate `I_current` using the function defined in Problem #1.1 `intTR` and append the value of `I_current` to `I_list`. We then check for the error of the 2nd integral estimate <b>without</b> explicitly doubling our function variable `N`. If the error for the `2nd integral` is less than the target accuracy value, then `I_current` is probably a good estimate for your integral.<sub>[1]</sub> \n",
    "    \n",
    "    The error is checked using the `erroraTR` function, which takes the arguments `I_prev` and `I_current` as the previous and current estimates of the integral, and then calculates the value using the equation for $\\epsilon$ given above.\n",
    "    \n",
    "    If our value `I_1` is not good enough, then the program goes through the `while` loops with the condition that `error` is greater than the instance variable `accuracy`. This will then proceed to calculate and store the next integral estimates and errors<sub>[2]</sub> using the equations above, until `error` which is the error associated for the current estimate $I_i$ is less than the target accuracy given by the instance variable `accuracy`.\n",
    "    \n",
    "    Lastly, the `output` function just prints out each calculated integral and their associated error (only $I_1$ has no associated error since it cannot be calculated using the equation above). This is done by using a string variable `myString` and looping  through the lists `I_list` and `error_list`.\n",
    "    \n",
    "    <sub>[1] This was really only done to account for line polynomials. In such a case, the trapezoidal function will be able to perfectly/almost perfectly calculate the integral at the first try using the normal trapezoidal rule.</sub>\n",
    "    \n",
    "    <sub>[2] Each integral estimate is saved to `I_list` and each error is saved to `error_list`.</sub>\n",
    "    \n",
    "    I then create an instance of the class `aintTR` named `myaintTR` with the arguments `0`, `1`, `function2` and `10**-6`, and then call the functions `approximate` and `output`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "GABjwlXCR97i"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adaptive Trapezoidal Rule for 1 slice(s): 0.147979484546652\n",
      "Adaptive Trapezoidal Rule for 2 slice(s): 0.3252319078064746 ; Error = 0.05908414108660753\n",
      "Adaptive Trapezoidal Rule for 4 slice(s): 0.5122828507233315 ; Error = 0.06235031430561895\n",
      "Adaptive Trapezoidal Rule for 8 slice(s): 0.4029974484782483 ; Error = 0.03642846741502771\n",
      "Adaptive Trapezoidal Rule for 16 slice(s): 0.43010336929474696 ; Error = 0.009035306938832883\n",
      "Adaptive Trapezoidal Rule for 32 slice(s): 0.4484146657874699 ; Error = 0.0061037654975743165\n",
      "Adaptive Trapezoidal Rule for 64 slice(s): 0.4539129312153758 ; Error = 0.001832755142635293\n",
      "Adaptive Trapezoidal Rule for 128 slice(s): 0.45534850437280205 ; Error = 0.000478524385808754\n",
      "Adaptive Trapezoidal Rule for 256 slice(s): 0.455711266453241 ; Error = 0.00012092069347964991\n",
      "Adaptive Trapezoidal Rule for 512 slice(s): 0.45580219965166413 ; Error = 3.0311066141042176e-05\n",
      "Adaptive Trapezoidal Rule for 1024 slice(s): 0.45582494813241997 ; Error = 7.582826918613635e-06\n",
      "Adaptive Trapezoidal Rule for 2048 slice(s): 0.4558306362016466 ; Error = 1.8960230755427077e-06\n",
      "Adaptive Trapezoidal Rule for 4096 slice(s): 0.455832058278271 ; Error = 4.7402554146739695e-07\n"
     ]
    }
   ],
   "source": [
    "from numpy import sin,sqrt, inf\n",
    "\n",
    "def function2(x):\n",
    "    y = sin(sqrt(100*x))**2\n",
    "    return y\n",
    "\n",
    "class aintTR:\n",
    "    def __init__(self, lowerB, upperB, function, accuracy):\n",
    "        self.lowerB = lowerB\n",
    "        self.upperB = upperB\n",
    "        self.function = function\n",
    "        self.accuracy = accuracy\n",
    "        self.I_list = []\n",
    "        self.error_list = []\n",
    "    \n",
    "    def erroraTR(self, I_prev, I_current):\n",
    "        return abs((3**-1)*(I_current - I_prev))\n",
    "    \n",
    "    def approximate(self):\n",
    "        #reset I_list and error_list\n",
    "        self.I_list = []\n",
    "        self.error_list = []\n",
    "        \n",
    "        #first integral\n",
    "        N = 1\n",
    "        I_current = intTR(self.lowerB,self.upperB,N,self.function)\n",
    "        self.I_list.append(I_current)\n",
    "        \n",
    "        #error check for second integral\n",
    "        error = self.erroraTR(I_current, intTR(self.lowerB,self.upperB,N*2,self.function))\n",
    "        \n",
    "        while(error>self.accuracy):\n",
    "            N *= 2\n",
    "            I_prev = I_current\n",
    "            I_current = 0.5*I_current+intTR(self.lowerB,self.upperB,N,self.function,True)\n",
    "            self.I_list.append(I_current)\n",
    "            error = self.erroraTR(I_prev,I_current)\n",
    "            self.error_list.append(error)\n",
    "    \n",
    "    def output(self):\n",
    "        looplimit = len(self.I_list)\n",
    "        N = 1\n",
    "        for i in range(0,looplimit):\n",
    "            myString = \"Adaptive Trapezoidal Rule for \" + str(N) + \" slice(s): \" + str(self.I_list[i])\n",
    "            if i > 0:\n",
    "                myString = myString + \" ; Error = \"  + str(self.error_list[i-1])\n",
    "            N *= 2\n",
    "            print(myString)\n",
    "\n",
    "myaintTR = aintTR(0,1,function2,10**-6)\n",
    "myaintTR.approximate()\n",
    "myaintTR.output()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Here I've defined a class called `intRI` with class methods used to calculate the integral of a function using the Romberg integration method. The class has the instance variables `lowerB`, `upperB`, `function`, and `accuracy`. The variables `lowerB` and `upperB` represent the lower and upper bounds of integration, respectively. The `function` variable stands in for the integrand. The `accuracy` variable stands in for the accuracy we are aiming for when calculating using the Romberg integration method. The `I_list` and `error_list` variables keep track of the integral estimates and error values, which will be useful for graphing Error vs N slices.\n",
    "\n",
    "    The Romberg integration method is summarized in this diagram (taken from Newman):\n",
    "    \n",
    "    <img src = \"RImethod.png\" width = \"300\">\n",
    "    \n",
    "    The first entry of each row $\\left(R_{i,1}\\right)$ is calculated using the normal trapezoidal rule (which has already been coded in Problem #1.1). Subsequent integral estimates allow us to calculate values which have been normally been \"shaved off\" by the trapezoidal function, and therefore allow us to get more accurate estimates of the integral. The $i$ in each row's $I_i$ stands for the doubling of `N` slices starting from `N` = 1, and so each row can be thought of having $2^{i-1}$ slices. The next entries $R_{i,m+1}$ in each row are then calculated using the equation\n",
    "    $$R_{i,m+1} = R_{i,m} + \\frac{1}{4^m - 1}\\left( R_{i,m} - R_{i-1,m}\\right)$$\n",
    "    and the error associated with $R_{i,m}$ is\n",
    "    $$\\epsilon = \\frac{1}{4^m - 1}\\left( R_{i,m} - R_{i-1,m}\\right)$$\n",
    "    This also means that each row can only be made more accurate up to $i-1$ times, since it is dependent on the previous row's values. Also, the Romberg integration method does not give the associated error with the last term in each row. Note that we can also rewrite $R_{i,m+1}$ as\n",
    "    $$R_{i,m+1} = R_{i,m} + \\epsilon$$\n",
    "    \n",
    "    <h3>intRI Class Methods</h3>\n",
    "    \n",
    "    I'll first be explaining the `approximate` function which calculates the integral to a desired accuracy set by the `accuracy` instance variable. Firstly, the instance variables `I_list` and `error_list` are reset to ensure that we are storing freshly calculated integral estimates and error values, respectively. Then, we calculate the first integral estimate `I_current` for the first row $I_1$ using the `intTR` function defined in Problem #1.1 and append `I_current` to the list of integral estimates `I_list`. We then check for the error of the 2nd integral estimate <b>without</b> explicitly doubling our function variable `N`. If the error for the `2nd integral` is less than the target accuracy value, then `I_current` is probably a good estimate for your integral.<sub>[1]</sub> \n",
    "    \n",
    "    The error is checked using the `errorRI` function, which takes the arguments `R_a`, `R_b`, and `m` where `R_a` and `R_b` are integral estimates with the same $m$ for which we can calculate the error value using the equation for $\\epsilon$. In the diagram, `R_a` and `R_b` would be in the same column, with `R_b` being the integral estimate with the larger `N`.\n",
    "    \n",
    "    If the absolute value of the error of the 2nd integral estimate is larger than the desired accuracy set by the `accuracy` instance variable, then we will continually loop and calculate further integral estimates and errors<sub>[2]</sub> using the equations above until the desired accuracy is reached.\n",
    "    \n",
    "    Lastly, the `output` function just prints out each calculated integral and their associated error (the last integral estimate per row has no associated error value). This is done by using a string variable `myString` and looping  through the lists `I_list` and `error_list`.\n",
    "    \n",
    "    <sub>[1] Once again, this was really only done to account for line polynomials.</sub>\n",
    "    \n",
    "    <sub>[2] The integral estimates and error values are stored into the instance variables `I_list` and `error_list`</sub>\n",
    "    \n",
    "    I then create an instance of the class `intRI` named `myintRI` with the arguments `0`, `1`, `function2` and `10**-6`, and then call the functions `approximate` and `output`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I_1: Romberg Integration for 1 slice(s): 0.147979484546652\n",
      "I_2: Romberg Integration for 2 slice(s): 0.3252319078064746 ; Error = 0.05908414108660753\n",
      "I_2: Romberg Integration for 2 slice(s): 0.38431604889308213\n",
      "I_3: Romberg Integration for 4 slice(s): 0.5122828507233315 ; Error = 0.06235031430561895\n",
      "I_3: Romberg Integration for 4 slice(s): 0.5746331650289505 ; Error = 0.012687807742391222\n",
      "I_3: Romberg Integration for 4 slice(s): 0.5873209727713417\n",
      "I_4: Romberg Integration for 8 slice(s): 0.40299744847824825 ; Error = -0.036428467415027734\n",
      "I_4: Romberg Integration for 8 slice(s): 0.3665689810632205 ; Error = -0.01387094559771533\n",
      "I_4: Romberg Integration for 8 slice(s): 0.35269803546550516 ; Error = -0.0037241736080291504\n",
      "I_4: Romberg Integration for 8 slice(s): 0.34897386185747603\n",
      "I_5: Romberg Integration for 16 slice(s): 0.4301033692947469 ; Error = 0.009035306938832883\n",
      "I_5: Romberg Integration for 16 slice(s): 0.43913867623357977 ; Error = 0.004837979678023951\n",
      "I_5: Romberg Integration for 16 slice(s): 0.44397665591160373 ; Error = 0.001448866991207914\n",
      "I_5: Romberg Integration for 16 slice(s): 0.4454255229028116 ; Error = 0.0003782418080209239\n",
      "I_5: Romberg Integration for 16 slice(s): 0.44580376471083255\n",
      "I_6: Romberg Integration for 32 slice(s): 0.4484146657874699 ; Error = 0.006103765497574335\n",
      "I_6: Romberg Integration for 32 slice(s): 0.45451843128504427 ; Error = 0.0010253170034309662\n",
      "I_6: Romberg Integration for 32 slice(s): 0.4555437482884752 ; Error = 0.0001836046409027219\n",
      "I_6: Romberg Integration for 32 slice(s): 0.45572735292937794 ; Error = 4.039933343751497e-05\n",
      "I_6: Romberg Integration for 32 slice(s): 0.45576775226281546 ; Error = 9.739968281508223e-06\n",
      "I_6: Romberg Integration for 32 slice(s): 0.455777492231097\n",
      "I_7: Romberg Integration for 64 slice(s): 0.4539129312153758 ; Error = 0.001832755142635293\n",
      "I_7: Romberg Integration for 64 slice(s): 0.4557456863580111 ; Error = 8.181700486445582e-05\n",
      "I_7: Romberg Integration for 64 slice(s): 0.4558275033628756 ; Error = 4.504048800005895e-06\n",
      "I_7: Romberg Integration for 64 slice(s): 0.45583200741167557 ; Error = 4.1040973450050226e-07\n"
     ]
    }
   ],
   "source": [
    "class intRI:\n",
    "    def __init__(self, lowerB, upperB, function, accuracy):\n",
    "        self.lowerB = lowerB\n",
    "        self.upperB = upperB\n",
    "        self.function = function\n",
    "        self.accuracy = accuracy\n",
    "        self.I_list = []\n",
    "        self.error_list = []\n",
    "        \n",
    "    def errorRI(self,R_a,R_b,m):\n",
    "        return ((4**m - 1)**-1)*(R_b-R_a)\n",
    "\n",
    "    def approximate(self):\n",
    "        self.I_list = []\n",
    "        self.error_list = []\n",
    "        #first integral: R_1,1\n",
    "        N = 1\n",
    "        I_current = intTR(self.lowerB, self.upperB, N, self.function)\n",
    "        newList = [I_current]\n",
    "        self.I_list.append(newList)\n",
    "        \n",
    "        I_next = intTR(self.lowerB, self.upperB, N*2, self.function)\n",
    "        \n",
    "        #error check for R_2,1\n",
    "        error = self.errorRI(I_current,I_next,1)\n",
    "              \n",
    "        while(abs(error) > self.accuracy):\n",
    "            N *= 2\n",
    "            #initializing R_i,1\n",
    "            I_first = intTR(self.lowerB,self. upperB, N, self.function)\n",
    "            newList = [I_first]\n",
    "            looplimit = len(self.I_list)\n",
    "            \n",
    "            for m in range(0,looplimit):\n",
    "                error = self.errorRI(self.I_list[looplimit-1][m],newList[m],m+1)\n",
    "                self.error_list.append(error)\n",
    "                if (abs(error) > self.accuracy):\n",
    "                    I_next = newList[m] + error\n",
    "                    newList.append(I_next)\n",
    "                else:\n",
    "                    break\n",
    "            self.I_list.append(newList)\n",
    "\n",
    "    def output(self):\n",
    "        looplimit = len(self.I_list)\n",
    "        erroriterator = 0\n",
    "        N = 1\n",
    "        for i in range (0,looplimit):\n",
    "            newlooplimit = len(self.I_list[i])\n",
    "            for j in range(0,newlooplimit):\n",
    "                myString = \"I_\" + str(i+1) + \": Romberg Integration for \" + str(N) + \" slice(s): \" + str(self.I_list[i][j])\n",
    "                if (j != i):\n",
    "                    myString = myString + \" ; Error = \" + str(self.error_list[erroriterator])\n",
    "                    erroriterator += 1\n",
    "                print(myString)\n",
    "            N *= 2\n",
    "myintRI = intRI(0,1,function2,10**-6)              \n",
    "myintRI.approximate()\n",
    "myintRI.output()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the results from both calculations by creating a plot of error estimates vs.  𝑁 . Annotate your chart completely (include proper axis labels and legend). Describe how the trends for the two techniques compare\n",
    "\n",
    "3. Here, I've graphed the truncation errors vs N of both the adaptive trapezoidal rule and the Romberg Integration Method in calculating the integral of $$\\int_0^1 \\sin^2 \\sqrt{100x} \\mbox{d}x $$ to an accuracy of $10^{-6}$. I first calculated appropriate $N$'s as x values for the errors calculated in Problem #2.1 and Problem #2.2. I then use the class instance `myaintTR`'s `error_list` as y values for the adaptive trapezoidal rule. For the Romberg Integration Method, I've only included the last truncation errors per $N$ for its y values. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEaCAYAAADg2nttAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABI6klEQVR4nO3dd3gU5fbA8e8hlFCCwCVeFZCgohIgJIAUqaJIkS4KiEpRVK6KFcWG2L0iKijXckVAuUYQFLH97ChNBRQVUQQVBBsx0osEOL8/3knchJRN2Ux293yeZ55kZ2dnzs6Wd+ct5xVVxRhjjMlPOb8DMMYYU/ZZYWGMMaZAVlgYY4wpkBUWxhhjCmSFhTHGmAJZYWGMMaZAVlj4QERmiMjdJbzP4SKyuCT3GQoiMkFEZvlw3HYisk5EdolIv9I+fkkTkQ4istbvOKJJzs+Y9146zs+YSpMVFiEkIgtFZKuIVPI7lkDFKVhEpI2I7BaRarnc97mIXFH8CEPiTuAxVa2mqvMD7/A+9JnLIRHZG3B7qD/hZiciKiInZN5W1UWqelIIjpPgHWtXjmVQSR8r1LwfZSoirQLWnSAieQ4uE5H2IrJURLaLyJ8iskRETsltW++99EMoYi+LrLAIERFJADoACvTxN5qSo6ofA5uBgYHrRaQJkAik+hFXEOoDX+d2h/ehr6aq1YCfgN4B6/6XuZ2IlC+lWMuCGoHnRVVn57aRiMTkuF2oc1QK5/RPIKireBGpDrwGPArUAuoAdwB/hSy6MGKFRehcCHwMzACG5XJ/bRF5R0R2isiHIlIfQJyHRWSLiOwQka+8L2JE5AgReVZE0kRko4jcKiKHvYYBvw7LB6xbKCIXi0gj4AmgrfeLcZt3fyUReVBEfhKR30XkCRGpnMdzm+k9v5zP9w1VTReRySKyyYt/pYh0yG0nItJZRDbnWLdBRM7w/i8nIuNE5HsRSReROSJSK4+YEJFRIrLe+0W4QESO8dZ/DxwHvOo956Cu9DLjE5EbReQ3YHpuV2WBv/q9X7NTReR177X9RESOD9i2sfe6/+md55u99a1EZJmIbBORX0XkMRGp6N33kffwLzJ/5ec8dyLSyHuNt4nI1yLSJ+C+fGMqDG9fj4vIGyKyGzjNe81uFJEvgd0iUl5E+nhxbPPiahSwj8O2z3GMx0XkwRzrXhGRa73/bxSRn73nslZETs8n5JlAkoh0CuLpnQigqqmqelBV96rq26r6ZR7nIvB1rywik7zP5XYRWZz5+RF3Nb7UOxdfiEjngH0MF5EfvOfyo5SRK9lcqaotIViA9cC/gBZABvDPgPtmADuBjkAlYDKw2LuvG7ASqAEI0Ag42rvvWeAVIA5IAL4DLvLuGx6wjwTcFU35gGMuBC7OuW3A/Q8DC3C/qOKAV4H78nhu9YADQD3vdjnc1UY/7/b5wD+A8sB1wG9ArHffBGCW939nYHOOfW8AzvD+vwpX4Nb1ztOTQGoeMXUB/gCae9s+CnyU234LeN0Cj9/Ze57/9vZZOY9zp8AJAa9tOtDKe/7/A17w7osDfvXOSax3u7V3XwugjfeYBOAb4OrcjpHz3AEVcO+3m4GK3rnYCZxUUEy5PP/D3js57p8BbAfaea97rHfOVnnvi8q4L93dQFcvthu8+CoGnOOs7XM5RkdgEyDe7ZrAXuAY4CTvvmMC4j0+n1jvBsbw92fjBEDz2L66d55mAj2Amjnuz/ba53jdp+I+Y3WAGOBU3HumjrfPnt756urdjgeqAjsCXqejgcZ+f3fl+dnwO4BIXID2uAKitnf7W+CagPtnBH5YgWrAQe/D0wVXCLQBygVsEwPsBxID1l0KLPT+z3oj5/aBJ5/CAlco7Q780AFtgR/zeY7vAjd7/3cF0oAKeWy7FWjm/T+B4AuLb4DTA+472juvh32RAdOAB3Kc0wwgIed+C3jtAo/f2TvnsQH3Zzt33rqchcXTAff1BL71/h8CfB7ke+hq4OXcjpHz3OGqO3/L8X5JBSYUFFMux81872zLsTQK2NezuZyzkQG3bwPmBNwuB/wMdM5t+1xiEFx1YEfv9ijgfe//E4AtwBl5vd9yfM7uxn1p/4QrAPIsLLzHNPIetxn3Q2EB3g+9nK995mviPb+9eO/xHPu7EXgux7q3cLUNVb1zeza5FJplbbFqqNAYBrytqn94t5/n8KqoTZn/qOouXN3qMar6PvAY7pfKFhF5Slxdam3cr7SNAfvYiPvlUlzxQBVgpXepvA34P299XmYCF3j/X4Ar/DIAROR6EfnGuxzfBhzhxV9Y9YGXA2L6Bleo/jOXbY8h4Nx45zSd4p+fNFXdV8jH/Bbw/x5cwQXux8D3uT1ARE4UkddE5DcR2QHcS/Dn7Bhgk6oeCliX872RV0x5qa2qNQKWbwLu25TL9oHrcr4Wh7z76+SxfTbqvlFfwBWuAOfhroZQ1fW4gnQC7vPxQmZ1Yz77+wu4y1vyparfqOpwVa0LNPGeyyMFPKw27gort9e2PnBO5nvYex+3x9UW7AYGAZcBv3rVhCcXFKNfrLAoYV495blAJ++D/xtwDdBMRJoFbFov4DHVcNU/vwCo6hRVbYFrMD4RGIurYsnAvfkyHYv7xZbTbu9vlYB1RwX8n7M3yB+4X0aNA74cjlDX4JuXl4C6InIaMABXeCCufeIG7xzUVNUauGoLySPOrBjFNZYGFlCbgB45vrRiVTW35/wLAedGRKriqsJy27Ywcp6rnDEfRfA24dpOcvM47gq0oapWx1Up5XbOcvMLUE+yt1/l9d4oCbn1Jgpcl/O1ENz7/ec8ts9NKjBQXFtea2Be1gNVn1fV9t4xFFdNWJDpuKrdAUFsm3mcb3FXGU0K2PQPYB+QWzvQJtyVReB7uKqq3u8d4y1V7Yq7av4W+G+w8ZU2KyxKXj/cr99EINlbGgGLyN4o3FNcN72KuF88H6vqJhE5RURai0gF3BfTPuCQqh4E5gD3iEic9yG6FjhszIKqpuE+mOeLSIyIjCT7G/l33Bd9RW/7Q7g36cMiciSAiNQRkW55PUnvV9Fc3Idwo6qu8O6Kw12+pwHlRWQ8ri44N98BsSJylvd8b8VVGWR6wnu+9b2Y4kWkbx77SgVGiEiyuAbse4FPVHVDXs+hiL4AGnvHicX9wg3Wa8DRInK1uA4FcSLS2rsvDld/vcv7dTk6x2N/J++C5hPc1cINIlLBa0Dtjft17oc5wFkicrr3ul6H61G0NNgdqOrnuC/hp4G3VHUbgIicJCJdvNd4H+5HzqE8d/T3/g4At+OqhXIlIieLyHUiUte7XQ93dfNxAfs+BDwDPCQix3ifubZejLOA3iLSzVsfK65zQl0R+aeI9PV+2PwF7ArmufjFCouSNwyYrqo/qepvmQuuamloQM+P53Fv3j9xjZvne+ur4764t+Iu5dOBid59V+IKkB+Axd4+nskjjlG4K5J0oDHZP6jv47qR/iYimVVlN+IaIT/2qkHexTUm5mcm7tfdswHr3sJVYX3nxb+PPKocVHU7rhPA07jCbTeurjjTZFyd8dsishP3oW2dcz/evt7F1ZXPwzUiHw8MLiD+QlPV73BjNt4F1uFeh2AfuxPXvtMbVy20DjjNu/t6XHXLTtzrn7Or6gRgpleVcW6O/e739tkD9wX7H+BC75dxUW2T7OMsrg32gaq6Fvd+ftSLpzeuO/L+QsbwPK5t4vmAdZWA+739/gYcCdwU5P5Sce+NvOzEvb8+EdfT62NgNa6wK8j1wFfActxn+t+4NqRNQF/clWIa7rMwFvfdWw73g+8X7zGdOPxHQpmR2dvAGGOMyZNdWRhjjCmQFRbGGGMKZIWFMcaYAoVFYSEix4nINBGZ63csxhgTjUJeWIjIM+LyHK3Osb67uLwu60VkXH77UNUfVPWi0EZqjDEmL6WRRXMGrttoVvdKb/DVVFw3ws3AchFZgEtpcV+Ox49U1S2FOWDt2rU1ISGhGCEbY0z0Wbly5R+qmmvmhpAXFqr6kbh03YFaAevVywUvIi8AfVX1PqBXcY+ZkJDAihUrCt7QGGNMFhHZmNd9frVZ1CH7QK3N5JPDR0T+ISJPACkikusAHBG5RERWiMiKtLS0ko3WGGOiXFhM5qKq6bhkW/lt85SI/Ar0rlixYovSicwYY6KDX1cWPxOQSA83X0Gxk56p6quqeskRRxxR3F0ZY4wJ4NeVxXKgoYg0wBUSg3F5cYpFRHoDvU844YQCtzWmKDIyMti8eTP79hU2a7kxZUdsbCx169alQoUKQT8m5IWFiKTiJmqpLW4ayNtVdZqIXIFLOhcDPKOquc6PXBiq+irwasuWLUcVd1/G5Gbz5s3ExcWRkJCAy7xtTHhRVdLT09m8eTMNGjQI+nGl0RtqSB7r3wDeCPXxCy0tDTZsgIQEiM9v7p/iqV4ddu78+3ZcHOzYEbLDmRKyb98+KyhMWBMR/vGPf1DYjkBh0cAdrOJWQ42ITWXqXxexn4pUZD+XV5rG9H25lnUFUoWMDNi7F/bt+3vJvB1YUMDht03ZZQWFCXdFeQ+HRbqPYBWrgTstjXv/uppVNON/nMcjXM0//9rI2Mv3cOWVMGoUnH8+DBwIvXrBGWdAu3bQogU0bgzHHw916kCtWlClCsTEQKVKUKMGHHWUu1A5+WRISYG2bXMP4eWXYc+e4pwBEy3mz5+PiPDtt3lPWdG5c+cijzeaP38+a9asybo9fvx43n333SLtK9Pll19OcnIyiYmJVK5cmeTkZJKTk5k71/8sPgsWLOD+++/P9b5q1QqagTbvbWJiYkhOTqZJkyb07t2bbdu25bufhQsX0qtXsYeahYRdWWTasIFdxNGCz6jGbrrxFr/zT2KnKbFVoHJliI39e6lc2S01a/59O+f9+d3u2fPwEAYMcAVNjx7u/169XHWVMTmlpqbSvn17UlNTueOOO0p8//Pnz6dXr14kJiYCcOeddxZ7n1OnTgVgw4YN9OrVi1WrVmW7/8CBA5Qv789XUp8+fejTp0+J77dy5cpZz3PYsGFMnTqVW265pcSPUypUNeKWFi1aaKFt2aI7qKp/UkMPgR4C3U0l1S1bCr+vIMTFqbrKKrdUq6b67ruqo0erHnWUW1exomrPnqrTpqmmpYUkDFNIa9as8TsE3blzpx5zzDG6du1aPfHEE7PW79mzRwcNGqQnn3yy9uvXT1u1aqXLly9XVdXLLrtMW7RooYmJiTp+/Pisx9SvX1/Hjh2rTZo00VNOOUXXrVunS5Ys0Zo1a2pCQoI2a9ZM169fr8OGDdMXX3xR33zzTR04cGDW4z/44AM966yzVFX1rbfe0jZt2mhKSooOHDhQd+7cmWv8P/74ozZu3Djr8e3bt9fevXtrw4YNVVW1b9++2rx5c01MTNQnn3wy63FVq1bVq6++WhMTE7VLly66xftsrl+/Xrt166bNmzfX9u3b6zfffKOqqs2aNctaYmNjdeHChZqenq59+/bVpk2bauvWrfWLL75QVdXp06fr5ZdfrqqqP/zwg7Zp00abNGmit9xyi1atWjXrvHfp0kVTUlK0SZMmOn/+/Gyx5SZw/eOPP66jR49WVdVOnTplvTZpaWlav379w87nrl27dMSIEXrKKadocnJytuOVhNzey8AKzeN71fcv9lAsRSosVHV4ped1N5V1J1X1EOg2jlBdtapI+yqOgwdVlyxRvfZa1YQE9yqVK6d62mmqjz2munlzqYdkPIUtLHL+KIiLK34Ms2bN0pEjR6qqatu2bXXFihWqqjpp0iQdMWKEqqp+8cUXGhMTk/WFlJ6erqqqBw4c0E6dOmV9SdavX1/vvvtuVVWdOXNm1hdVZuGQKfN2RkaG1qtXT3ft2qWqrhB67rnnNC0tTTt06JC1/v7779c77rgj1/hzFhZVqlTRH374Iev+zFj37NmjjRs31j/++ENVVQGdNWuWqqrecccdWV/uXbp00e+++05VVT/++GM97bTTsh1vwYIF2r59e92/f79eccUVOmHCBFVVfe+997RZs2aqmr2w6N27t86cOVNVVR977LGsL/yMjAzdvn27qrov+OOPP14PHTqkqgUXFgcOHNCBAwfqm2++qarBFRY33XSTPvfcc6qqunXrVm3YsGHW+S0JhS0sIqrNQkR6i8hT27dvL9Ljp+8bQpUtG6n26QfI4sUcUTcOOnWCRYtKONL8lSsHp54KkybBDz/AypVw003w229wxRVQt667/8EH3f2m7ApFR4bU1FQGD3bTiw8ePJjU1FQAPvroI84/303lnpSURFJSUtZj5syZQ/PmzUlJSeHrr7/O1h4xZMiQrL/Lli3L99jly5ene/fuvPrqqxw4cIDXX3+dvn378vHHH7NmzRratWtHcnIyM2fOZOPGPNMMZdOqVatsXTinTJlCs2bNaNOmDZs2bWLdunUAlCtXjkGDBgFw/vnns3jxYnbt2sXSpUs555xzSE5O5tJLL+XXX/+eZnvdunWMHTuWOXPmUKFCBRYvXswFF1wAQJcuXUhPT2dHjm6IS5YsyTonmduC+2F98803k5SUxBlnnMHPP//M77//nu9z27t3L8nJyRx11FH8/vvvdO3aNahzAvD2229z//33k5ycTOfOndm3bx8//fRT0I8vaRHVZqElMc4iPv7vLrNLlsCZZ7rlxRddI0IpE4Hmzd1y993wzTfw0ktuGTvWLcnJro1jwABITHSPMZHpzz//5P333+err75CRDh48CAiwsSJE/N8zI8//siDDz7I8uXLqVmzJsOHD882qDCwZ0wwvWQGDx7MY489Rq1atWjZsiVxcXGoKl27ds0quAqjatWqWf8vXLiQd999l2XLllGlSpWsL8nciAiHDh2iRo0ah7V/AOzatYtzzz2X//73vxx99NGFiim38/C///2PtLQ0Vq5cSYUKFUhISChwcGZmm8WePXvo1q0bU6dOZcyYMZQvX55Dhw4B5LkPVWXevHmcdNJJhYo9VCLqyqLEHXusu6po0gT69YNnny3wIaHWqBHccou72vjhB3f1UbUq3H67C7NRI7j5ZlixwlV8mMgyd+5cLrjgAjZu3MiGDRvYtGkTDRo0YNGiRXTs2JHnn38egNWrV/Pll18CsGPHDqpWrcoRRxzB77//zptvvpltn7Nnz87629brqhcXF8fOPC6DOnXqxGeffcZ///vfrCucNm3asGTJEtavXw/A7t27+e677wr9/LZv307NmjWpUqUK3377LR9//HHWfYcOHcrqOfX888/Tvn17qlevToMGDXjxxRcB9wX7xRdfADBy5EhGjBhBhw4dsvbRoUMH/ve//wGuYKpduzbVc/QiadeuHS+88AJA1raZsR155JFUqFCBDz74IOgrJ4AqVaowZcoUJk2axIEDB0hISGDlypUAefYG69atG48++qhrLwA+//zzoI8XClZYFCQ+Ht5/31VHDRsGjzzid0RZGjSAa6+FxYvh55/hP/+BevXggQfglFNcd91rrnH3x8W5K47MxXpZlY64uPxvF1Zqair9+/fPtu7ss88mNTWV0aNHs2vXLho1asT48eNp0cLl02zWrBkpKSmcfPLJnHfeebRr1y7b47du3UpSUhKTJ0/m4YcfBtzVw8SJE0lJSeH777/Ptn1MTAy9evXizTffzOrmGR8fz4wZMxgyZAhJSUm0bds23269eenevTsHDhygUaNGjBs3jjZt2mTdV7VqVT799FOaNGnC+++/z/jx4wH3hT5t2jSaNWtG48aNeeWVV9i4cSNz587lmWeeyeqiu2LFCiZMmMDKlStJSkpi3LhxzJw587AYJk+ezNSpU2natCk///x3yrqhQ4eyYsUKmjZtyrPPPsvJJ59cqOeWkpJCUlISqampXH/99Tz++OOkpKTwxx9/5Lr9bbfdRkZGBklJSTRu3JjbbrutUMcraaIR9PMzoOvsqMx6zhKzbx8MHerqf265Be66q8zW96Snw6uvulDffhv++iv37SLopS8133zzDY0aNfI7jBKTOfdL7dq1/Q6lQNWqVWPXrl1+hxExcnsvi8hKVW2Z2/YRdWWhocw6GxsLc+bAxRfDPffA6NFw8GDJH6cE/OMfMHw4LFjgspd4V9TGGFNkEVVYhFxMDDz1FIwbB08+CYMH5/2zvYyIiwOvA8lhvPY1E8U2bNgQFlcVgF1V+MwKi8ISgfvuc/1W586F3r0hDN7EudWVn3MO7N5d+rEYY8KPFRZFdd11MH26a/w+/XTXUFCG7djx99CwQ4fg4Ydh/nzo0AE2bSrw4caYKBdRhUVxB+UV2vDhMG8efPGF+9bdvLl0jltMInD11a4RfP16aNUKPv3U76iMMWVZRBUWIW3gzkvfvvDWW66gaNcO1q4tvWMXU8+esGyZS3LYqZM1hBtj8hZRhYVvOnWCDz903Wvbt3cj5sJE48buquKUU2DIEDe4zxq+y7bCpr0OxoYNG2jSpEnxgyuh49x7770hjwVg1apVvPHG33Ow5ZeqvLBEJCv9CrisuvHx8QWmIM8Z04QJE3jwwQeLHEdxH5/JCouSkpLiRr9VrQqdO8MHH/gdUdBq14Z334URI+DOO13vKZtXo+zKTCGxevVqatWqlZX6208HDhwo0f2VZGGRX2w5v5j79OnDuHHjSuS4VatWZfXq1ezduxeAd955hzp16hT4uJwxlRVWWJSkhg1dPqljj4Xu3d2ouDBRsSJMm+Y6ec2bBx07ulHhpgSkpcHy5e5vCWvbtm3WKONVq1bRpk0bkpKS6N+/P1u3bgXcJEjXXHMNLVu2pFGjRixfvpwBAwbQsGFDbr311qx9HThwgKFDh9KoUSMGDhzIHu8Xw8qVK+nUqRMtWrSgW7duWYn6OnfuzNVXX03Lli2ZPHkyy5cvJykpieTkZMaOHVvgFcSMGTMYMGAA3bt3p2HDhtxwww0AjBs3LisB39ChQwGYNWsWrVq1ykoWeNAb4zRt2jROPPFEWrVqxahRo7jiiisAGD58OJdddhmtW7fmhhtu4NNPP6Vt27akpKRw6qmnsnbtWvbv38/48eOZPXs2ycnJzJ49mxkzZmTtY8OGDXTp0oWkpCROP/30rCR+w4cPZ8yYMZx66qkcd9xx+U7e1LNnT15//XXAjb7PTFAILiXKyJEjadWqFSkpKbzyyiu5xgSwZs0aOnfuzHHHHceUKVOy9vHQQw/RpEkTmjRpwiMB2SXuueceTjzxRNq3b8/akqoazysdbTgvRU1RXmLS01XbtHF5xZ9+2t9YiuDVV938GkcfreplUTaeQs9n8fzzqpUrqx5xhPv7/PPFjiGvtNdNmzbVhQsXqqrqbbfdpldddZWqunTYN9xwg6qqPvLII3r00UfrL7/8ovv27dM6deroH3/8oT/++KMCunjxYlVVHTFihE6cOFH379+vbdu2zZo74oUXXshKg96pU6es+RlUVRs3bqxLly5VVdUbb7wxKw15oMD05NOnT9cGDRrotm3bdO/evXrsscfqTz/9lO05qrpz3qtXL92/f7+qqo4ePVpnzpypP//8s9avX1/T09N1//792r59+6w048OGDdOzzjpLDxw4oKqq27dv14yMDFVVfeedd3TAgAFZMWQ+JuftXr166YwZM1RVddq0adq3b9+sfQ8cOFAPHjyoX3/9tR5//PF5vk5ffPGFnn322bp3715t1qxZUCnIc8Z0++23a9u2bXXfvn2alpamtWrV0v379+uKFSu0SZMmumvXLt25c6cmJibqZ599lrV+9+7dun37dj3++ON14sSJh8VX2BTlEZV1trhzcJeYWrVcvc7ZZ7sR3+np4P1qCge9esHSpdCnj+vkNXMmnHuu31GFobQ0uOgiN/G6VxXBRRe5OXkzMxsXQeav7p9//plGjRrRtWtXtm/fzrZt2+jUqRPgZmU755xzsh6TOQtc06ZNady4cVYW1uOOO45NmzZRo0YN6tWrl5U36vzzz2fKlCl0796d1atXZ6XWPnjwYLYMrpkpw7dt28bOnTuzEhGed955vPbaawU+l9NPP53MDimJiYls3LiRevXqZdvmvffeY+XKlZxyyilZz//II4/k008/pVOnTtSqVQuAc845J1vywnPOOYeYmBjAJQEcNmwY69atQ0TIyMgoMLZly5bxklc7cMEFF2Rd+QD069ePcuXKkZiYmG+a8qSkJDZs2EBqaio9c0yP+fbbb7NgwYKs9oT8UpCfddZZVKpUiUqVKnHkkUfy+++/s3jxYvr375+VtXfAgAEsWrSIQ4cO0b9/f6pUqQJQYjMARlQ1lPrRGyovVau6fBuDBsGNN7pc4mGUjKlpU/jkEzfH+KBBMGFCWIVfNmzY4Or3AlWo4NYXQ2abxcaNG1HVoNosKlWqBLg5ITL/z7ydWaefMy23iKCqNG7cmFWrVrFq1Sq++uor3n777axtAtOLF0VgLDExMbm2L6gqw4YNy4ph7dq1TJgwocB9B8Z22223cdppp7F69WpeffXVAlOLFyZuLeCD0adPH66//vpsVVCZj5s3b17W8/rpp5/yzDsWzHkKtYgqLMqcihXhf/+Df/3LNQZcdBH48CIX1ZFHwnvvuWS7d9zhsptk/kA2QUhIgP37s6/LyHDrS0Bg2uuqVatSs2ZNFnkTdT333HNZVxnB+umnn7ImP8pMAX7SSSeRlpaWtT4jI4Ovv/76sMfWqFGDuLg4PvnkE4CsFN9FVaFChaxf/6effjpz585ly5YtgJvTY+PGjZxyyil8+OGHbN26lQMHDjBv3rw897d9+/asxuUZM2Zkrc8vFfupp56aLVV5YKrzwhg5ciS33347TZs2zbY+rxTk+cUUqEOHDsyfP589e/awe/duXn75ZTp06EDHjh2ZP38+e/fuZefOnbz66qtFijsnKyxCLSYGHnsMxo93I74HDnRdbMNEpUou7AcecPM/dewIv/zid1RhIj7e9RqoXNnlhK9c2d0uRhVUToFpr2fOnMnYsWNJSkpi1apVWSm8g3XSSScxdepUGjVqxNatWxk9ejQVK1Zk7ty53HjjjTRr1ozk5GSWLl2a6+OnTZvGqFGjSE5OZvfu3RTnCv+SSy4hKSmJoUOHkpiYyN13382ZZ55JUlISXbt25ddff6VOnTrcfPPNtGrVinbt2pGQkJDnMW+44QZuuukmUlJSsv0qP+2001izZk22xuRMjz76KNOnTycpKYnnnnuOyZMnF+m51K1blzFjxhy2Pq8U5PnFFKh58+YMHz6cVq1a0bp1ay6++GJSUlJo3rw5gwYNolmzZvTo0SOr+q7Y8mrMCOfF9wbuvEyZ4jJudO6s6s3lG05eecU1fB9zjKo37XPUKXQDt6rqli2qn37q/kawnTt3Zv1/33336ZgxY0rtmBkZGdqrVy996aWXQn7MSBHVc3CXeVdeCbNmufEYp50Ga9aErEtlKPTp43oGly/vGr7z6TFoAsXHu1GPJXhFURa9/vrrWYMFFy1alK1bbqhMmDAh65gNGjSgX79+IT9mtIqoyY8ytWzZUlesWOF3GHl74w3o39/VX1er5toxpk1zQ6jDwO+/u/m+ly51g/huvbXMzgNV4iJt8iMTvaJ68qOwkVmHqAo7d7pW44suCpsrjH/+0zV8X3CBa4o57zxr+DYm0oVFYSEi/UTkvyIyW0TO9DueYtuwwTV2Bjp0CHLMdVyWxca68Rf33w+zZ7v0WN7A3ogXiVfjJroU5T0c8sJCRJ4RkS0isjrH+u4islZE1otIvslYVHW+qo4CLgPymPctjOTWpfKvv9wcGWGUY0PEDSF56SXX/NKqFXi9/yJWbGws6enpVmCYsKWqpKenExsbW6jHhbzNQkQ6AruAZ1W1ibcuBvgO6ApsBpYDQ4AY4L4cuxipqlu8x00C/qeqn+V3zDLfZgGQmuqqnipUcG0Xw4bBc8+5n+wzZrhh1GFk1SrXAJ6e7p7GgAF+RxQaGRkZbN68udiDuozxU2xsLHXr1qVChQrZ1ufXZlEqDdwikgC8FlBYtAUmqGo37/ZNAKqas6DIfLwA9wPvqOq7BR0vLAoLcG0UGza4K434eDcXxuDB7pv36qtdHU/AyM2y7rffoF8/N/K7YsXsF09xcW62PmNM2VUWG7jrAIGTeW721uXlSuAMYKCIXJbbBiJyiYisEJEVaWHSUHxYl8qTTnKzEV15JTzyCJx6Kqxb52uIhXHUUbBwIQwdengtWxADUo0xZVhYNHCr6hRVbaGql6nqE3ls85SqtlTVlvHh3J89NhamTHETZP/4IzRv7lKGhInYWFcNlZupU90MtF52aWNMGPGrsPgZCEwtWddbVyylPgd3KPXt675Zk5Ph/PPdzES7dvkdVVDyGnNxxRXu6dSqBT16wL33wkcfWbdbY8KBX4XFcqChiDQQkYrAYGCBT7GUXfXquRn3xo93/VRbtnTtGWEgLi777WrVXPPMrFlu7OGmTXDLLa7L7RFHuBq3G25wiXrT030J2RiTj9LoDZUKdAZqA78Dt6vqNBHpCTyC6wH1jKreU1LHDJsG7sL44AN3hZGeDpMmuUy2YT5s+s8/3SjwxYvdsnz5320djRq5lCLt27slISHsn64xZZ7vvaFKW0QWFuB6Tw0f7tKF9OvnUoR4E79Egn37YMUKWLTIFR5LlkBmjeIxx/xdcLRvD0lJLqGvMabkRE1hETBT3qh1YdSLqFAOHYLJk91ouKOOguefd9+eEejQIfj667+vPBYtctVX4Kq5Tj3VPfV7783e7mHddI0pmqgpLDJF7JVFoBUr3JiMH390MxPddFNU/NT+6ae/C4/Fi2H16txn8IvAt7UxIRc1hUVUXFkE2rEDLrvMjQY/7TTXenzMMX5HVaq2bs29Ju7QIWvjMKawyuKgvJDQsjQHd2moXt2NwXjmGTdsulkzePNNv6MqVTVr5r7+7LOtV5UxJSmiCouoJOLGYKxY4a4qevaE668/fAh1BMvZTbdSJXjtNVd2LlzoS0jGRJyIKiwialBeYTVqBB9/7LrUTpoE7dqFVcrz4tixw7VRZC779rmsKVWqQJcubjxHRobfURoT3iKqsIi6aqicKld2OTXmzYP16yElBV54we+ofNGiBXz2GYwc6XpLdegAP/zgd1TGhK+IKiyMZ8AAN9K7aVM3XPrii2H3bjdOI4zm/C6uatXg6afd5EzffutSjcya5XdUxoQnKywiVf368OGHrg7mmWdcRtt69aBrV3dfaqrfEZaac891abaaNXNTwV5wgY3DMKawIqqwiOo2i9yULw933w0vvuhm4PvrLzckOszm/C4J9eu7jCl33OHGMSYnuyYeY0xwIqqwiPo2i7wce+zhXYYqVHCZ/aJI+fIuJ+NHH7lxGJmjvy1lujEFi6jCwuQhIQEOHMi+bv9+tz4KtWvnmnQGDnS1dGecAZs3+x2VMWWbFRbRID7eJR2sXNn1J4XsM/RFoRo1XLPN9OmuzT8pCV5+2e+ojCm7rLCIFkOGwMaNbpTatde6rHxR3jVIxCXx/fxzOO4414nssstgzx6/IzOm7LHcUNHowAGXS2rVKjcYoWFDvyPy3f79cNtt8MADbnxjaqrrPWVMNLHcUCa78uVdTqkKFdwVRxSlBslLxYrw73/DO+/Atm3QqpXLBB9Bv6WMKZaIKixMIRx7rGvHWLkSbr7Z72jKjDPOcGMyunWDq6+Gs86CLVv8jsoY/1lhEc3694fRo10uqSjLVpuf+Hh45RV47DF4/33X+P3WW35HZYy/rLCIdpMmubQgw4bBr7/6HU2ZIQKXX+6S+dauDd27w3XXuXGNxkQjKyyiXeXKLtngrl1w4YVutJrJ0qSJ61p7+eXw0EPQti2sXet3VMaUvogqLCzdRxElJrrW3Hffdd2BTDaVK7sqqQUL3LSuJ5/srjwyl+rV/Y7QmNCLqK6zmaJiDu6SpgqDBsFLL7nJrdu08TuiMumXX6BOncPXR+DHyESh/LrOli/ggT8CCoj397BNvPWPqOqU4gZqfCQCTz0Fn37qutN+/rkb5myyibIpzo3Jkm9hoaoNSisQUwZk5sDo0AEuvdS1ZYj4HVVYePdd1+3WmEhVYJuFiMSIyAelEYwpA9q2hbvugjlz3DgMc5icCXzLlXPjMiZNsuooE7kKLCxU9SBwSERsWHS0uPFGOP10GDMG1qzxO5oyJ+ec39u3uyEr118P559vuaVMZAq2N9Qu4CsRmSYiUzKXUAZmfFSuHDz3nJuXdPBgN1mSyVO1am5+qXvucbV47du7nI3GRJJgC4uXgNuAj4CVAYuJVEcfDTNnwldfuZ/MJl8iLmvKq6/C999Dy5Yuwa8xkSKowkJVZwKp/F1IPO+tKxUi0khEnhCRuSIyurSOG/V69HDpzP/zH5vsIUhnneU6lNWu7Rq8H33U2jFMZAiqsBCRzsA6YCrwH+A7EekY5GOfEZEtIrI6x/ruIrJWRNaLyLj89qGq36jqZcC5QLtgjmtKyL33QvPmbs7un37yO5qwcNJJ8MknruAYMwZGjoR9+/yOypjiCbYaahJwpqp2UtWOQDfg4SAfOwPoHrhCRGJwBU8PIBEYIiKJItJURF7LsRzpPaYP8DrwRpDHNSWhUiXXhTYjA4YOPXx6VpOr6tXdxdjtt8OMGdCxo03dasJbsIVFBVXNyoijqt8BFYJ5oKp+BPyZY3UrYL2q/qCq+4EXgL6q+pWq9sqxbPH2s0BVewBDg4zZlJSGDV1V1OLFcPfdfkcTNsqVgwkTXKHxzTfQooU7hcaEo2ALi5Ui8rSIdPaW/wLFyadRB9gUcHuzty5X3jGniMiT5HFlISKXiMgKEVmRlpZWjNBMri64wC133QUffuh3NGGlXz9XLXXEEW6CwieesHYME36CLSwuA9YAY7xlDVBqDc2qulBVx6jqpao6NY9tnlLVlqraMj4+vrRCiy5Tp7rJqocOhfR0v6MJK4mJruH7zDPdFCKXXmrpzk14CWoEN/CFqj6kqgO85WFVLc5b/WegXsDtut66YrGssyEWF+faL7Zsca229vO4UGrUcJlrb7kF/vtfd5VhU4iYcBHsCO61InJsCR53OdBQRBqISEVgMLCguDu1ObhLQYsWbrLqBQvclYYplJgY1+zz4ovw5ZfudH78sd9RGVOwYKuhagJfi8h7IrIgcwnmgSKSCiwDThKRzSJykaoeAK4A3gK+Aeao6tdFeQI5jmVXFqXh6quhZ083WO+LL/yOJiwNHAjLlrm5Mjp1sjRcpuwLaj4LEemU23pVLZMtnTafRSlIS4NmzVwf0ZUroWpVvyMKS3/+6TKqvPMO/Otf8PDDULGi31GZaJXffBbBtlk8qaof5lxKPNJisiuLUhQfD7NmwXffuZFnpkhq1YI33oCxY13v5DPOgN9/9zsqYw7nV5tFSFibRSnr0gVuugmeecY1fJsiKV/ezWb7/POwYoXLK2UXxqasCXmbRWmyKwsfTJjg5sC45BL44Qe/owlrQ4bA0qWuEbx9e3j2Wb8jMuZv1mZhim/DBkhOdkmRFi+GCkEN7jd5+OMPOPdc+OADuOoqmDjRTqkpHUVusxCRkyGrUPg4R3uFDSkyTkICPP20G3V2661+RxP2ateGt992nc4mT3Y9pkT+XqpX9ztCE40KqoZ6PuD/ZTnu+08Jx1JsVg3lo4EDXVXUAw+4bzpTLOXLu55RM2fCwYPZ79u505+YTHQrqLCQPP7P7bbvrIHbZw8/7PJaXHihdekpIRde6HcExjgFFRaax/+53TbRrkoV1ytq+3YYNswVGMuXuzEZpkQdOuR3BCbalC/g/rreXNsS8D/e7TyzxJoo1rQpPPSQG2FWt64brLd/vxuiPGSI39GFpbi4w6ueBg1y82TYWEhTWvLtDSUiw/J7cGlOrRoMEekN9D7hhBNGrVu3zu9woteWLW4O78Cfv5Urw8aNbjCfKTJVVxbfcAMkJcErr8CxZX4ElAkX+fWGCqrrbLixrrM+W77cDUXesePvddWrw7vvwimn+BdXBHnzTZcmJDbWTa506ql+R2QiQbHSfRhTaAkJbhrWQBkZbr0pET16uGy11au7VOczZvgdkYl0VliYkhcf79ooYmPd7ZgYd9uqoEpUo0ZuBr6OHWHECLj2Wpsi3YSOFRYmNIYMgZ9+gr593aCB00/3O6KIVKuWq5K68krXc7lXL9i2ze+oTCQKqrAQkXgRuVlEnhKRZzKXUAdXWDYor4yJj4f773fzh9pESSFTvjxMmQJPPQXvvw9t2rhkwMaUpGBzQy0FFgErgazxpKo6L3ShFZ01cJcxffvCkiXuSqNKFb+jiWgffQRnn+2qo2bPdnN+GxOskmjgrqKqN6rqHFWdl7mUYIwmko0dC+npMH2635FEvI4dXWe0evVcI/jkyTZVuikZwRYWr4lIz5BGYiJXu3bQurUbIJAz0ZEpcQkJLtV5nz4uGeGoUa4m0JjiCLawuApXYOwTkZ3esqPARxkDLlXq2LFuvouXXvI7mqhQrRrMm+eSAE+b5voXbNnid1QmnAVVWKhqnKqWU9VY7/84VbVEySZ4/frBCSe4yRmsXqRUlCsHd93l0nWtXOnGQ65a5XdUJlwF3XVWRPqIyIPe0iuUQZkIFBPjBgIsXw6LFvkdTVQZNMjNSXXwoKsRtIs7UxTBdp29H1cVtcZbrhKR+0IZWFFY19kybvhwN7PPxIl+RxJ1WrRw5XTTpq631J132gWeKZxgryx6Al1V9RlVfQboDpwVurCKxuazKOMqV4YrroDXXoM1a/yOJuocfTQsXOjmyLj9dnfFsXu331GZcFGYEdw1Av63b2NTNP/6l0sDMmmS35FEpdhYl0dq4kSYOxc6dHDDX4wpSLCFxX3A5yIyQ0Rm4gbn3RO6sEzEio93iYxmzYJff/U7mqgkAtdf7y7wvv/eNXwvXep3VKasC7Y3VCrQBngJmAe0VdXZoQzMRLBrr3VZaKdMKXhbEzI9e7rMtXFxlrnWFCzfwkJETvb+NgeOBjZ7yzHeOmMK74QTYMAAePzxw6eAM6WqUSP49FNXHTViBFx3nY2bNLkr6MriWu/vpFyWB0MYl4l0Y8e6ubqnTfM7kqgXmLn2oYdcYkIRt1S30VTGE2wiwVhV3VfQurLCEgmGiY4d3VSr69dDhQp+R2NwBURO1sU2epREIsHcmr9KtUlMRKqKyAobEBhBrr/edcV58UW/IzHGFKCgNoujRKQFUFlEUkSkubd0BoLKNe3NfbFFRFbnWN9dRNaKyHoRGRfErm4E5gRzTBMmevWCk0+2FCBl3NSp9vIYKF/A/d2A4UBd4KGA9TuBm4M8xgzgMeDZzBUiEgNMBbriGsyXi8gCIAbXTTfQSKAZbuR4bJDHNOGgXDnXojpqFLz3Hpxxht8RRb24uOx9DmJi3DjK1atd5zWrLYxewbZZnF2c+StEJAF4TVWbeLfbAhNUtZt3+yYAVc01hYiI3ANUBRKBvUB/VT2UY5tLgEsAjj322BYbN24sarimNO3b53JqN2sGb73ldzQmh4MH4ZZb4N//hs6d3UC+f/zD76hMqOTXZlHQlQXgZsQTkbOAxgT8ulfVO4sYUx1gU8DtzUDrfI5/C4CIDAf+yFlQeNs8BTwFroG7iHGZ0hYbC2PGuG+kL7+EpCS/IzIBYmLczLiNG8PFF0OrVvDqq5CY6HdkprQFm0jwCWAQcCUgwDlA/RDGlStVnaGqr+V1vyUSDFOjR0PVqvCg9cYuqy64wOWV2r3bzfH9xht+R2RKW7C9oU5V1QuBrap6B9AWOLEYx/0ZqBdwu663rlgskWCYqlnT/WxNTYVNmwre3viibVs3gO/4413fhEmTrOE7mgRbWOz1/u4RkWOADNyI7qJaDjQUkQYiUhEYDCwoxv4Au7IIa1df7b55Jk/2OxKTj2OPdXNjDBjgej6PHGlTtkaLwszBXQOYCHwGbABSg3mgiKQCy4CTRGSziFykqgeAK4C3gG+AOar6dSFjP4xdWYSxhAQ45xx46ik3stuUWVWrwpw5MH68yyfVpYtN2RoNguoNle0BIpWAWFUtc59oEekN9D7hhBNGrVu3zu9wTGF99pmbpeff/4YbbvA7GhOE2bPdnFZHHgkLFrhObSZ8FXsEt4hc7l1ZoKp/AeVE5F8lF2LJsCuLMNe8ufuZOnky7N/vdzQmCIMGuVlyDxxwU7bOn+93RCZUgq2GGqWq2zJvqOpWYFRIIjLRbexY+OUX19htwkLLlm7K1sRE6N8f7r3XGr4jUbCFRYzI3ynGvBHYFUMTUtFZA3cE6NYNmjRx3WjtGydsHHMMfPghnHeeGzIzdCjs3Vvw40z4CLaw+D9gtoicLiKn4xq3/y90YRWNVUNFgMxp3Favhv8rc28xk4/Kld0EiPfe6y4MO3VyF4kmMgSb7qMccClwurfqHeBpVS2T06RYivIwt38/HHccnHgivP++39GYIpg/H84/H444Al55xVVVmbKv2A3cqnpIVR9X1YHe8mRZLCisGipCVKwIV10FH3wAK1f6HY0pgn79YMkSN5FShw6u15QJb8H2hmonIu+IyHci8oOI/CgiP4Q6uMKyaqgIcsklLgWqpQAJW82auYbvFi1g8GA3LuPQYVndTLgIts1iGi5FeXvgFKCl99eY0DjiCLj0Ujcx0oYNfkdjiujII132+REj4K673LjL3bv9jsoURbCFxXZVfVNVt6hqeuYS0siKwKqhIsxVV7kG74cf9jsSUwyVKrmp1idNcm0Z7du7CRJNeAm2gft+3MRELwFZmWBU9bPQhVZ01sAdQS68EObNcwkGa9XyOxpTTG+8AUOGuJ5TL7/skhOasqMk5uBujat6uheY5C1WmWxC7/rrYc8eePxxvyMxJaBnT1i2DKpVc5MpPftsgQ8xZUShc0OFA7uyiDDdu8OqVa7tItZm1o0E6emu/eKDD7Kvj4uDHTv8icmUwEx5IjI+t/XFmCnPmOCNHevm5541y817YcLeP/7hZtGtmCMPROD836ZsCbYaanfAchDoASSEKKYiswbuCNWlC6SkuG601vcyYlSo4HcEpjCCHZQ3KWC5B+gMHBfSyIrAxllEKBF3dbF2LbyW56y6JkJ89JHfEZjcBHtlkVMV3FSoxpSOgQPdNG0TJ/odiSlBcXHZb5cr52ocp03zJx6Tt2BHcH8lIl96y9fAWuCRkEZmTKAKFeCaa9ycnh9/7Hc0poTs2OGSC2cu6elw2mmuaeqaa9w8GaZsCHacRf2AmweA372pUcsk6w0VoXbtgnr1XBvGvHl+R2NC5MABuO46mDLFZayfPdsN6DehV6xxFt7cFW+p6kZv+bksFxQmglWrBqNHu9Fc69f7HY0JkfLl3WSJTz7pUoW0aQM2S7L/CiwsvOyya0Xk2FKIx5j8XXmlq5J66CG/IzEhdskl8O67kJYGrVtbtnq/5VtYiMgA79+awNci8p6ILMhcQh9e4VjX2Shw9NFwwQUwfbr7FjERrVMn+PRT97KfeaYN5PdTvm0WIvKZqjYXkU653a+qH4YssmKwNosI9803bsLn22+HCRP8jsaUgh073JStr78O//oXPPKIjdMIhZKY/OjD3JaSDdOYIDVqBL16wWOPubxRJuJVr+5m3Lv+evjPf6BHD/jzT7+jii4FFRYnB3SZPWwplQiNyc3Ysa6f5YwZfkdiSklMjBtmM306LFrk2jG+/dbvqKJHQdVQXwM987pfVTeGIqjismqoKKDqusmkp7uR3TExfkdkStGSJdC/v5uuffZs18XWFF9xqqH2B3SZPWwJQazGBCczBcj337sZdUxUadfOTdlav75Lez55svv9YEKnoMJiSalEYUxR9O8Pxx/v6ibsmyLq1K/vrjD69IGrr3Zdbffv9zuqyFVQYfGAiBwbxFK9VKI1JlBMDFx7LXzyiUsDYqJOtWpuMP/NN8PTT0PXrvDHH35HFZkKarP4AFBA8tmHAjNUNWRzXolIZ+Au4GvgBVVdmN/21mYRRfbscQkGTz0VFpS5oT+mFD3/PIwcCccc494KTZr4HVH4KfLkR6p6Wgkc/BmgF7BFVZsErO8OTMbN7f20qt6fXyjALiAW2FzcmEwEqVIFLr8c7rwT5sxxWeji4/2OyvjgvPNcrWS/fm5u79RU18PalIyipigvjBlA98AVXr6pqbhJlBKBISKSKCJNReS1HMuRwCJV7QHcCNxRCjGbcHLMMe7v+ee7iuzUVH/jMb5p3do1fJ94omvLsOaskhPywkJVPwJyDp9pBaxX1R9UdT/wAtBXVb9S1V45li2qmjk92lagUm7HEZFLRGSFiKxIszQQ0SMtzeWyBsjIgL174aKLLBVIFKtb143DGDgQbrgBRoyAv/7yO6rwVxpXFrmpA2wKuL3ZW5crERkgIk8CzwGP5baNqj6lqi1VtWW8VUNEjw0bDp/IuUIFt95ErSpV4IUXXDaYmTNdVvvff/c7qvDmV2FRKKr6kqpeqqqD8mvctkSCUSgh4fD+kn/95dabqFaunEsfNmcOfP45tGoFX3zhd1Thy6/C4megXsDtut66YrE5uKNQfLybg7NyZdePEiAlxRq5TZZzznHVUgcPQnKyG8+ZuVS3Tv9B86uwWA40FJEGIlIRGAwUu9+jXVlEqSFDYONGN+HBNde4aVc/+8zvqEwZ0qKFa/jOaefO0o8lXAU1rWqxDiCSCnQGagO/A7er6jQR6YmbxzsGeEZV7ympY9o4iyi2bRscd5zrFvPmm35HY8oYyWXEmPWW+luxU5QXh6oOUdWjVbWCqtZV1Wne+jdU9URVPb6kCgq7sjDUqAHjxsH//R989JHf0Zgw8OuvfkcQHsKigTtY1mZhALjiCje12k032c9Gk01c3OHrTjkFVq4s/VjCTUQVFsYArt/k+PGwdKmbWs0Yz44d7vdD5rJqlUsx1qEDvPii39GVbRFVWFg1lMly0UUu98Mtt8ChQwVvb6JSs2Zuju+UFDj3XLjjDnu75CWiCgurhjJZKlRw+aK+/NKNzjImD//8p+tId+GFbhDf4ME2W29uIqqwMCabwYMhKclVSWVk+B2NKcMqVXIz9E6cCHPnumqpzZayNJuIKiysGspkU64c3HOPm01v2jS/ozFlnAhcf71Lb/7dd27E96ef+h1V2RFRhYVVQ5nDnHWWm+vizjutbsEEpVcvWLYMYmOhY0c3T4aJsMLCmMOIwH33uc70j+Wag9KYwzRp4q4qWreGoUPh1lut4TuiCgurhjK56tgRuneH++93I7yNCULt2vDOO3Dxxa42c+BA2LXL76j8E1GFhVVDmTzdey9s3QoPPuh3JCaMVKwITz0FjzwCr7wC7du7NGTRKKIKC2PylNmR/pFHbGIDUygicNVVbnznjz+6hu+lS/2OqvRZYWGix113wb59rk7BmELq3t0lNK5e3U31/uyzfkdUuqywMNHjxBPdHJtPPGEz6ZkiadQIPvnEVUcNGwY33ujmyYgGEVVYWAO3KdDtt7vxFxMm+B2JCVO1armkxqNHwwMPQL9+0TEvRkQVFtbAbQpUt67LSvvcc7Bmjd/RmDBVoQL85z8wdaqbNuXUU117RiSLqMLCmKCMGwdVq7rO88YUw7/+5a4yNm92qc4jeQoVKyxM9Kld2+V1ePlly+dgiu2MM9zbqHZt93+kZpaxwsJEp2uugfh4uPlmvyMxEaBhQ9dT6rTT3CC+a66BAwf8jqpkWWFholNcnCso3nvPLcYUU40abizGmDFuOE/v3hBJfW1EI3DayZYtW+qKFSv8DsOUdfv2ue60Rx/tfhaK+B2RiRBPPQWXXpp9XVycm6mvLBORlaraMrf7IurKwrrOmkKJjXVdaD/9FObP9zsaE0EuueTwdeHevdauLEx0O3AAmjZ1Yy++/NJNyGxMCcjtQrWsf91GzZWFMYVWvrxLA7JmDcya5Xc0JsJdcUX4NnxbYWHM2WdDixZudPdff/kdjYkQcXHZb1eo4Abx9ejhEiCHGyssjBFxKcw3bnQtk8aUgB07XLVT5rJ/PzzzDHz4oZtUae1avyMsHCssjAHo2hU6d4a7747uGW5MSI0YAe+/7+bgat0a3n7b74iCZ4WFMfD39KtbtsDkyX5HYyJY+/auA96xx0LPnvDoo2W/4RussDDmb23aQJ8+MHEi/Pmn39GYCJaQAEuWwFlnuUF8l10GGRl+R5W/sCgsRKSciNwjIo+KyDC/4zER7J57XGXzv//tdyQmwsXFufRk48a5prIzz4T0dL+jylvICwsReUZEtojI6hzru4vIWhFZLyLjCthNX6AukAFsDlWsxtCkCQwdClOmwC+/+B2NiXDlyrnaz+eeg2XL3JStZTVzfmlcWcwAugeuEJEYYCrQA0gEhohIoog0FZHXcixHAicBS1X1WmB0KcRsotkdd7jO8Hfd5XckJkqcfz4sXAi7d7va0Dfe8Duiw4W8sFDVj4CcFcCtgPWq+oOq7gdeAPqq6leq2ivHsgV3NZHZMznXSQxF5BIRWSEiK9LS0kL1dEw0OO44l6/h6afh++/9jsZEiTZtYPlyOP54l4TwoYfKVsO3X20WdYBNAbc3e+vy8hLQTUQeBXKdXkRVn1LVlqraMj4+vuQiNdHp1lvdKKrx4/2OxESRevVg8WLo3x+uuw4uuqjsjBMNiwZuVd2jqhep6pWqOjWv7SyRoCkxRx8NV10FqakuZ5QxpaRqVZgzB267DaZPdxMqbdnid1T+FRY/A/UCbtf11hWLzcFtStQNN8ARR8Att/gdiYky5crBnXfCCy/AihWu4furr3yOyafjLgcaikgDEakIDAYWFHendmVhSlTNmq7AeO011ynemFI2aBAsWuTGYJx6Kiwo9rdk0ZVG19lUYBlwkohsFpGLVPUAcAXwFvANMEdVvy7usezKwpS4MWPgn/90s+qVpdZGEzVatnQjvk8+Gfr1g/vv9+etGFHzWYhIb6D3CSecMGrdunV+h2MixdSpLrf0m29C9+4Fb29MCOzZAyNHwuzZcMEFbiBfbGzJHiO/+SwiqrDIZJMfmRK1f7/7WVejhqtALhcW/UJMBFJ1uS7Hj3ddbV9+GY46quT2b5MfGVMcFSu6gXqffw5z5/odjYliIq6X1Ny5rpNeq1bubVkaIqqwsAZuEzLnnQeNG7tParhOdWYixtlnu/EYqi6L7bx5oT9mRBUW1sBtQiYmxiUZ/O47mDHD72iMISXFjfhOSoKBA112mlC2KlibhTHBUnX9FzdvhnXrSr510Zgi2LcPRo06fAr5uDiXQLkwoqbNwqqhTEhlTr+6eTM88ID7WWd5yIzPYmPh2WcPX79zZ8kex64sjCmspCQ3nLZ6dTdaato0GDLE76hMlBM5fF1hv96j5srCmJBLS3PtFuCu8ffuddne7ArDRDgrLIwpjA0bDm+rqFDBrTfGR3Fx+d8urogqLKzNwoRcQoIbpBcoI8OtN8ZHO3a4aqfMpbCN2wWJqMLCus6akIuPd20UlSu7NovKld1tm0PFRLjyfgdgTNgZMsRNMrBhg7uisILCRAErLIwpivh4KyRMVImoaihrszDGmNCIqMLC2iyMMSY0IqqwMMYYExpWWBhjjCmQFRbGGGMKFJG5oUQkDdhYjF0cARSnlbwwjw9m2/y2Kex9wayrDfxRQEwlyc536Z3v4p7rwu4jVOe7MOvtfAe/TX1Vzb2bn6rakmMBniqtxwezbX7bFPa+YNYBK+x8R+b5Lu65LivnuzDr7XyXTExWDZW7V0vx8cFsm982hb0v2HWlyc536SmJY5eF812Y9Xa+i76/LBFZDWWKR0RWaB5pik3Js/Nduux8F41dWZjcPOV3AFHGznfpsvNdBHZlYYwxpkB2ZWGMMaZAVlgYY4wpkBUWxhhjCmSFhSmQiBwnItNEZK7fsUQDEeknIv8Vkdkicqbf8UQyEWkkIk+IyFwRGe13PGWZFRZRSkSeEZEtIrI6x/ruIrJWRNaLyDgAVf1BVS/yJ9LIUMjzPV9VRwGXAYP8iDecFfJcf6OqlwHnAu38iDdcWGERvWYA3QNXiEgMMBXoASQCQ0QksfRDi0gzKPz5vtW73xTODApxrkWkD/A68EbphhlerLCIUqr6EfBnjtWtgPXelcR+4AWgb6kHF4EKc77F+Tfwpqp+VtqxhrvCvrdVdYGq9gCGlm6k4cUKCxOoDrAp4PZmoI6I/ENEngBSROQmf0KLSLmeb+BK4AxgoIhc5kdgESiv93ZnEZkiIk9iVxb5sjm4TYFUNR1Xf25KgapOAab4HUc0UNWFwEKfwwgLdmVhAv0M1Au4XddbZ0LDznfpsXNdTFZYmEDLgYYi0kBEKgKDgQU+xxTJ7HyXHjvXxWSFRZQSkVRgGXCSiGwWkYtU9QBwBfAW8A0wR1W/9jPOSGHnu/TYuQ4NSyRojDGmQHZlYYwxpkBWWBhjjCmQFRbGGGMKZIWFMcaYAllhYYwxpkBWWBhjjCmQFRbGBElEVEQmBdy+XkQmBPnYhMyU2SLSUkQsnYcJK1ZYGBO8v4ABIlK7ODtR1RWqOqaEYjKmVFhhYUzwDgBPAdfkt5GIdBKRVd7yuYjE5bi/s4i85v1fTUSmi8hXIvKliJztrT9TRJaJyGci8qKIVPPW3y8ia7xtHwzN0zTmcJZ11pjCmQp8KSIP5LPN9cDlqrrE+5Lfl8+2twHbVbUpgIjU9K5cbgXOUNXdInIjcK2ITAX6AyerqopIjZJ4QsYEw64sjCkEVd0BPAvkV420BHhIRMYANby8RHk5g4DZ8FR1K9AGN5vbEhFZBQwD6gPbcQXPNBEZAOwpxlMxplCssDCm8B4BLgKq5nanqt4PXAxUxn3hn1zI/Qvwjqome0uiqmYmw2sFzAV6Af9X1CdgTGFZYWFMIanqn8AcXIFxGBE5XlW/UtV/41Jj51dYvANcHvDYmsDHQDsROcFbV1VETvSqtI5Q1Tdw7SbNSuQJGRMEKyyMKZpJQF69oq4WkdUi8iWQAbyZz37uBmp6238BnKaqacBwINXbxzJcgRMHvOatWwxcWzJPxZiCWYpyY4wxBbIrC2OMMQWywsIYY0yBrLAwxhhTICssjDHGFMgKC2OMMQWywsIYY0yBrLAwxhhTICssjDHGFOj/AayGaeIlyw1EAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from numpy import abs\n",
    "\n",
    "x_N_aintTR = []\n",
    "x_N_intRI = []\n",
    "\n",
    "#calculate x values for the adaptive trapezoidal rule\n",
    "for i in range(1,len(myaintTR.error_list)+1):\n",
    "    x_N_aintTR.append(2**i)\n",
    "\n",
    "#get only the last errors calculated per N in the Romberg integration method\n",
    "y_error_intRI = []\n",
    "erroriterator = 0\n",
    "looplimit = len(myintRI.I_list)\n",
    "for i in range (0,looplimit):\n",
    "    newlooplimit =  len(myintRI.I_list[i])\n",
    "    for j in range(0,newlooplimit):\n",
    "        if (j!=i):\n",
    "            erroriterator += 1\n",
    "        if (j == newlooplimit-1 and j!=0):\n",
    "            y_error_intRI.append(abs(myintRI.error_list[erroriterator-1]))\n",
    "            \n",
    "#calculate x values for the Romberg integration method\n",
    "for i in range(1, len(y_error_intRI) + 1):\n",
    "    x_N_intRI.append(2**i)\n",
    "\n",
    "#create the graph\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "#add graph values for Adaptive Trapezoidal Rule\n",
    "ax.scatter(x_N_aintTR, myaintTR.error_list, s = 20,c = 'b', marker = \"s\", label = \"Adaptive Trapezoidal Rule\")\n",
    "ax.plot(x_N_aintTR, myaintTR.error_list, c = 'b')\n",
    "\n",
    "#add graph values for Romberg Integration method\n",
    "ax.scatter(x_N_intRI, y_error_intRI, s = 20, c = 'r', marker = \"o\", label = \"Romberg Integration Method\")\n",
    "ax.plot(x_N_intRI, y_error_intRI, c = 'r')\n",
    "\n",
    "ax.set_xlabel(\"N slices\")\n",
    "ax.set_ylabel(\"|Truncation Error|\")\n",
    "\n",
    "#display graph in log scale\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "plt.legend(loc = 'upper right')\n",
    "plt.title(\"Absolute Value of Truncation Error vs N Slices\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "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": "markdown",
   "metadata": {},
   "source": [
    "1. Here I copy the `gaussxw` function from Appendix E of the Newman text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "A2vxCvzvVe60"
   },
   "outputs": [],
   "source": [
    "from numpy import ones,copy,cos,tan,pi,linspace\n",
    "\n",
    "def gaussxw(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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Here, I copied the `gaussxwab` function and renamed it to `gaussmap`which is found in Appendix E of Newman and maps the calculated weights and x values from the `gaussxw` function to new appropriate weight and x values given a range `a` to `b` and `N` number of slices. I then calculate the integral using the equation $$I(a,b) \\approx \\sum_{k=1}^{N}w_k^{'}f\\left(x_k^{'}\\right) $$, where $w_k^{'}$ and $x_k^{'}$ are calculated using the $gaussmap$ function. I then use a `for` loop to add all the necessary terms to `intGQ1` and print out my results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Integral for Problem 1 using the Gaussian quadrature method:  1450.000000000002\n"
     ]
    }
   ],
   "source": [
    "def gaussmap(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",
    "\n",
    "x_prime, w_prime = gaussmap(3,0,5)\n",
    "\n",
    "intGQ1 = 0\n",
    "for i in range(0,len(x_prime)):\n",
    "    intGQ1 += w_prime[i]*function1(x_prime[i])\n",
    "print(\"Integral for Problem 1 using the Gaussian quadrature method: \", intGQ1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "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": "markdown",
   "metadata": {},
   "source": [
    "1. Here, I've written a function named `function3` to stand in for the integrand of the equation above. I then define a function named `period` with the parameter `a`. The `period` function maps the weights and appropriate x values using the `gaussmap` function for $N = 20$ slices, from 0 to an argument `a`. Afterwards, the integral is calculated using the Gaussian quadrature method and summed to `intGQ2`. The function then returns `intGQ2` multiplied by the constant $\\sqrt{8m}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "p28JF-OQYock"
   },
   "outputs": [],
   "source": [
    "from numpy import sqrt, linspace\n",
    "\n",
    "def function3(x,a):\n",
    "    return sqrt(a**4 - x**4)**-1\n",
    "\n",
    "def period(a):\n",
    "    m = 1\n",
    "    x_prime, w_prime = gaussmap(20,0,a)\n",
    "    intGQ2 = 0\n",
    "    for i in range(0,len(x_prime)):\n",
    "        intGQ2 += w_prime[i]*function3(x_prime[i],a)\n",
    "    return sqrt(8*m)*intGQ2\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. I then calculate $16$ equally spaced values for the amplitude $a$ from $0$ to $2$ using the `linspace` function to set as x values for the graph. I pass each x value into the `period` function to calculate the associated period for `a`. I then construct the graph of period vs amplitude.\n",
    "\n",
    "    From the graph, it can be clearly seen that as amplitude increases, the period of the Anharmonic Oscillator decreases. This is interesting, because in a simple harmonic oscillator, the period does not depend at all on the amplitude of the wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\hydroforth\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in double_scalars\n",
      "  after removing the cwd from sys.path.\n",
      "c:\\users\\hydroforth\\appdata\\local\\programs\\python\\python37\\lib\\site-packages\\ipykernel_launcher.py:11: RuntimeWarning: invalid value encountered in double_scalars\n",
      "  # This is added back by InteractiveShellApp.init_path()\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Period (T)')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAt8ElEQVR4nO3deZxV9X3/8ddn9n2GTWBmWBQURRhASUzURNRU4xaxNf1pbKpJrLHZW2urNYttTGJrliYxiVGTYjaNSVxojVtcYhPcICyDILKEbdhGYGBghlk/vz/OuXAY7mwwd5m57+fjcR/3rPd87vee+znnfr/nfo+5OyIikjmyUh2AiIgklxK/iEiGUeIXEckwSvwiIhlGiV9EJMMo8YuIZBgl/hQyMzezyamOo7/M7F/N7P5UxxETjcfMJoblmhOOv2hm16c2woFnZuPNbJ+ZZR/FunPMbHMi4kqV6HfJzO4xsy8M4GvPM7M7Bur10oES/wAIk8tuM8tPdSzJ4O5fdfejSqZmdqmZvWZm+81sp5n93MyqUxVPl9jWm9n7jvV1enj9OWbWGSbsRjNbZWYfOZrXcveN7l7i7h0DHedAMrPbw6R8RrK26e43uvuXw+0PuYPcQFDiP0ZmNhF4D+DAB1IUQ04qtttfZnYl8Avgv4CRwKlAC/AHMxuWwtCOWT8+gy3uXgKUAf8C3GdmUxO0rZQyMwP+FtgVPkuaUOI/dn8LvALMA66Nzgh/In7PzJ4Iz/BeNbNJXdZ/n5mtNrOGcFkL151kZs+HZ8Vvh2fGFZHXXm9m/2Jmy4D9ZjY5PLP6iJltCn+B3Ghm7zCzZeHr3x1ZP8vMPm9mG8xsh5n9xMzKw3mx6pJrzWxjuP3bIuvebmY/i4yfbWYLwm1sMrPruhZS+L6+Adzh7r9w92Z33wZcD+wD/iFcbrKZ/d7M9oTb/WXkNU41s2fNbJeZbTezf40XT3d6KlMz+ykwHvif8Iz8n8PpHzCzN8L39qKZndLDZ9DnhOyBx4DdwNTw87jFzNaG8T1sZsPD7cQ+j4+Z2UbgeTuySqvSzOaHZbPGzP4uEmdhuC/uNrMVwDt6KKMfmNnXu0x73Mz+MRz+FzOrs0O/WM7v4W2+BxgLfAa4yszyIq95nZn90cy+FZbtOjM7M5y+Kdwnr40sP8+CKpxnw23/3swmdPMe5pnZHWZWDDwJVIaf6b6wnA6rurEuvwrMbJaZ/Snczi+Bgi6vf6mZLQnjXmBmNT2UQXpydz2O4QGsAT4BnA60AaMj8+YBO4F3AjnAz4GHIvMd+F+ggiDp1APvD+dNBv4CyAdGAS8B/xVZdz2wBBgHFAITw9e7h2BHvQA4ADwGHAdUATuAc8L1PxrGfgJQAjwC/DScF3ut+8LXnkFwZn5KOP924Gfh8ASgEbgayAVGADPjlNPJ4WseH2fevwEvh8MPArcRnJQUAGeH00uBrcBN4fRS4Iw48cRizwnHXwSu70eZvi8yfhKwP1wnF/jnsMzy4n0G4bRlwIe62VfmAJvD4SzgCoJ9ZgrwWYITiOowvh8CD3Z5Tz8Birt83rH3+RLw/bBsZhLsS+eF8+4E/g8YHsa6PBZHnBjfC2wCLBwfBjQDlWGcm4DKSFyTevhu/Ah4OCy7ncBfReZdB7QDHwGygTuAjcD3wvd/AcF+VRL5LjWG8eUD3wb+0OW7NDmy7B1dy7zL9/KObj6XPGADwYlILnBl+BnFXm8WwffojDDua8P9ID/VuahfeSvVAQzmB3B2uFOMDMffBP4hMn8ecH9k/GLgzci4Eya2cPxh4JZutjUXWBwZXw98NDIeSwRVkWk7gf8XGf8N8Llw+DngE5F5U8L3khN5rerI/NeAq8Lh2zmUaG8FHu1jWTlQEGfejcDqcPgnwL3RbYfTr46+/y7zovHEYj8i8fexTKOJ/wvAw5HxLKAOmBPvM+hDGcwBOoEGguqPJZEyXQmcH1l2bJzP44Q4n3cOQTLvAEoj878GzAuH1xGeUITjN9B94jeCBPzecPzvgOfD4ckESe99QG4v77UI2AvMDcd/CDwemX9d7DMPx6eH7yd64rST8CSC4LsUPWkqCd/zuMh3aSAS/3uBLYQHvnDagsjr/QD4cpfXW0V4QjVYHqrqOTbXAs+4+9vh+C/oUt0DbIsMNxHssL3ON7PRZvZQ+LN6L/AzgnrxqE1xYtoeGW6OMx7bfiXBmU3MBoIkMrofsUOQdNbGmd5VrIzGxpk3NjL/nwmSz2thFctH+7mdbvWxTKMOKyN37yQo86rIMvE+g55scfcKdx/u7jPd/aFw+gTg0bD6oIHgQNDB4Z9Hd9uqBHa5e2Nk2oZInJVd1o1+7ofxIJM9RHCgBfgQwS9V3H0N8DmCA+2OsCwru3mpKwjO6H8bjv8cuMjMRkWW6bpv4u7d7a9E34O77yM4eHa3/aNVCdSF5RATLa8JwE2xzyn8rMYlII6EUuI/SmZWCPw1cI6ZbTOzbQQ/D2eY2YwB2MRXCc5iprt7GfA3BAkxyo9Yq++2EOzEMeMJvqjb4y/erU1A13aLeFYBm4EPRieaWRbwVwS/QHD3be7+d+5eCXwc+L4Fl+ltIqiWOha9lWnX8jysjMJ2inEEZ/3drXO0NgEXhQeF2KPA3fuyrS3AcDMrjUwbH4lzaxh3dF5PHgSuDOvQzyD4pRgEELTPnE1QLg78RzevcS1B0t4Yfjd+RVB18qFett2Tg+/BzEoIqq629LJOvDLbT/CLJGZMZHgrUBV+1jHR8toEfKXL51Tk7g/26R2kCSX+ozeX4IxsKkGd6kzgFIK61IG4gqGUoNFzj5lVATcPwGtGPQj8g5kdH36Jvgr80t3b+/k6PydooP5rM8sxsxFmNrPrQuEZ1D8BnzezD5lZgZmNAe4nuMLlWwBm9kE7dHnnboIvbidBW8hYM/ucmeWbWan1/xLB3sp0O4cfXB4GLjGz880sl6B9oYXgp/9Auwf4SqzB0sxGmdnlfVnR3TeFMX0tLNca4GMEv2ggeB+3mtmwsGw/3cvrLSb4BXY/8LS7N4QxTTGz8yy4bPkAwRl5Z9f1w7I9H7iUQ9+NGQQHiWP5blxswYUEecCXgVfC996T7cAICy9cCC0JX2t4uA9+LjLvZYIToM+YWa6Z/SVBG13MfcCNZnaGBYrN7JIuB920p8R/9K4F/tuD66m3xR7A3cA1duyX3P0bcBqwB3iCoPF1IP0Y+ClBo+CfCb7IPSaEeNx9I0HbxU0cqreO+4vH3X8JfJjgl9FOYAVBQ+VZ7r4zXOwdwKtmtg+YD3zW3deF1Rh/AVxGUAW1Gji3n+H2VqZfIzgwNZjZP7n7KoJfBd8lSISXAZe5e2t3Gwirp67pZ1wQNFbOB54xs0aCht7+HNiuJqj33wI8CnzJ3X8Xzvs3guqKPwPPEHzuvfkFQV3+LyLT8gkait8m+AyOI2jj6erDwBJ3f6bLd+M7QI2ZTevH++oa05cI9rPTCT6bHrn7mwQnOevCz7WS4P0vJWijeQb4ZWT5VuAvCdogdgH/j8h+4u4LCdo97iY4MVkTLjuoxFruRUTSlpnNI2iA/XyqYxkKdMYvIpJhlPhFRDKMqnpERDKMzvhFRDLMoOjsaeTIkT5x4sRUhyEiMqgsWrTobXcf1XX6oEj8EydOZOHChakOQ0RkUDGzuP/SVlWPiEiGUeIXEckwSvwiIhlGiV9EJMMo8YuIZJhBcVXP0XhscR13Pb2KLQ3NVFYUcvOFU5g7q6r3FUVEhrghmfgfW1zHrY/U0tzWAUBdQzO3PlILoOQvIhlvSFb13PX0qoNJP6a5rYO7nl6VoohERNLHkEz8Wxqa+zVdRCSTDMnEX1lR2K/pIiKZZEgm/psvnEJhbvZh0wpzs7n5wikpikhEJH0MycbdWAPufz71Jlv2HKAoL5uvXjFdDbsiIgzRM34Ikv+CW8/nPSeOZOKIYiV9EZHQkE38MdOrynlreyMHulzlIyKSqYZ84q+prqC901mxdW+qQxERSQsZkPjLAVi2qSG1gYiIpIkhn/jHlhcwsiSfZXV7Uh2KiEhaGPKJ38yoqS5n2WYlfhERyIDED0F1z9r6fexraU91KCIiKZcxid8dlqu6R0QkUxJ/BQC1qu4REcmMxD+yJJ+qikKWbm5IdSgiIimXEYkfgj9y1aqqR0QkcYnfzMaZ2QtmtsLM3jCzz4bTbzezOjNbEj4uTlQMUTXjytmws4mGptZkbE5EJG0l8oy/HbjJ3acC7wI+aWZTw3nfcveZ4eO3CYzhoJqqCgBd1ikiGS9hid/dt7r7n8LhRmAlkLKe0qZXBf/gVXWPiGS6pNTxm9lEYBbwajjpU2a2zMx+bGbDkhFDeVEuE0cUsVRdN4hIhkt44jezEuA3wOfcfS/wA2ASMBPYCnyjm/VuMLOFZrawvr5+QGKpqa7QGb+IZLyEJn4zyyVI+j9390cA3H27u3e4eydwH/DOeOu6+73uPtvdZ48aNWpA4qmpLmfrngPsaDwwIK8nIjIYJfKqHgN+BKx0929Gpo+NLHYFsDxRMXSlP3KJiCT21otnAR8Gas1sSTjtX4GrzWwm4MB64OMJjOEwp1aWkWWwdPMezj9ldLI2KyKSVhKW+N39D4DFmZWUyzfjKc7PYfJxJdTqH7wiksEy5p+7MTXVFSzbvAd3T3UoIiIpkYGJv5yd+1upa2hOdSgiIimRgYm/AlADr4hkroxL/CePKSUny1iqxC8iGSrjEn9BbjYnjy2ltq4h1aGIiKRExiV+gOlVQQNvZ6caeEUk82Rk4p9RXU7jgXY27GpKdSgiIkmXkYk/1sC7TNfzi0gGysjEf+LoEvJzstQ3v4hkpIxM/LnZWZxaWaYzfhHJSBmZ+CGo7llet5f2js5UhyIiklQZnPjLaW7rYG39/lSHIiKSVBmd+AGWqrpHRDJMxib+E0aWUJyXra4bRCTjZGziz8oyplWVq4FXRDJOxiZ+gBnjKli5tZHWdjXwikjmyOjEX1NdTmtHJ6u2NaY6FBGRpMnsxF9VAcAyddgmIhkkoxP/uOGFVBTlsmyTGnhFJHNkdOI3M6ZXleuSThHJKBmd+AFmVFewesc+mls7Uh2KiEhSZHzin15dTkens2KrqntEJDNkfOKfcbCLZiV+EckMGZ/4R5flM6o0X4lfRDJGxid+M2NGtf7BKyKZI+MTPwRdNK97ez+NB9pSHYqISMIp8RM08LrD8rq9qQ5FRCThlPiBmqqgi2ZV94hIJlDiB0aU5FNVUagGXhHJCEr8oRnjytVnj4hkhIQlfjMbZ2YvmNkKM3vDzD4bTh9uZs+a2erweViiYuiP6VUVbNrVzK79rakORUQkoRJ5xt8O3OTuU4F3AZ80s6nALcBz7n4i8Fw4nnIzwlsx1tapukdEhraEJX533+rufwqHG4GVQBVwOfBAuNgDwNxExdAf08LEv2xTQ2oDERFJsKTU8ZvZRGAW8Cow2t23hrO2AaO7WecGM1toZgvr6+sTHmNZQS4njCxmmc74RWSIS3jiN7MS4DfA59z9sAvl3d0Bj7eeu9/r7rPdffaoUaMSHSYQ3JFLl3SKyFCX0MRvZrkESf/n7v5IOHm7mY0N548FdiQyhv6YXl3B9r0tbN97INWhiIgkTCKv6jHgR8BKd/9mZNZ84Npw+Frg8UTF0F+xBl5dzy8iQ1kiz/jPAj4MnGdmS8LHxcCdwF+Y2WrgfeF4WphaWUaW6R+8IjK05STqhd39D4B1M/v8RG33WBTl5XDS6FKd8YvIkKZ/7nYxvSpo4A3anUVEhh4l/i5qxlWwu6mNzbubUx2KiEhCKPF3oQZeERnqlPi7mDKmlNxsU4dtIjJkKfF3kZ+TzSljy1i2SWf8IjI0KfHHMb2qnOV1e+jsVAOviAw9SvxxzKiuoLGlnT/v3J/qUEREBpwSfxzTq3UrRhEZupT44zjxuBIKcrN0ZY+IDElK/HHkZGdxamW5Er+IDElK/N2oqS7njS17aO/oTHUoIiIDSom/GzOqKzjQ1snqHftSHYqIyIBS4u9GrIG3VtU9IjLEKPF34/gRxZTm57BUV/aIyBDT526ZzWwYUAk0A+vdfUhXfmdlGdOqyqnVPXhFZIjp8YzfzMrN7F/NrBZ4Bfgh8DCwwcx+ZWbnJiPIVKkZV87KrXtpae9IdSgiIgOmtzP+XwM/Ad7j7g3RGWZ2OvBhMzvB3X+UoPhSqqaqgrYOZ9W2RmqqK1IdjojIgOgt8V/v7hvizXD3RcCigQ8pfdSEDbxLN+9R4heRIaO3xt1HkxJFmqoeVsiwolyWbWpIdSgiIgOmt8Tf3T1zM4KZUVNdoQZeERlSeqvqqTKz73Q3090/M8DxpJ0Z1eXc/UI9Ta3tFOUl7N70IiJJ01sma2aI1+P3Znp1BZ0OK7bsZfbE4akOR0TkmPWW+He6+wNJiSRNRRt4lfhFZCjorY6/NSlRpLHRZQWMLsunVv/gFZEhorfEf1VPMy1QPYDxpKWa6gp10SwiQ0ZvVT13mVkW8DhBXX89UABMBs4Fzge+BGxOZJCpVlNVzrMrtrP3QBtlBbmpDkdE5Jj0mPjd/YNmNhW4BvgoMBZoAlYCvwW+4u4HEh5litWMqwBg+eY9nDl5ZGqDERE5Rr1en+juK4DbkhBL2ppedaiBV4lfRAY7dcvcB8OL8xg3vJDauoZUhyIicsyU+PuoprqCpZvUwCsig1/CEr+Z/djMdpjZ8si0282szsyWhI+LE7X9gVZTVU5dQzM797WkOhQRkWPSYx2/mZ3W03x3/1MPs+cBdxN06xz1LXf/ep+iSyOx3jmX1e3h3CnHpTYYEZFj0Fvj7jfC5wJgNrCUoOO2GmAh8O7uVnT3l8xs4gDEmBamVZVhFtyDV4lfRAazHqt63P1cdz8X2Aqc5u6z3f10YBZQd5Tb/JSZLQurgoZ1t5CZ3WBmC81sYX19/VFuauCUFuRywshilukfvCIyyPW1jn+Ku9fGRtx9OXDKUWzvB8AkYCbBweQb3S3o7veGB5rZo0aNOopNDTz9g1dEhoK+Jv5lZna/mc0JH/cBy/q7MXff7u4d4Y3a7wPe2d/XSKWa6nJ2NLawbc+Q/8+aiAxhfU38HwHeAD4bPlaE0/rFzMZGRq8Alne3bDo61FNnQ2oDERE5Bn26s4i7HzCz7wG/AxxY5e5tPa1jZg8Cc4CRZraZoE+fOWY2M3yN9cDHjzryFJg6tpzsLKN28x4uPHVMqsMRETkqfUr8ZjYHeIAgWRswzsyudfeXulvH3a+OM/lH/Q8xfRTmZXPS6FKd8YvIoNbXewl+A7jA3VcBmNlJwIPA6YkKLF3VVJXz9IptuDtmGX1LYhEZpPpax58bS/oA7v4WkJH9E3e609DUxvG3/paz7nyexxYf7VWtIiKp0dcz/oVmdj/ws3D8GoI/cGWUxxbXMX/ploPjdQ3N3PpIcJXr3FlVqQpLRKRf+nrG//cEV/J8JnysCKdllLueXkVLe+dh05rbOrjr6VXdrCEikn76elVPC/DN8JGxtjQ092u6iEg66q2Ttofd/a/NrJbgEszDuHtNwiJLQ5UVhdTFSfKVFYUpiEZE5Oj0dsb/2fD50kQHMhjcfOEUbn2klua2joPT8rKzuPnCKSmMSkSkf3q75+5WM8sG5oWdtWW0WAPuXU+voq6hmewsY1hxLhdPH9vLmiIi6aMv99ztMLNOMyt394zvoWzurKqDB4AX3tzBR+a9zn3/t45Pnjs5xZGJiPRNXy/n3AfUmtmzwP7YRHf/TEKiGiTOPfk4Lpo2hu88t5rLaioZP6Io1SGJiPSqr5dzPgJ8AXgJWBR5ZLwvXjaVnCzjS/OX435E+7eISNrp6+WcD5hZITA++g9egbHlhfzDX5zEHU+s5Ok3tvH+aarvF5H01qczfjO7DFgCPBWOzzSz+QmMa1C57syJnDK2jNvnr2BfS3uqwxER6VFfq3puJ7hpSgOAuy8BTkhIRINQTnYWX7liGtsbD/CtZ99KdTgiIj3qa+Jvi3NFT2fcJTPUaeOHcfU7xzNvwXre2JLxFz+JSBrra+J/w8w+BGSb2Ylm9l1gQQLjGpT+5cKTqSjM5fOPLaezUw29IpKe+pr4Pw2cCrQAvwD2AJ9LUEyDVnlRLrddcgqLNzbw0OubUh2OiEhcvfXVUwDcCEwGaoF3u7taL3twxawqHl64iTufXMkFp45mZEl+qkMSETlMb2f8DwCzCZL+RcDXEx7RIGdm3DF3Os1tHXz1iZWpDkdE5Ai9Jf6p7v437v5D4ErgvUmIadCbfFwJN54ziUcW17Fg7dupDkdE5DC9Jf622ICqePrnk+dOZvzwIj7/2HJa2jt6X0FEJEl6S/wzzGxv+GgEamLDZrY3GQEOVgW52fz75aeyrn4/9720LtXhiIgc1GPid/dsdy8LH6XunhMZLktWkIPVnCnHccn0sXz3+TVs3NmU6nBERIC+X84pR+kLlwaduH3hcXXiJiLpQYk/wcaUF3DTBVP4/Vv1PLl8W6rDERFR4k+Gv333BE6tLOPf/ucNdeImIimnxJ8EQSdu09nR2MI3n1EnbiKSWkr8STJzXAXXnDGeeQv+zPI6deImIqmjxJ9EN194MsOL87jtseV0qBM3EUmRhCV+M/uxme0ws+WRacPN7FkzWx0+D0vU9tNReWEun79kKks3NfDgaxtTHY6IZKhEnvHPA97fZdotwHPufiLwXDieUS6fWcmZk0bwH0+9SX1jS6rDEZEMlLDE7+4vAbu6TL6coOM3wue5idp+ujIzvjx3Gi1tnXz1t+rETUSSL9l1/KPdfWs4vA0YneTtp4VJo0q48ZwTeHRxHQvWqBM3EUmulDXuevA31m5bOM3sBjNbaGYL6+vrkxhZcnzi3MlMGFHE5x9XJ24iklzJTvzbzWwsQPi8o7sF3f1ed5/t7rNHjRqVtACTJejEbRrr6vdz7+/ViZuIJE+yE/984Npw+Frg8SRvP62cc9IoLqkZy3dfWMOGnftTHY6IZIgeb714LMzsQWAOMNLMNgNfAu4EHjazjwEbgL9O1PYHiy9eOpXfr6rn4z9dxN7mNrbuOUBlRSE3XziFubOqUh2eiAxBCUv87n51N7POT9Q2B6PRZQVcMHU0jyyuOzitrqGZWx+pBVDyF5EBp3/upoFX1u08YlpzWwd3Pb0qBdGIyFCnxJ8Gtu45EHf6lobmJEciIplAiT8NVFYU9mu6iMixUOJPAzdfOIXC3Owjpp81eUQKohGRoS5hjbvSd7EG3LueXsWWhmbGlBVQUZzLwws3k5eTxRcvPZW8HB2jRWRgKPGnibmzqg67gqej07nr6VXc8/u1vLV9H9+/5jRGluSnMEIRGSp0GpmmsrOMWy46mW9fNZNlmxv4wHf/oBu4iMiAUOJPc5fPrOLXN54JwJX3LGD+0i0pjkhEBjsl/kFgWlU58z99NjVVFXzmwcX8x1Nv6g5eInLUlPgHiZEl+fzs+jO45ozx/ODFtVz/wOvsPdCW6rBEZBBS4h9E8nKy+MoV07lj7jT+b/XbzP3eH1lbvy/VYYnIIKPEPwj9zbsm8Iu/exd7mtqYe/cfeeHNbnu3FhE5ghL/IPXO44cz/9NnM2FkER994HW+/+IagnvbiIj0TIl/EKuqKORXHz+Ty2oq+c+nVvHpBxfT3Kq7eYlIz5T4B7nCvGy+fdVMbrnoZJ6o3cqV9yxg8+6mVIclImlMiX8IMDNuPGcSP772HWzc1cQH7v4jr8bp6llEBMAGQ73w7NmzfeHChakOY1BYV7+P63+ykI07m5g7q4qX1+5kS0Oz7uolkoHMbJG7z+46XWf8Q8wJo0p47JNncdLoUn69aDN1Dc04h+7q9VjkTl8ikpmU+IegsoJcGppaj5iuu3qJCCjxD1nd3dWrTnf1Esl4SvxDVE9377rq3pdZsOZtXfcvkqGU+IeoeHf1KsjJYu7MStbV7+dD97/Klfe8zIurdugAIJJhdCOWIarrXb2iV/UcaOvgV4s2c8+La7nuv1+nprqcT507mfedMpqsLEtx5CKSaLqcM4O1tnfy6OLNfO+FtWzc1cTJY0r51HmTuWjaWLJ1ABAZ9Lq7nFOJX2jv6OR/lm3h7ufXsLZ+P5NGFfOp8yZzWU0lOdmqDRQZrJT4pVcdnc6Ty7dy9/NreHNbIxNGFPGJOZO4Yla1bvYuMggp8UufdXY6v1u5ne8+v4bauj1UVRRy4zkn8MHZ4yjo0mAsIulLiV/6zd35/Vv1fPf5NSzasJvRZfnc8N5JlObn8O3nVqsrCJE0113i11U90i0zY86U4zjnpFG8vHYn33l+NV/+3xWHLRPrCgJQ8hcZJFRxK70yM86cPJKHbng3I0vyjpjf3NbBfzz1ZgoiE5GjkZLEb2brzazWzJaYmepwBpGd+47sAwiCLiKuf+B1frNoM3uadBN4kXSWyqqec9397RRuX45CZUVh3P5+ivOyWbFlL79buYOcrOAXwkXTxnDB1NGMKMlPQaQi0h3V8Uu/3HzhFG59pJbmtkO3eCzMzeYrV0zn8pmVLN28hyeXb+Wp5du49ZFabnu0lnceP5yLpo3l/dPGMLqsIIXRiwik6KoeM/szsBtw4Ifufm+cZW4AbgAYP3786Rs2bEhukNKtxxbXxe0KIsrdWbm1kaeWb+XJ5dtYvWMfAKdPGMZF08Zw4aljGDe8KBXhi2SMtLqc08yq3L3OzI4DngU+7e4vdbe8Lucc/NbsaOTJ2m08uXwbK7buBWB6VTnvnzaGi6aNCW4g04cDioj0XVol/sMCMLsd2OfuX+9uGSX+oWXDzv08tTw4CCzZ1ADA2PIC6htbaO88tD8W5mbztb+cruQvcpTS5taLZlZsZqWxYeACYHmy45DUmTCimI+fM4nHPnkWC245jy9eOpWd+1oPS/oQXCb6lSdW0t7RmaJIRYampJ/xm9kJwKPhaA7wC3f/Sk/r6Ix/6Dv+lifobk8szstm1vhhzJ44jNkThjNrfAXF+bouQaQ3afPPXXdfB8xI9nYlvXV3meiwolwuralk4YbdfPu51bhDdpZxythSZk8YfvBgMKZcVwuJ9JVOmyQtdHeZ6JcuO/VgHf/eA20s3tjAovW7eH39bn75+ibmLVgPQFVFIe+YOIzTJw5n9oRhnDS69OA9BdRoLHK4lDfu9oWqejJDfxN0W0cnK7fu5fX1u1m0ITgY1De2AFBakMNp44dRnJfN71buoDXSTqBGY8kUaXtVT18o8UtfuDubdjWzMDwILNqwi7e274u77PDiPB77xFlUDyvU7SZlyFLil4w08ZYnepxflJfNiaNLmTK6hCljypgyupQpY0oZWZKHmQ4IMrilTeOuSDJVddNoPKokn5suOIk3tzXy1vZGnlu5g4cXbj44f3hx3sGDQOxx0uhSSsKridRuIIOZEr8Mad01Gt92ySlHJOr6xhbe2t4YHAy2NfLm9kYeXriJptZD61YPK6S8MJdV2xoP/u9A9ySQwUZVPTLkHcvZeWens3l3M6u2N7Jq217e3NbIU8u3HfFnM4C87CwuqRnL+OFFjB9exIQRRYwfUcSoknxVG0lKqI5fZID09GezqopCtuxpJvq1KszNZvzwIsaFB4MJI8Lh4UVUDys67Eb2qkKSgaQ6fpEB0t2fzaoqCvnjLefR0t7B5t3NbNzVxMadTWzc1cSGnU1s3LWfP6yp50DboUtLzaCyvJDxw4twdxZu2H1YFdItv1kGqApJBpYSv0g/ddducPOFUwDIz8lm0qgSJo0qOWJdd6e+seXgwWDDriY27Wpiw879LNnUQNcapAPtndz08FIeeHk9Y8sLGFNWyJjyfMaUF4bjBRxXlk9+TnaPMeuXhEQp8Yv0UyxhHk0iNTOOKyvguLICZk8cfti847u59LTDneK8HFZta+TFVfWHNTbHjCzJY0x4IBhTXsDY8kJGlxUwtryA5XV7+Nazb3GgPfilocZoUeIXOQpzZ1UNeNLsqQrpZ9efAQS/GBpb2tm+5wBb9xxg254DbNsbG25m8+5mFm3Yze5e7nvc3NbBbY/WsnFXEyNL8hlZksfI0nxGFuczsjSPory+pQb9khiclPhF0kRvVUgQ/GIoK8ilrCCXE0eXdvtaB9o6Dh4Urrr3lbjL7G/t4JvPvhV3XlFeNiNK8sKDQnhgCIdj05duakjILwkdTBJPiV8kTRxLFVJXBbnZTBxZzMSRxd3+ia2qopDn/+kcdu1v5e3GVt7e1xI+Dg3v3NfKpl1NLN64m137W49og+iqua2DWx5Zxstrd1JRnMuwojyGFeVSXhg8DyvOo6Iol4rCvMOuZop5bHHdYQc/VUslhi7nFBniuiZTOLqO6jo6nd1NrQcPCNfc/2q3yx5Xmk9DU9thneN1VZKfExwEioIDREVRHs+v3M7+OG0Yx5XmM/9TZ1NWmENhbvZR/S8iE39J6HJOkQw1UL8ksrPsYHUPdN8dRuyyVnenqbWD3U2tNDS1sbupld1NbTRExqPTN+5qipv0AXY0tvCurz0HQE6WUVaYS1lBTvicS1lhTvgcf/rr63fx7d+tTkgD92A8oOiMX0SOykD9kog6687nu70hz80XnszeA23sbW4Ln9sj4+0Hp0f/J9GbnCxjxrgKSvJzDj0KcijOz6E0znB0ueL8HH5bu3XAyyBmIA4oOuMXkQE1kG0SMX25IU9vWto7aDx4IAie//bHr8Vdtr3TKczNpqG5jc27m9jf0sG+lnb2tbQf9XuItXP83+q3Kc7Ppigvh+K8bArzsinOz6EoL5vivByK8oPn2DJFecFzXk5Wwts6lPhF5KgN9GWtA3Ewyc/JJr8k+2CVFPRcLRW7VDaqs9Npautg34F29rW00XigPTwoxIaDg8PXn4l/VdSBtk5eWbeTptZ29rd20Nre918hudlGe6fTtTKmua2Du55epcQvIkNPIv4j0ZdLZaOysuxglQ50fz/nB1/b1GM7R0x7RydNbR00tXSwv7X90HNrcEBpam2nqbWDptYO9re08/0X18bd3pY42zoaSvwiMuQloloK+n5AycnOoiw7i7KC3D697uNLtsQ9oFRWFB5TvAfjGZBXERFJc4n4JZHqA8rRUuIXETkGg+mAEqPELyKShhJxQIk58j/TIiIypCnxi4hkGCV+EZEMo8QvIpJhlPhFRDLMoOikzczqgQ2pjgMYCbyd6iB6oRiPXbrHB4pxoAz1GCe4+6iuEwdF4k8XZrYwXk936UQxHrt0jw8U40DJ1BhV1SMikmGU+EVEMowSf//cm+oA+kAxHrt0jw8U40DJyBhVxy8ikmF0xi8ikmGU+EVEMowSf8jM3m9mq8xsjZndEmf+P5rZCjNbZmbPmdmEyLwOM1sSPuanMMbrzKw+Esv1kXnXmtnq8HFtiuL7ViS2t8ysITIvWWX4YzPbYWbLu5lvZvad8D0sM7PTIvOSUYa9xXdNGFetmS0wsxmReevD6UvMbGEi4utjjHPMbE/k8/xiZF6P+0gSY7w5Et/ycP8bHs5LeDma2TgzeyHMKW+Y2WfjLJO4fdHdM/4BZANrgROAPGApMLXLMucCReHw3wO/jMzblyYxXgfcHWfd4cC68HlYODws2fF1Wf7TwI+TWYbhdt4LnAYs72b+xcCTgAHvAl5NVhn2Mb4zY9sFLorFF46vB0amQRnOAf73WPeRRMbYZdnLgOeTWY7AWOC0cLgUeCvO9zlh+6LO+APvBNa4+zp3bwUeAi6PLuDuL7h7Uzj6ClCdbjH24ELgWXff5e67gWeB96c4vquBBwc4hl65+0vArh4WuRz4iQdeASrMbCzJKcNe43P3BeH2ITX7YV/KsDvHsg/3Sz9jTPq+6O5b3f1P4XAjsBLo2vl+wvZFJf5AFbApMr6ZIz+EqI8RHIljCsxsoZm9YmZzExAf9D3Gvwp/Fv7azMb1c91kxEdYTXY88HxkcjLKsC+6ex/JKMP+6rofOvCMmS0ysxtSFFPMu81sqZk9aWanhtPSrgzNrIggaf4mMjmp5WhmE4FZwKtdZiVsX9QduPrJzP4GmA2cE5k8wd3rzOwE4Hkzq3X3tSkI73+AB929xcw+DjwAnJeCOHpzFfBrd++ITEuXMhwUzOxcgsR/dmTy2WEZHgc8a2Zvhme+yfYngs9zn5ldDDwGnJiCOPriMuCP7h79dZC0cjSzEoKDzufcfW8ithGPzvgDdcC4yHh1OO0wZvY+4DbgA+7eEpvu7nXh8zrgRYKjd9JjdPedkbjuB07v67rJiC/iKrr8tE5SGfZFd+8jGWXYJ2ZWQ/D5Xu7uO2PTI2W4A3iUoGol6dx9r7vvC4d/C+Sa2UjSqAwjetoXE1qOZpZLkPR/7u6PxFkkcftiIhswBsuD4JfPOoLqh1ij06ldlplF0DB1Ypfpw4D8cHgksJoENFj1McaxkeErgFf8UGPQn8NYh4XDw5MdX7jcyQSNZ5bsMoxsbyLdN0xewuENaq8lqwz7GN94YA1wZpfpxUBpZHgB8P4UleGY2OdLkDQ3huXZp30kGTGG88sJ2gGKk12OYXn8BPivHpZJ2L6YkAIfjA+CFvS3CJL7beG0fyc4uwf4HbAdWBI+5ofTzwRqw524FvhYCmP8GvBGGMsLwMmRdT8aJow1wEdSEV84fjtwZ5f1klmGDwJbgTaCutGPATcCN4bzDfhe+B5qgdlJLsPe4rsf2B3ZDxeG008Iy29puA/clsIy/FRkP3yFyEEq3j6SihjDZa4DHuqyXlLKkaCKzoFlkc/y4mTti+qyQUQkw6iOX0Qkwyjxi4hkGCV+EZEMo8QvIpJhlPhFRDKMEr8MamY218zczE4ewNecY2b/Gw5/INaLZLitqUfxei+aWb9ulm1m/2Vm7+1lma+bWTr+M1vSnBK/DHZXA38Inwecu8939zvD0blAvxN/f5nZCOBd3ns3Ad8FEta1sQxdSvwyaIX9nJxN8OecqyLT55jZ783scTNbZ2Z3hv3Yvxb2sz4pXG6emd0Tdg73lpldGmcb15nZ3WZ2JvAB4K6wn/ZJ0TN5MxtpZuvD4UIze8jMVprZo0Bh5PUuMLOXzexPZvar8D109VfAU5F1vmhmr4f9xt9rZgbg7huAEWY25ljLUjKLEr8MZpcDT7n7W8BOMzs9Mm8Gwb8gTwE+DJzk7u8k+OfrpyPLTSToVuAS4B4zK4i3IXdfAMwHbnb3md5zB3J/DzS5+ynAlwj7TAr7q/k88D53Pw1YCPxjnPXPAhZFxu9293e4+zSCg0j0APWncHmRPlPil8HsaoI+3Qmfo9U9r3vQ53kLwV/enwmn1xIk+5iH3b3T3VcT9CMzEG0F7wV+BuDuywj+lg9BfytTgT+a2RLgWmBCnPXHAvWR8XPN7FUzqyXobfXUyLwdQOUAxCwZRN0yy6AU3ibvPGC6mTnB3Z3czG4OF2mJLN4ZGe/k8P2+a58l/enDpJ1DJ09xfyl0DZvgBhq9tUc0x14v/AXyfYJ+WjaZ2e1dtlUQLi/SZzrjl8HqSuCn7j7B3Se6+ziCXgrf08/X+aCZZYX1/icAq3pYtpHgNnkx6znU9fWVkekvAR8CMLNpQE04/RXgLDObHM4rNrOT4mxnJTA5HI4l+bfD9oAruyx7EhD3vrIi3VHil8HqaoK+0qN+Q/+v7tkIvEbQ/e2N7n6gh2UfAm42s8XhgeLrwN+b2WKC7qRjfgCUmNlKgt5JFwG4ez1Bj5APmtky4GXiVy09QXDfWty9AbiPILk/DbweWyjsz30yQVuBSJ+pd07JWGY2j+Cm4L9OdSxdmdkfgEvDxN/dMlcQ3LD7C0kLTIYEnfGLpKebCG660pMc4BtJiEWGGJ3xi4hkGJ3xi4hkGCV+EZEMo8QvIpJhlPhFRDKMEr+ISIb5/643DCd2LwQDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "x_a = linspace(0,2,16)\n",
    "y_period = []\n",
    "for i in range(0,16):\n",
    "    y_period.append(period(x_a[i]))\n",
    "\n",
    "plt.scatter(x_a,y_period)\n",
    "plt.plot(x_a,y_period)\n",
    "plt.title(\"Anharmonic Oscillator: Period vs Amplitude\")\n",
    "plt.xlabel(\"Amplitude (a)\")\n",
    "plt.ylabel(\"Period (T)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "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",
   "metadata": {},
   "source": [
    "1. Here I've defined a function called `function4` to stand in for the mathematical function $1 + \\frac{1}{2}\\tanh 2x$. I then defined a function called `centralDiff` which takes in the arguments `x`, `h`, and `function`. The variable `x` is the point upon which we will be taking the central difference about `h` for a given function `function`. The `centralDiff` function calculates the central difference using the equation $$\\frac{df}{dx} = \\frac{f \\left(x + \\frac{h}{2} \\right) - f \\left( x - \\frac{h}{2} \\right)}{h}$$ \n",
    "    and then returns it.\n",
    "    \n",
    "    I set my `h` close to zero with a value of $0.00001$. I then create a set of $100$ equally spaced values from $-2$ to $2$ under the list `x_set`. I then calculate the derivative of the function above<sub>[1]</sub> for ever value in `x_set` and append it to the list `y_cd`.\n",
    "    \n",
    "    <sub>[1] To calculate the derivative, I pass on each value of `x_set` with `h` and `function4` to the `centralDiff` function.</sub>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "Jiu5YAx6aMmq"
   },
   "outputs": [],
   "source": [
    "from numpy import tanh\n",
    "\n",
    "def function4(x):\n",
    "    y = 1 + 0.5*tanh(2*x)\n",
    "    return y\n",
    "\n",
    "def centralDiff(x,h,function):\n",
    "    return (function(x+(h*0.5)) - function(x-(h/2)))/h\n",
    "\n",
    "h = 0.00001\n",
    "x_set = linspace(-2,2,100)\n",
    "y_cd = []\n",
    "for i in range(0,len(x_set)):\n",
    "    y_cd.append(centralDiff(x_set[i],h,function4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Here, I've defined a function called `analfunction4` to stand in as the analytic formula of the derivative of the mathematical function above, i.e. $1 - \\tanh^2(2x)$. I then pass on each value of `x_set` to `analfunction4` and append each returned result to the list `y_theo`. I then graphed the  numerical answers calculated using the central diff as red dots, and the analytically computed answers as the blue line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEWCAYAAAAkUJMMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABEYklEQVR4nO3dd3wUZf7A8c83hQSQDiI9NE86UhUEAkJo51KEA0QUC+oBevbDO0/96Xmeeh52z4oFpCqwUiQghI4QpIMKCEhQkC6YkPr8/pgJLiGbRjaz5ft+vfLK7uzszHeeKd+dZ555RowxKKWUUk4JczoApZRSoU0TkVJKKUdpIlJKKeUoTURKKaUcpYlIKaWUozQRKaWUclRIJSIRGSki8U7HkU1ESovIFyJyWkRmOh1PUYjIfhHpeQnfv2CdiEhnEdktImdFZKCIVBeRFSJyRkReKp6og4NdRg18MN0uIvJdcU83EIhIgojcWcTv/k1E3rvE+ceIiBGRiEuZTkkqjpiLlIhE5CYRSbR3hJ9FZKGIXFfUIEqKMWaKMSbO6Tg8DAGqA1WMMUNzG0FErhSRmSJyzE5YW0XkQREJv5QZl8QGLyIfikianUTOiMh2EXlORCpkj5PLOnkaeN0Yc5kxZg5wF3AMKG+MechXsfpSQcqhKOwy+qEY4jMi0shjuiuNMX+41On6kohcZh9/Fjo0/1gRSfIcZoz5lzGmSEmsEPPdLyIp9rJn/9X05TxLQqETkYg8CLwM/AvrIFoXeBMYUKyRFTM//YVRD/jeGJOR24ci0hD4GjgItDDGVACGAu2Acr4OrpjK7AVjTDmgGnAbcA2wWkTKehm/HrAjx/udpgh3XvvZOi9sOXjlZ8vllBuBVKCXiFzhdDAl7Ab7R0j230+F+bJfbj/GmAL/ARWAs8DQPMaJwkpUP9l/LwNR9mexQBLwKPAL8DMwEOgHfA+cAP7mMa2ngFnAdOAM8A3QyuPzCcBe+7OdwCCPz0YDq4GJwHHgn/awVfbnYn/2C/ArsA1o7rGcHwNHgQPA40CYx3RXAf8BTgL7gL55lEcTIAE4hXWAddnD/w9IA9LtMr0jl+9OBubns06uAdbY098CxHp8lgA8Y5fDGSAeqGp/9iNg7HmfBa71UmYNgaX2+2PAFKCixzz2Az29xPYh8M8cw8rZ6328Z3nar/cCWUCKHdNUu3zS7Pc9sX48Za/348AMoLL9/Rh7me6wl2+FPfx2YJe9vhYB9TziMcA9wG67DN8AxOPzMfZ3s7exNvbwmsBn9jayD7gvj3WUbzkUMM5xdpz7PIY1AjoCh4Fwj/EHAVvt1x2Atfby/Qy8DpSyP1thT+c3u4yHYe+n9ud/BWbliP0V4FWPfeV9e7qH7G0mPJcyqGmv18oew67G2qYi7eVYDpy2h03PZ7tfCjyLdUx4OMdn+4GHga329KYD0fZnlYB59no7ab+unWOfuRMohXU8auHx2eVAMtaPoxSsbTV7/6mJdbya7DH+dfy+bx4ERtvD+wObsI47B4GnPL4TY6+PCC/LvZ9c9jcKdtz9q72dfGLHOhPrGHMG6/h3JfAY1jHxIBDnbb6ey5ozZqwfWtn7zA/A3XmtS2NMoRNRHyDDWyHZ4zwNrLNXWjV7RTzjUSAZwBNYG98Ye4P4FGvHbGav4PoeC5uOVYUVibVx7QMi7c+H2htAGNYO9BtQw+MAlwHcC0QApbnwoNcb2AhUxEpKTTy++zEw144pBitJ3uEx3XQ79nDgz/aKl1zKIhLYA/wNa8PuYa+cP+RcmV7K8jBwWx6f18I6GPezy6CX/b6ax061197AStvv/+1tg/dSZo3s6UbZ63MF8HJ+O4a3A7BH+U73mOeqPDb4C6YB/AVr+6ptx/Q2MDXHMn0MlLXjH2Cvgyb2Mj0OrMlxgJ9nbwd1sbbHPh7b1yGgvb2NNMI6CIVhbTtP2Ou1AdYO1/sSyqEgcS4GKgOlPYY1sl/vBXp5jD8TmGC/bov1gyXCLqNdwP05pt3I430svyeielgH33L2+3CspHON/X62vQ7KYu3z6/Fy4MFKHmM83r8I/M9+PRX4u1220cB1eWz39bCSQFPgIeyEm2MbWo91bKhsL+899mdVsM6mymDt3zOBOR7fTQDutF+/CTyfY9v7ImcZeTk418Pa10dgHQeqAK09vtvCXtaWwBFgoLf9Mpdlyy0RFeS4+zzWPlPajvUc1nEwAmtb3Gevg+xj87489kvPZb0gZqxE2xBrn+mGtf208bY+jSl8IhoJHM5nnL1AP4/3vYH9HgWSgv2Lyd4QDNDRY/yNHivlKWCdx2dhWDtBFy/z3gwM8DjA/Zjj89H8noh6YCWYa7DPdjx2tDSgqcewu4EEj2ns8fisjL0MV+QSTxesZOI5/anYv4DIPxGlYx8UvXz+V+CTHMMWAbd67FSPe3w2FvjS2wafW5nlMs+BwKb8dgz7sw/J/QD8b2BxznXiZYO/YBpYB5XrPd7XsMsp+yBrgAYeny/E42zT3oaSsc827PGv8/h8Br8fwBcBf8kl/o65bFuPAZMuoRwKEmePHN/3TET/BD7w2K9+w+OMKsf37gdm5zYd8/t+muTxfhVwi/26F7DXfl0dq3qstMe4I4BlXuZ7J7DUfi1Yv7q72u8/Bt7B4+wkj23wcWCz/boWkAlcnWMbutnj/QvYCS+XabUGTnq8T+D3RNQR68xa7PeJwJ9yKyN72FP8fnB+zLOM81mel4GJxst+mWPc/VhnYKfsvzn28PyOu2nYZ4UesS72eH+DPd2cx+aKXvZLz2XNL+Y55LIfef4V9hrRcaBqPnWMNbGqs7IdsIedn4YxJtN+nWL/P+LxeQpwmcf7g9kvjDFZWKeYNQFE5BYR2Swip0TkFNAcqJrbd3MyxizFqqJ4A/hFRN4RkfL29yNzWYZaHu8Pe0wn2X7pGXO2msBBO25v08rLcawDrTf1gKHZy2+XwXU5vnPY43Wylzg9XVBmdqu1aSJySER+xTqVr5r7VwusFla1R1HUA2Z7LO8urANRdY9xDuYY/xWP8U9gHQRzXZ9cWEZ1sHbw3GKomaPc/5YjhoLwLIeCxOl1e8aqVRgsIlHAYOAbY8wBON/gZZ6IHLbX4b8o3Dr8FCvBANxkv8+OORL42SPut7F+lefmM+BaEakBdMU6q1lpf/Yo1vKuF5EdInJ7HvHcglVFjDHmEFaV3q05xsl1nYpIGRF5W0QO2GWxAqiYW+MfY8zX9ndjReQqrDNidx5xefK27SAiHUVkmYgcFZHTWFXDhVkfA40xFe2/gfaw/I67R40x53JMJ+dx91gux+b8jhcXEZG+IrJORE7Y20Q/8lm+wiaitVi/gAbmMc5PWBtotrr2sKKqk/1CRMKwqmR+EpF6wLvAeKxWZxWB7VgbczaT14SNMa8aY9pineJfCTyCVT+dnssyHCpC7D8Bdey4izKtJVjVCN4cxDojqujxV9YY8+8CTNtb2eQc/i97WAtjTHngZi4s40IRkcuwrvWszG9cLw5iXZPzXOZo+4CUzeQY/+4c45c2xqwp4Lwaehm+L8c0yxlj+hV0IXIph4LE6XV7NsbsxDr49OXCZAHwFvAt0Nheh3+jcOtwJtbBuDbWtafsaR/EOh5U9Yi5vDGmmZcYT2JdpxxmxzjN2D+ZjTGHjTFjjDE1sWog3vRsyZdNRDoBjYHH7MR6GOvM5aYCXoR/CPgDVi1MeayECN7L4yOsbX4U1rWy7IN5nscWvG87YJWfG6hjrAZI/8tj/gWV33E3v3jz8xtW7U+2XBuI2D+EPsO6hl7dPi4vIJ/lK1QiMsacxqoXf0OsezzKiEiknQFfsEebCjwuItVEpKo9/uTCzCeHtiIy2N7I7sfa8Ndh1UkbrDp9ROQ2rDOiAhGR9vYvk0isQj4HZNm/CGYAz4pIOTvhPVjEZcj+RfWoXU6xWKfA0wr4/SeBTiLyYnbLIBFpJCKTRaSiHdMNItJbRMJFJNpuVlq7ANM+ivWLNL/7UMphnbKfFpFaWMm60EQkSkTaYp2mnwQmFWU6WDvts/Z6wd7OBuQz/mMi0swev4KI5NpUPhfvAQ+LSFuxNLLnux44IyJ/FetesHARaS4i7fObYB7lcClxZvsU6zpGV6zkka0c1oXxs/Yv+z/n+N4R8tgOjDFHsaqsJmEl4F328J+xEstLIlJeRMJEpKGIdMsnxluwrvueT5YiMtRjuz2JtW9nXfx1bsW6VtYUq1qtNdZ+XxorCeenHNav/VMiUhlrH8vLZKzkezNW9WG2I0AV8d4EfwrQU0T+JCIRIlJFRFp7xHDCGHNORDpgJeVLVdzH3Zw2A8Pt41g7rPWXm1JY16GOAhki0hfI95aZQjffNsa8hHVgftye2UGss5I59ij/xKpL3YrVEuMbe1hRzcX6BXUS61fJYGNMuv0L8CWss7QjWBf/VhdiuuWxzqhOYv2SPI518RSsi/W/YV2AXoW1w3xQ2MCNMWlYiacv1pnWm1h17d8W8Pt7sVqzxQA77NP4z7DK94wx5iDWRe6/8fu6eIQCrFe7SvFZrCbEp0TkGi+j/h/QBqv10Xzg84LE7uFRETmDVb4fY10D7GSM+a2Q08n2CtavyXh7uuuwfhHnyhgzG+si7TS7KmY7BTtgYYyZiVVGn2JdeJ6D1eorE/gj1kFwH9a6fQ+rBZk3eZbDpcTpYSrWxeGlxphjHsMfxjrYncHa5qfn+N5TwEf2dvAnL9P+FOsM7tMcw2/BOvjsxNqXZpF3dbIb64zmsDFmi8fw9sDXInLWHucvJsc9UiISDfwJeM0+g8r+24fVEixn9VxuXsZKWsewtp0v8xrZ3se+wUqMKz2Gf4tV3j/Y5VYzx/d+xKqSegirmnUz0Mr+eCzwtL09PIH1w/dSFfdxN6d/YJ3hncQ6JuTcDgAwxpwB7sNappNY212+1ZnZF+H8kog8hXUR9WanY1FKhSYR+QD4yRjzuNOxBCv/u7FJKaX8hIjEYDX+uNrhUIJaSPU1p5RSBSUiz2BVkb5oV/8pH/HrqjmllFLBT8+IlFJKOSrgrhFVrVrVxMTEOB2GUkoFlI0bNx4zxlRzOo7cBFwiiomJITEx0ekwlFIqoIjIgfzHcoZWzSmllHKUJiKllFKO0kSklFLKUQF3jSg36enpJCUlce5czs5lVSiIjo6mdu3aREZGOh2KUqoIgiIRJSUlUa5cOWJiYhC51E5sVSAxxnD8+HGSkpKoX7++0+EopYrAZ1VzIvKBiPwiItu9fC4i8qqI7BGRrSLSpqjzOnfuHFWqVNEkFIJEhCpVqujZsFIBzJdnRB9iPXjuYy+f98XqhbcxVu/Jb5FHL8r50SQUuhxf9243xMdDnN3bffZrl8vZuJQKED5LRMaYFXaHgd4MAD62H4y1TkQqikgN+xknSgUGt5svhn7MhrSW8NY2EOG6zO+JmzQJpk7VZKRUATjZaq4WFz76OAkvj9AWkbtEJFFEEo8ePVoiwRXFnDlzEBG+/bZAjxvyavTo0cyaNSvPcf71r39d8L5Tp05FmtdTTz3Ff/7zn1yH16pVi9atW9O6dWsmTJhQpOkXxf79+2nevMDPOHRMZiY8+lRpXGmzeIYn+GfWYzyT+Td6E8/TyQ9hFsU7HaJSASEgmm8bY94xxrQzxrSrVs0ve6gAYOrUqVx33XVMnTrV5/PKmYjWrCnIk68L54EHHmDz5s1s3ryZf/+7IE8fh4yMjGKPwx+dmTafgY228eKmXoyNeJs0IskqVZpzpcpzKx/yJE8zYtldpNzzgFV1p5TyyslEdAio4/G+tj0sIJ09e5ZVq1bx/vvvM23a708CT0hIIDY2liFDhnDVVVcxcuRIsns8f/rpp2nfvj3NmzfnrrvuImdP6EuXLmXgwIHn3y9evJhBgwYxYcIEUlJSaN26NSNHjgTgsssuOz/e888/T4sWLWjVqtX5M5l3332X9u3b06pVK2688UaSk5MLvYzGGB555BGaN29OixYtmD59+vll7NKlCy6Xi6ZNm5KQkEC3bt0YMGAADRo0YMKECUyZMoUOHTrQokUL9u7dC1x85ue5DNn2799Ply5daNOmDW3atPFJwi2sjNlf0HPk5Szc34Q3Iu/njUd/JHLc3TBzJlEzJzNpbCIv9FjIjF3N6fe2i6zhN2kyUioPTjbfdgPjRWQaViOF08Vxfej++2Hz5kudyoVat4aXX857nLlz59KnTx+uvPJKqlSpwsaNG2nbti0AmzZtYseOHdSsWZPOnTuzevVqrrvuOsaPH88TTzwBwKhRo5g3bx433HDD+Wl2796dsWPHcvToUapVq8akSZO4/fbbueGGG3j99dfZnMuCLly4kLlz5/L1119TpkwZTpw4AcDgwYMZM2YMAI8//jjvv/8+9957b57LNHHiRCZPth57//zzz3P27Fk2b97Mli1bOHbsGO3bt6dr164AfPPNN2zfvp369euTkJDAli1b2LVrF5UrV6ZBgwbceeedrF+/nldeeYXXXnuNl/MrUNvll1/O4sWLiY6OZvfu3YwYMcLxvgZffSmd9VntmcpwhqdPh9Pj4PXXz38uLhePjB9PxaWfcxfv8n7KCMbEx+v1IqW88GXz7anAWuAPIpIkIneIyD0ico89ygLgB2AP8C7Wc9wD1tSpUxk+fDgAw4cPv6B6rkOHDtSuXZuwsDBat27N/v37AVi2bBkdO3akRYsWLF26lB07dlwwTRFh1KhRTJ48mVOnTrF27Vr69u2bZxxLlizhtttuo0yZMgBUrlwZgO3bt9OlSxdatGjBlClTLppXbjyr5nr37s2qVasYMWIE4eHhVK9enW7durFhw4bzy+h5H0/79u2pUaMGUVFRNGzYkDi7RVmLFi3OL39BpKenM2bMGFq0aMHQoUPZuXNngb/rCwcPwhMbXfwxbAHDmA5lyvzeWs5TXBx3lv6UbiTwV57nlw5/LPlglQoQvmw1NyKfzw0wrrjnW8Af2sXqxIkTLF26lG3btiEiZGZmIiK8+OKLAERFRZ0fNzw8nIyMDM6dO8fYsWNJTEykTp06PPXUU7neC3Pbbbdxww03EB0dzdChQ4mIKNoqGz16NHPmzKFVq1Z8+OGHJCQkFGk63pQtW/aC957LHBYWdv59WFjY+etIERERZGVlAZCVlUVaWtpF0504cSLVq1dny5YtZGVlER0dXaxxF9Zf/gJZEsFr/4tEtozz3kzb5UKmTeWtaStpNbMrj3zVh49uKfl4lQoEAdFYwd/NmjWLUaNGceDAAfbv38/BgwepX78+K1eu9Pqd7KRTtWpVzp4967WVXM2aNalZsyb//Oc/ue22284Pj4yMJD09/aLxe/XqxaRJk85fA8qumjtz5gw1atQgPT2dKVOmFGk5u3TpwvTp08nMzOTo0aOsWLGCDh06FGlaYD3SY+PGjQC43e5cl+f06dPUqFGDsLAwPvnkEzIzM4s8v0s17x9fM3s2PDFkJzFjelnVcXlVt7lcNPn0Hzz8SBgffwwJg17Ra0VK5UITUTGYOnUqgwYNumDYjTfemGfruYoVKzJmzBiaN29O7969ad++vddxR44cSZ06dWjSpMn5YXfddRctW7Y831ghW58+fXC5XLRr147WrVufb5r9zDPP0LFjRzp37sxVV11VlMVk0KBBtGzZklatWtGjRw9eeOEFrrjiiiJNC2DMmDEsX76cVq1asXbt2ovOqgDGjh3LRx99RKtWrfj2229zHackZMz+gvuevZym7ODBWZ0KlVAev3o+MbKf8XOu14YLSuXGGBNQf23btjU57dy586JhwWTcuHHmvffeczoMv+brbWBm73cNGDMHlzFgzLhxBf/yuHHmE0YaMOZL4gr3XaWKCZBo/OAYntufnhH5ubZt27J161Zuvvlmp0MJaa8eHER92ccfmee9gYI3cXEMLT2f6hzm1bAHCvddpUJAUPS+Hcyyr6Eo52zeDCt3VuGl238mvPSfC9+PnMtF1LSPuOef3/B/G/qxu4nVwaJSyqJnRErl47XXrJOg219qnn8DBW9cLu6e24/ISHjjjeKPUalApolIqTwcOwZTpsCtt0LFipc2rRo14E9/gg8+gDNniiU8pYKCJiKl8vDugztJTYXxTZcWy/Tuu89KQh/1n6Gt55SyaSJSyovM2W7enFyeniym6V9vKJbE0eGwmw5hG3h9ZUvM8BGajJRCE1GxOXz4MMOHD6dhw4a0bduWfv368f333xdpWh9++CE//fRTob9XkEc6NG7cmMGDB1/QVc6dd955/v3MmTNp0qQJ3bt3B2DEiBG0bNmSiRMnFmlZAtnyj/aTZGpzN29DcrL1wLtLFR/P3Vlv8R1XkZjStHimqVSA00RUDIwxDBo0iNjYWPbu3cvGjRt57rnnOHLkSJGml1ciKmrPAtn9xu3evZthw4bRo0cPsp/t9N5779G0aVMA3n//fd59912WLVvG4cOH2bBhA1u3buWBBx4o0HyC6TEQ01JcXMYZ+jO/8E22vYmLY1DpRUSSxrSIUdqUWyk0ERWLZcuWERkZyT333HN+WKtWrejSpQsAL774Iu3bt6dly5Y8+eSTgPV4gyZNmjBmzBiaNWtGXFwcKSkpzJo1i8TEREaOHEnr1q1JSUkhJiaGv/71r7Rp04aZM2de8iMdhg0bRlxcHJ9++ikAsbGxJCYm8vTTT7Nq1SruuOMOHnnkEeLi4jh06BCtW7dm5cqV7N27lz59+tC2bVu6dOly/gGAo0eP5p577qFjx448+uijeY5333330alTJxo0aHBBt0a5PbrC23RKQno6fJYYg6vraUqPu6P4nrbqclFp2lv0jvmOGRXGkPVH7ZFbKcfvqC3sX7H1rDB3rnWH+9y5hf9uDq+88oq5//77c/1s0aJFZsyYMSYrK8tkZmaa/v37m+XLl5t9+/aZ8PBws2nTJmOMMUOHDjWffPKJMcaYbt26mQ0bNpyfRr169czzzz9//v2xY8fOv/773/9uXn31VWOMMU8++aR58cUXL4oht+ETJ04099xzz0Xz83y9b98+06xZs/Pf6dGjh/n++++NMcasW7fOdO/e3RhjzK233mr69+9vMjIy8h1vyJAhJjMz0+zYscM0bNjQGGPMggULzLXXXmt+++03Y4wxx48fz3M6uSnunhUWLrQ6UHC7i3Wy502ebE1/1SrfTF+pnPDjnhVC84ZWtxtGjLDq/SdNKr5fu7mIj48nPj6eq6++GrAeoLd7927q1q1L/fr1ad26NWD1oJDX4xGGDRt2/vX27dt5/PHHOXXqFGfPnqV3796FjsvaLgvu7NmzrFmzhqFDh54flpqaev710KFDCQ8Pz3e8gQMHEhYWRtOmTc9XXeb26Ir8puNr06ZBhQq+qzlzuSA6GqZPh86dfTMPpQJFaCai+HgrCcHvF6EvIRE1a9bMa+/Zxhgee+wx7r777guG79+//6LHQ6SkpHidh2dnn8XxSIdNmzbRrl27Ao+flZVFxYoVc30Yn2d8+Y3nucx5JcP8puNLqakwZw4MGgQe4RarcuWgXz+YORMmToTwcN/MR6lAEJrXiOLirIvPUCwXoXv06EFqairvvPPO+WFbt25l5cqV9O7dmw8++ICzZ88CcOjQIX755Zc8p1euXDnO5HHH46U+0uGzzz4jPj6eESPyfGTUBcqXL0/9+vWZOXMmYCWRLVu2FHk8T7k9uqIo0ykuixbB6dPgcRLqE8OGweHDsGKFb+ejlL8LzUTkclnVcePGFUu1nIgwe/ZslixZQsOGDWnWrBmPPfYYV1xxBXFxcdx0001ce+21tGjRgiFDhuSZZOD3i//ZjRVyKsojHSZOnHi++fbkyZNZunQp1apVK9RyTpkyhffff59WrVrRrFkz5s6de0njZfP26IrCTqe4THspiSrRZ7k++Qufzqd/fygbncH0e1fp/UQqpElhrxU4rV27diYxMfGCYbt27brgWT0q9BTXNpAycx7V/hTLSKbwdpkHfXr9ELebEYNTWZIZy8+lGxIxbbLv5qVCnohsNMYUvD6+BIXmGZFSXiz96CC/cRk38lnx3cTqTXw8QzKncYxqrE1ppTe3qpCliUgpD/My+1KWs3RjefHdxOpNXBy9Sq8mkjTmRQzSm1tVyAqaVnPGGETE6TCUA4qretkYmL8jhl4dfyaq3ZjCP3eosFwuyk+DLmP3MT/8Tp53lffdvJTyY0FxRhQdHc3x48eL7YCkAocxhuPHjxMdHX3J09q+HQ4ehP531ij6c4cKy+Wi/4N/YMeP5TlwwPezU8ofBcUZUe3atUlKSjrfd5oKLdHR0dSuXfuSpzN/vvW/X79LnlSh9O8PDz1kzX/s2JKdt1L+ICgSUWRkJPXr13c6DBXg5s2DNm2gZs2Sne+VV0KjRtb8NRGpUBQUVXNKXarjx2HtWuvspKSJWPNdtuz3Dj+UCiWaiJTC6k0hK8uZRATWfM+dg6XF8yBYpQKKJiKlgPlvJ1Gt9Bna/+xMDwddu1q9LMyfsFJ7WVAhRxORCnkZs79g4Yoy9E35nLCRzjy+O2qRm17pC5i3I0YfIa5CjiYiFfLWT/6ek1S2nsTq694UvImPp3+mmyTqsCOlvvayoEKKJiIV8hZH9UfI4nq+8n1vCt7ExREXvdKKJ7K/9rKgQopPE5GI9BGR70Rkj4hMyOXzuiKyTEQ2ichWESnhOziUgiUHr6Jto9NUGTfCt52c5sXlou70F7my4hGWtHxQOz9VIcVniUhEwoE3gL5AU2CEiDTNMdrjwAxjzNXAcOBNX8WjVG7OnIF166DnkEol15uCNy4XvUZWZ/m31UlLcy4MpUqaL8+IOgB7jDE/GGPSgGnAgBzjGCC7g60KwE8+jEepiyxfDhkZ0KuX05FYevaE336zkqNSocKXiagWcNDjfZI9zNNTwM0ikgQsAO71YTxKXWTxYoiOhk6dnI7E0r07hIVZcSkVKpxurDAC+NAYUxvoB3wiIhfFJCJ3iUiiiCRqf3KqOC1ZYt3DUwx9phaLChWgQwcrLqVChS8T0SGgjsf72vYwT3cAMwCMMWuBaKBqzgkZY94xxrQzxrQr7OOtlfLm0CHYudN/quWy9eoF69fDqVNOR6JUyfBlItoANBaR+iJSCqsxQs679H4ErgcQkSZYiUhPeVSJyD7r6NnT2Thy6tXL6m5o2TKnI1GqZPgsERljMoDxwCJgF1bruB0i8rSIZDdNeggYIyJbgKnAaKMPFVIlZMkSqFYNWrZ0OpILdewIZctq9ZwKHT59DIQxZgFWIwTPYU94vN4JdPZlDErlxsx1s+Tzblzf9gxhYZf+LKPiVKoUxDY5zOIpYdB7nd5TpIKe040VlCp5bjc7hj/D4eQK9Pr6Wf/r183tpteWl9h9+nIODHvU/+JTqphpIlKhJz6epeeuBeD6tAX+169bfDzXpy8EYNm5a/wvPqWKmSYiFXri4lgW3pP6/EC9Msf8r1+3uDiald5HNX5hWXgv/4tPqWKmiUiFnKw/ulheug/dmxx2rm+5vLhcyLSpxDZMYmnFwZgb/Cw+pYqZTxsrKOWPtmyBk2dL0f1vncBfj/EuF90PwcyxsHcvNGrkdEBK+Y6eEamQk31/TvfuzsaRn+z49H4iFew0EamQs3QpNG4MtXL2fOhn/vAHuOIKTUQq+GkiUiElIwNWrPD/syEAESvOZctAb/NWwUwTkQop33xjPYOoRw+nIymY7t3h8GH47junI1HKdzQRqZCSXc0VG+toGAWm14lUKNBEpELKsmXQtClUr+50JAXTsCHUrm1d11IqWGkiUiEj/fMvWPVVKt1j9jkdSoFlXydKiE8la9y92t2PCkqaiFRocLvZcNNEfsuIovuSvwfUAb17xU0c+zWKHW8mwIgRARW7UgWhiUiFhvh4ElKvAaBbWnxA9d/W/dRsABKIheTkgIpdqYLQRKRCQ1wcCWHX04KtVC2TElD9t8UMaUc9OWAlojJlAip2pQpCE5EKCel9XayO7EZsixP+2b9cXlwuYrsLy6N7kzUlwGJXqgA0EamQkJgIyakRxD4VG5AH8thRdTl+7jJ2Ngq82JXKjyYiFRISEqz/Xbs6GkaRZd/3lL0cSgUTTUQqJCQkQIsWULWq05EUTUwM1KuniUgFJ01EKuilp8OqVYHTm4I3sbGwfDlkZTkdiVLFSxORCnqJiVar52BIRMeOwc6dTkeiVPHSRKSCXqBfH8qm14lUsNJEpIJeQgI0bx6414eyxcRA3bqaiFTw0USkglr651+wamkqsfV+cDqUYhEbC8sXa79zKrhoIlLBy+0m8ab/kpwRRezix4PiwB1bwep3bueby7TfORU0NBGp4OXRv1zXtMVB0Udb7Kk5gPY7p4KLJiIVvOz+5ZqzjWplkoOij7aYG9tSV37UfudUUNFEpIJWWh8XqyJj6d7ieOD1L+eFDHBZzyeK7qP9zqmgoYlIBa1A71/Om+x+53Y0DJ5lUqFNE5EKWsuWWf+7dXM2juKWfT9R9vIpFeg0EamglZAALVtClSpOR1K8YmKsP72fSAULnyYiEekjIt+JyB4RmeBlnD+JyE4R2SEin/oyHhU6UlNh9Wro3t3pSHyje3ftd04FD58lIhEJB94A+gJNgREi0jTHOI2Bx4DOxphmwP2+ikeFlg0bICUl8PuX8yY2Fk6cgG3bnI5EqUvnyzOiDsAeY8wPxpg0YBowIMc4Y4A3jDEnAYwxv/gwHhVCli0DkcDvX84bvU6kgokvE1Et4KDH+yR7mKcrgStFZLWIrBORPrlNSETuEpFEEUk8evSoj8JVwSQhAVq1gsqVnY7EN+rWhQYN9DqRCg5ON1aIABoDscAI4F0RqZhzJGPMO8aYdsaYdtWqVSvZCFXASU2FNWuCt1ouW/bziTIznY5EqUvjy0R0CKjj8b62PcxTEuA2xqQbY/YB32MlJqWK7Ov/ruLcOehe5munQ/Gp7t3h1CnYetO/tc85FdB8mYg2AI1FpL6IlAKGAzn3ljlYZ0OISFWsqrrg6CZZOcPtZtmTCQhZdJk4OKgP0LFpVj9zS2cc1Q5QVUDzWSIyxmQA44FFwC5ghjFmh4g8LSLZt4QvAo6LyE5gGfCIMea4r2JSISA+nqXpXWjDN1RK+SmoOwWt/Y2bxnzPMrprB6gqoPn0GpExZoEx5kpjTENjzLP2sCeMMW77tTHGPGiMaWqMaWGMmebLeFTw+61rX9ZyLdfzVfB3ChoXx/URy1lON9JLlw/uZVVBzenGCkoVq9UV+5NOKXrccFnQdHTqlctFjwdbc5ZyJD41L7iXVQU1TUQqqHz1FURGwnVTx4XEgbn7I+0B+Cqti8ORKFV0mohUUFm6FK65BsqWdTqSklG1qnW/1NKlTkeiVNFpIlJB4+RJ2LgRrr/e6UhK1vXXW/dNpaQ4HYlSRaOJSAWN5cvBGOjRw+lISlaPHr/fxKtUINJEpILGV3ZDuY4dnY6kZHXtCuHhWj2nApcmIhU0li6FLl2gVCmnIylZ5cpBhw5WIlYqEGkiUkHh5w8XsXMnXH/FDqdDccT118OGDYbTYx7WHhZUwNFEpAKf283Su6cD0GP6XSF5IO4RtZqsLGHFe99pdz8q4GgiUoEvPp6v0q6jEidofW5dSHZ1c23STKJJYQk9tbsfFXA0EamAZ3rFES+96ckSwstEh2RXN9H9etAtbBWL6RX8XRupoKOJSAW8XY1dHDK1iOueHvzd+njjchE3uga7aErSa7NDswxUwNJEpALe4sXW/16TRob0AbjX/c0BWCx6NqQCS76JSETuFZFKJRGMUkURHw9XXgn16jkdibOaN4crrtDLQyrwFOSMqDqwQURmiEgfERFfB6VUQaWmQkKCXhIBELHKYfFiyMpyOhqlCi7fRGSMeRzr8d3vA6OB3SLyLxFp6OPYlMrXmjVWIzFNRJa4ODh+HDZtcjoSpQquQNeIjDEGOGz/ZQCVgFki8oIPY1MqX/HxEBEBsbFOR+Ifeva0/mv1nAokBblG9BcR2Qi8AKwGWhhj/gy0BW70cXxK5Sk+Hjp1srq5UVC9OrRurYlIBZaCnBFVBgYbY3obY2YaY9IBjDFZwB99Gp1SeTh6FL75RqvlcoqLg9Wr4exZpyNRqmAKco3oSWPMAS+f7Sr+kJQqmCXPbwQgLnqFw5H4l7g4SE+H5cPf0q5+VEDQ+4hUYHK7iX95J5U4QZt/9NcDrofOJ+dRmmQWzU/XfudUQNBEpAJS1qLFLMiMozeLCE85qxdFPEQnfEl3lrGAfhjtd04FAE1EKiB902AIv1Cd/szXvtVyioujf+Ri9tKI3dEttWyU39NEpALS/LPdEDH0ub1W6PYv543LRb83+gMw/6YpWjbK74l1i1DgaNeunUlMTHQ6DOWwDh2sx2OvXet0JP6rWTOoUQOWLHE6EuUPRGSjMaad03HkRs+IVMA5cgQ2bID+/Z2OxL/17w8rVsCZM05HolTeNBGpgPPll9Z/TUR569/fasatZ0TK32kiUgFn/nyryql1a6cj8W+dOkGFClZ5KeXPNBGpgJKebrVG7tfP6m1aeRcZaTWYW7AAAuxSsAoxmohUQFnz4mpOn4b+1dY7HUpA6N8ffv4ZNg//t97YqvyWJiIVONxu5j+5nkjS6PnKDXpgLYA+sgiA+TPOai8Lym/5NBHZD9L7TkT2iMiEPMa7UUSMiPhl00LlJ+LjcWf0pRvLKZfyi/YYUADV139BB77Gjct6cJOWmfJDPktEIhIOvAH0BZoCI0SkaS7jlQP+Anztq1hUcNjVZDDfcRUDmaO9KRRUXBwDIxewgQ4kRTfSMlN+yZdnRB2APcaYH4wxacA0YEAu4z0DPA+c82EsKgjM+bUHAANvrai9KRSUy8Wgl7sBMPeWz7TMlF/yZSKqBRz0eJ9kDztPRNoAdYwxeTYwFZG7RCRRRBKPHj1a/JGqgDB7ttWjQq0Pn9UDaiFcNbYHV10Fs/e2dDoUpXLlWGMFEQkD/gs8lN+4xph3jDHtjDHtqlWr5vvglN9JSrJ6Uxg0yOlIAtOgQZCQACdOOB2JUhfzZSI6BNTxeF/bHpatHNAcSBCR/cA1gFsbLKjczJlj/R840MkoAtfAgZCZCfPmOR2JUhfzZSLaADQWkfoiUgoYDpxvO2qMOW2MqWqMiTHGxADrAJcxRns0VReZPRuuusr6U4XXrh3UqmWVo1L+xmeJyBiTAYwHFgG7gBnGmB0i8rSIaAW/KrDjx2H5cq2WuxRhYdZZ0aJFVitupfyJT68RGWMWGGOuNMY0NMY8aw97whhz0V11xphYPRtSuZn3zCYyM2FQlRVOhxLQBg2ClBRYNOQdvbFV+RXtWUH5N7eb2a8nUZuDtPtHXz2AXoKup7+gEieYvTBae1lQfkUTkfJrp79YwZeZvRjM50iK9gxwKSKXLmIgc5jLAM4lZ2pZKr+hiUj5tTllR5JKNCOYqr0pXKq4OIZHzeFXKrCw1EAtS+U3NBEpvzb126uJufw3Oo5tp70pXCqXix5Tx1Ct9Bmmtn9Jy1L5jQinA1DKm6NHraeLPvJIWeS5150OJyhEDLqBobfBBx+U48wZKFfO6YiU0jMi5cdmzbJuwhwxwulIgsuIEXDunLZVUP5DE5HyW1OnQpMm0KKF05EEl06doHZtq3yV8geaiJRfSkqClSutX+/6SPDiFRYGw4dbN7dq33PKH2giUn5p+hM7ABhe7SuHIwlOI0ZARgZ8duOnWkenHKeJSPkft5tPP0yjLYk0fsilB0ofuPqgm8aym08TaujNrcpxmoiU39n66Xa+MVczik/08dY+IovjGWU+JoHu7Eu+XMtYOUoTkfI7k5KHUYpUbmay3sTqK3FxjI6ejpDFhxFjtIyVo8QY43QMhdKuXTuTmKh9owartDSoWRN6XHWIGa2fsw6QeuOlb7jd9Lm3ETuTY9h3uAzh4U4HpHxJRDYaY/zyeW96Q6vyK198YT324fbHa0EfvYnVp1wubj8Hw4bBV1/pSZFyjlbNKb/ywQfWA9x69XI6ktAwYABUrmyVu1JO0USk/MahQ/DllzB6NFpNVEKiomDkSOvJrXpPkXKKJiLlNz5+bCdZWXBbnSVOhxJS7rjDujb36cAZ2oxbOUITkfILWXPcfDC5FLEso+GDA/SAWIJaHXDTRjbx3sorMcP1niJV8jQRKb+w6P0k9phG3MU7eu9QSYuP527zFltozeqUq7XsVYnTRKT8wms/30gN+ZkhzNJ7h0paXBwjS8+mIid5Lfx+LXtV4jQRKcd9/z0s3FidP484TeS4u/UBeCXN5aLstPe5s/VGPjODSWqjZa9KliYi5bg33oBSpeCu/14Fr7+uScgJLhdjP+9Jlgnjf/9zOhgVajQRKUf9+itMmmTdVFm9utPRhLb69a3fAG+/bT04T6mSoolIOeqjh7dx5gzc22qF06Eo4N574dgxmNb/E209p0qMJiLlmKw5bl5/L4prWEv7J/rqgc8P9Djrpqns5NWlzbQptyoxmoiUY2a//Qvfmyu5n5e1ybafkMXxPGD+yybasDils64TVSI0ESlHGAPPfj+EK+V7bbLtT+LiGFX6M2qRxLNh/9B1okqEJiLliIULYdMPFZkw/jfCx/1Zm2z7C5eLqGkf8ch161iR1YWVlXSdKN/T5xGpEmcMdO5sdXK6Zw9ERjodkcopORliYqBNG6sjWhX4/Pl5RHpGpEpcwrOrWbsW/tp3qyYhP1WmDDz4ICxaBIl/ekEbLSif8mkiEpE+IvKdiOwRkQm5fP6giOwUka0i8pWI1PNlPMoPuN08+2QaV/Azt38cqwc4PzY2ZgEVOcm/ZjaCEdqCTvmOzxKRiIQDbwB9gabACBFpmmO0TUA7Y0xLYBbwgq/iUf4hYdI+vsrqzsP8h+iUk9oqy4+VX7WA+3iV2QxmY/JVuq6Uz/jyjKgDsMcY84MxJg2YBgzwHMEYs8wYk2y/XQfU9mE8ymHGwKM7bqG2JDGWN7WlnL+Li+PB0v+jKkd5NOwlTC9dV8o3Inw47VrAQY/3SUDHPMa/A1jow3iUw2bOhA27KzHpvv2UzrzDSkLaUs5/uVxUmAb/eDGBv6wayqIo6ON0TCoo+azVnIgMAfoYY+60348COhpjxucy7s3AeKCbMSY1l8/vAu4CqFu3btsDBw74JGblO2lp0KQJlC0Lmzbpo8ADia674BCqreYOAXU83te2h11ARHoCfwdcuSUhAGPMO8aYdsaYdtWqVfNJsMq33h6/jR9+gBcGr9UDWYApVQqeew62bYPJvbUPOlX8fJmINgCNRaS+iJQChgMXbMEicjXwNlYS+sWHsSgHnfp0AU+/ewU9+IreL/TUA1kAGhrlpn1YIo9/FUvy8Nt1Hapi5bNEZIzJwKpuWwTsAmYYY3aIyNMikn1h4EXgMmCmiGwWEd26g9Dfn7uME1TmRR5BUrRPuUAki+N5KesBkqjDMykP6TpUxcqXjRUwxiwAFuQY9oTH656+nL9y3vr18NaOLtwb8RZtMjZpS7lAFRdHl0mTGJ08if/wMCP/sIrmTsekgob2rKB8JiMD7r4batQQnvkkBsaN0z7lApXLBVOn8uLt31KhXBb3TO9OVpbTQalgoYlI+cyrd29n82Z4ddQGyg/vp48BD3QuF1Xff54XX4li9Wr4oOcUvVakioUmIuUT+99dzBMfxNCfeQx+NVYPWEFkdGU3XcNW8uiyvhwefr+uW3XJNBGpYpeZCaOerE8YWbzOeG2gEGRkcTxvZ40hhdLclvIGZpGuW3VpNBGpYvfvf8OqnxvxRqkHieGANlAINnFxXFXmIP/hYb6kL6+n3+10RCrA+bTVnAoxbjfrP/mOJz9/iOHDw7h5uAsWR2tXPsHGbrgwdlE8C74+wiMfNqX7mWdpPqKFrmdVJPpgPFU83G7ODr+Tq1NWkyZRbJm8nYo39XM6KuVjRz5eRMtbr6Y6h1lfOpboaR9qMvJTodrFjwohWYsWc2vKW/xAAz42o6i4ZkH+X1IBr/r6L5jEaLbRkntS/qvXi1SRaCJSxeKZk+P4nBt5kUfoViZRrwmFirg4+pVZzpM8xUeM5uVz9zgdkQpAeo1IXRq3m9n/O8xTC+/ilu4HeaBJGvTWm1ZDhn296IlFi9n6zc88/GFzmtVdS9zRKXptUBWYXiNSRed2s2XYv+h8bgnNwnaxfPoRoof80emolEPOnoVrm/9K0oFM1nItV5U5qD1p+BG9RqSC0vczNhN3bi6VOMnsrAFEJ3zpdEjKQZddBu6uL1GKNHqxmAPJVfX+MVUgmohU4bnd/HjrP+j5xX0YwlhCT2qWOa3XhRT1h7RlUfRAzlCOnvIVR8Jrwvjx2vuCypNWzanCcbs5MvwvdE35kiNUJ+GOybSO/lavB6jfud2s+Wg3vb64l0YZ37HMdKNymVStpnOYVs2poLF/ViLXpcSTRG3m099KQtqZqfLkctHps4eY2/dtvjVX0o3l/JRcQavplFeaiFTBuN1sH/EsnWc/zHGqsISedC6zWavjlFc976jHwqhB7CeGzrKGPZn1tZpO5Uqr5lT+3G5W/elVXKkziJZU4u+YQfOo3Vodp/LndpM4+Vv6zhtHWMpvzKcf7crs0mo6B2jVnApYxsCbL6XQPXUhVTnGatPJSkJaHacKwuWi3YxHWTnov5QmmetYxUfJQ7SaTl1AE5Hy6tysedzRdC3jVgyjd9gS1tOB+mV+0eo4VWhXDWvFhtLd6MxqRvMR964bSfqf79NqOgVozwrKi80vJ3DLgzFsM815IuJfPPlIMmG/3qzVcapoXC6qTYNFX85lwjeleOnrLqzbGMHHH95Nk+noNhXi9IxIXSBj9hc8e80XdHiwM0dNFebTj//L+Dthv57S6jh1aVwuIt58lf+0m85nDGYf9bn63BomTjhC1rh79ewohGkiUueteX4lHW+szeNf38Bgmc32Um3px0J9sJ0qXnFxDC6ziO00J06W8OCuMXR+8yY2Dntek1GI0kSk+GnSIkb9YT2dJ3ThiKnGDIYyLWsYVXpeDePGaQsnVbzsjlKvGDeEuX3e4mNG8QMNaH9uJWPGl+KX2ydoQgox2nw7VLndHJmzlud39Oet9W3JIoyHwybyt4gXKJt20joL0gSkfM3thhEjOJ0cwf+FPc1rWWOJIpV7I/7Hw+PPUSX9sF6XLCb+3HxbE1EocbshPp59JobX3onifxl3kEoUo/iEf/AMDfkB+vWD+vV151clx94u2bePbxfs5f94kukMoyy/MY43GBf9AXUeHAqnT+t2eQn8ORFhjAmov7Zt2xpVeJmz55olUf3MAGYbIdNEkGZG8ZH5nkbGREQYA8aUKWPM3LlOh6pC1dy51jYIZntYCzOMqSaMDBNOuhkqM8xyupis0rqNFhWQaPzgGJ7bn54RBbnv3vyKTz5I55Mdbfjx3OVU5Sh38zZ/Dn+XWpk/WlVw99+vvzaVf8g+O6pQAV5+mf3J1XhTxvOuuYNTVKIBe7ml8Vpu7riHhkPb6PZaCP58RqSJKFjYO3BWzzg2/VCBOVOT+Xx/G3Yeq04YmcTJEm4Nn8yAjFmULhOmyUf5P4+k9NvEd/g8pQ8fyW0sNbEYwmgp2xjU+RcG1lxPq5uaIYI1vm7TudJEVIw0EXlwuzGL4tmd2YBVk3azOK0rS+jJMaoRRiZdWckgPmcIs6jJz3r9RwUuj+tIPy7YxkyGMptBrKEThjCqc5heYUvpmbWI66I30uCBAciv+kPLkyaiYhRyiSh7B7Tv4zkyZy3f1BnAxr0VSJy2hzWZHTnK5QBcwc/0YjG9WEwfvqQaxyAiAjIytBWcCg52KzuSkyE8nMOZVVlIXxbTiyX0PL8vVOcwnVhDu8gttL0xhrZlvqXqgM7W9u+5T4XQ/qCJqBgFVSLKkWSyX6elCwdmf8OejBh2f76N79Nj2CHN2WGant/RABrzPZ1YQ2dW0zn8a5pkbkNKlQIRSE3V6z8qOOW4jkRyMpQqRRZh7EhrxGrpwmpzLWu5lr00Ov+16hyhWb2zNDu0iCszdtI46iCNbu9KXXOAyL49rZFy2R+DZb8J2UQkIn2AV4Bw4D1jzL9zfB4FfAy0BY4Dw4wx+/Oapl8kIi8JJL/Xad17c+JsKY5+uZEjcgVHPlvF4bRKHAqryyFqkpRVkwPU4ydqYjzuNS7PaZqyk+Zspxk7aMM3tO51OeVXL7R2wpwJx3PeQbITKZWr3PZFjwR1Kqwym7Jasomr2UEztpduz86U+pyl3PlJhJFJTX4iRn6ktvmRmmFHqCU/cUVmEtWjTlN9VBzVMg9T+Y+diIwwRdr3/WFfDMlEJCLhwPdALyAJ2ACMMMbs9BhnLNDSGHOPiAwHBhljhuU13SInooIkD5cLM9dN2pdLSevWi9SMcFKXrCS1U3fOpYVxbuV6UqIqkTLjC5JTw0gOL0+ylOW3jFKcDa/AWcrxa2YZfg2ryK+U53RWOU5KZU6aSpygEmcon2toZfiNWhyiFoeoxwFi2E8M+2nEHhqH7+PyzJ8uPtOZOvWi2JVSttzOmuwfbGbiyxxOKc+esD+wO6sB+4nhAPU4QD2SqM0hanGO0rlOtjynqcwJKskpKnKaiuYE5cPOUl7OUC7zFJeFp3CZ/MZlGacoE55GGUmhdMavlI4ylB45mNIpJ4iKvZboyEyiVi+lVM+uREVmEbl0EdI77uKqQyi2pBaqieha4CljTG/7/WMAxpjnPMZZZI+zVkQigMNANZNHUEVKRG43HwxdyH/S7iWDSDKIIINw0okk3X6fRinSw6LIyAov/MLawsmgPL+e/6vAaSpxkoqcohInqcoxqnCcahylOkeoHn6c6pmHqBCZgoTZSSavqjXQxKNUYeW8JuSlai97vzORpTgllTicVokjYTU5klWVY1TlOFU4RlVOUolTVOQklThNBY89vjxZFP34EUE6pSINkenJRJJOxPmjVToRZBIhGTxpnmJYmXlFut7rz4nIl4+BqAUc9HifBHT0No4xJkNETgNVgGOeI4nIXcBdAHXr1i18JPHxVE07RAu2EU6mnYYy7TRkrfBSpFGqSgUij/5EKdKIIpVI0onmHFGkEkUqpUmx/sLTKJN5hjKRGZSWc1yWdpzLItOIkjQkLUcy8fY6ryTj+TrnxqYJSKnCcbku3G8833fseNF+J3FxVAIqxcfTxEuyyu21SU0lNbIcZ6UcZ9MiSY6sQDJlSU6P4FxYWVKySpFCaftoEsU5okmzjjykEmUdjSrUIO3Yr3b6sY5QmYSTQQSZJpzKnLBiiY8PqmOBL8+IhgB9jDF32u9HAR2NMeM9xtluj5Nkv99rj3Mst2lC0c+Izre0yS8x5LfRFTSBFOR1EG1ISgWtIl4T9nbdqliOQUF2RhQaVXNQ8HpXH9XPKqVCnMPHoFBNRBFYjRWuBw5hNVa4yRizw2OccUALj8YKg40xf8prun7Rak4ppQKMPycin10jsq/5jAcWYTXf/sAYs0NEnsbqfM8NvA98IiJ7gBPAcF/Fo5RSyj/5srECxpgFwIIcw57weH0OGOrLGJRSSvk3fUKrUkopR2kiUkop5ShNREoppRyliUgppZSjAq73bRE5Chwo4terkqPXBj+hcRWOxlV4/hqbxlU4lxJXPWNMteIMprgEXCK6FCKS6I/t6DWuwtG4Cs9fY9O4Csdf47pUWjWnlFLKUZqIlFJKOSrUEtE7TgfghcZVOBpX4flrbBpX4fhrXJckpK4RKaWU8j+hdkaklFLKz2giUkop5aigTkQi8qKIfCsiW0VktohU9DJeHxH5TkT2iMiEEohrqIjsEJEsEfHaFFNE9ovINhHZLCI+f/ZFIeIq6fKqLCKLRWS3/b+Sl/Ey7bLaLCJuH8aT5/KLSJSITLc//1pEYnwVSyHjGi0iRz3K6M4SiusDEfnFfhBmbp+LiLxqx71VRNr4SVyxInLao7yeyG28Yo6pjogsE5Gd9r74l1zGcaS8fMoYE7R/QBwQYb9+Hng+l3HCgb1AA6AUsAVo6uO4mgB/ABKAdnmMtx+oWoLllW9cDpXXC8AE+/WE3Naj/dnZEiijfJcfGAv8z349HJjuJ3GNBl4vqe3JY75dgTbAdi+f9wMWAgJcA3ztJ3HFAvNKuKxqAG3s1+WwnumWcz06Ul6+/AvqMyJjTLwxJsN+uw6onctoHYA9xpgfjDFpwDRggI/j2mWM+c6X8yiKAsZV4uVlT/8j+/VHwEAfzy8vBVl+z3hnAdeLiPhBXI4wxqzAet6YNwOAj41lHVBRRGr4QVwlzhjzszHmG/v1GWAXUCvHaI6Uly8FdSLK4XasXxE51QIOerxP4uIV7xQDxIvIRhG5y+lgbE6UV3VjzM/268NAdS/jRYtIooisE5GBPoqlIMt/fhz7h9BpoIqP4ilMXAA32tU5s0Skjo9jKih/3gevFZEtIrJQRJqV5IztKt2rga9zfOTP5VUkPn0wXkkQkSXAFbl89HdjzFx7nL8DGcAUf4qrAK4zxhwSkcuBxSLyrf0rzum4il1ecXm+McYYEfF2z0E9u7waAEtFZJsxZm9xxxrAvgCmGmNSReRurLO2Hg7H5M++wdqmzopIP2AO0LgkZiwilwGfAfcbY34tiXk6KeATkTGmZ16fi8ho4I/A9cauYM3hEOD5y7C2PcyncRVwGofs/7+IyGys6pdLSkTFEFeJl5eIHBGRGsaYn+0qiF+8TCO7vH4QkQSsX5PFnYgKsvzZ4ySJSARQAThezHEUOi5jjGcM72Fde/MHPtmmLpVnAjDGLBCRN0WkqjHGp52hikgkVhKaYoz5PJdR/LK8LkVQV82JSB/gUcBljEn2MtoGoLGI1BeRUlgXl33W4qqgRKSsiJTLfo3V8CLX1j0lzInycgO32q9vBS46cxORSiISZb+uCnQGdvogloIsv2e8Q4ClXn4ElWhcOa4juLCuP/gDN3CL3RrsGuC0R1WsY0TkiuxreyLSAet46dMfFPb83gd2GWP+62U0vyyvS+J0awlf/gF7sOpSN9t/2S2ZagILPMbrh9U6ZS9WFZWv4xqEVa+bChwBFuWMC6v10xb7b4e/xOVQeVUBvgJ2A0uAyvbwdsB79utOwDa7vLYBd/gwnouWH3ga6wcPQDQw097+1gMNfF1GBYzrOXtb2gIsA64qobimAj8D6fb2dQdwD3CP/bkAb9hxbyOPlqQlHNd4j/JaB3QqgZiuw7o2vNXjuNXPH8rLl3/axY9SSilHBXXVnFJKKf+niUgppZSjNBEppZRylCYipZRSjtJEpJRSylGaiJRSSjlKE5FSSilHaSJS6hKJSHu7I9Fou0eMHSLS3Om4lAoUekOrUsVARP6J1aNCaSDJGPOcwyEpFTA0ESlVDOz+3TYA57C6gsl0OCSlAoZWzSlVPKoAl2E9VTPa4ViUCih6RqRUMRARN9ZTUesDNYwx4x0OSamAEfDPI1LKaSJyC5BujPlURMKBNSLSwxiz1OnYlAoEekaklFLKUXqNSCmllKM0ESmllHKUJiKllFKO0kSklFLKUZqIlFJKOUoTkVJKKUdpIlJKKeWo/weQS5FG9spLSgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def analfunction4(x):\n",
    "    y = 1 - tanh(2*x)**2\n",
    "    return y\n",
    "y_theo = []\n",
    "for i in range(0,len(x_set)):\n",
    "    y_theo.append(analfunction4(x_set[i]))\n",
    "\n",
    "#create the graph\n",
    "fig3 = plt.figure()\n",
    "ax3 = fig3.add_subplot(111)\n",
    "\n",
    "ax3.set_xlabel(\"x\")\n",
    "ax3.set_ylabel(\"y\")\n",
    "ax3.scatter(x_set,y_cd, s = 10, c = \"r\", label = \"Central Difference\")\n",
    "ax3.plot(x_set,y_theo, c = 'b', label = \"Analytical Formula\")\n",
    "    \n",
    "plt.legend(loc = \"upper left\")\n",
    "plt.title(\"Comparison of Central Difference Derivative vs Analytical Formula\")\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyPSuJbs/vagEh3+73Nx1fhV",
   "name": "lab_assignment_module2.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}