{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "#**Lab 3**\n",
        "\n",
        "**Name: (your name)**\n",
        "\n"
      ],
      "metadata": {
        "id": "6VpZ01tyznHo"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Some remarks before you start:\n",
        "* Import the Python packages below each time you start a session.\n",
        "* Before submitting your report, remove all ouputs by clicking on **Edit🠦Clear all ouputs**. Then download your report by clicking on **File🠦Download🠦.ipynb**. This is the file you will submit on Canvas.\n",
        "\n",
        "In this lab, you will practice:\n",
        "* define and graph functions\n",
        "* use Bisection method and Newton-Raphson method to find roots\n",
        "* use Gradient Descent/Ascent method to find local min/max\n",
        "* solve problems of projectile motion and mortgage rate\n",
        "\n",
        "\\begin{array}{|c| c|}\n",
        " \\hline\n",
        " \\text{Problems} & \\text{Points}\\\\\n",
        " \\hline\\hline\n",
        " 1, 3, 5, 7 & 3 \\\\\n",
        " \\hline\n",
        " 2, 4, 6, 8, 9 & 4 \\\\\n",
        " \\hline\n",
        " \\text{Readability of your report} & 3 \\\\\n",
        " \\hline\n",
        " \\text{Total: 9} & \\text{Total: 35} \\\\\n",
        " \\hline\n",
        "\\end{array}\n",
        "\n"
      ],
      "metadata": {
        "id": "3cVoL2qvzufn"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##**I. Import necessary Python packages**\n",
        "Execute the following code each time you start a session."
      ],
      "metadata": {
        "id": "bQagzUIHzzb_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from matplotlib.pyplot import*\n",
        "from autograd.numpy import*\n",
        "from autograd import elementwise_grad as diff"
      ],
      "metadata": {
        "id": "sGBvQoGSzoih"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "## **II. Define and graph functions**\n",
        "\n",
        "To define a function on Python, you will use the keyword **def**. For example, the code below defines the function $$f(x)=\\cos^2(x)-\\cos(x^2)$$"
      ],
      "metadata": {
        "id": "yGFzT_PXz431"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def f(x):\n",
        "  return cos(x)**2-cos(x**2)"
      ],
      "metadata": {
        "id": "kxIif3jP1xeP"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "To graph the function, you first need to generate an array of $x$-values using the command **linspace**. For example,"
      ],
      "metadata": {
        "id": "QqnmIFJyDe1-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "x=linspace(0,6,16)  # Generates 16 evenly-spaced points between 0 and 6\n",
        "print(x)"
      ],
      "metadata": {
        "id": "_fLg6iO5z5sG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "*Note that 16 is the number of points, not the number of subintervals. The number of subintervals in this case is 15.*\n",
        "\n",
        "At these values of $x$, the command **plot** will connect the points $(x,f(x))$ together by straight lines."
      ],
      "metadata": {
        "id": "PZSxjPbZFM2L"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plot(x,f(x))"
      ],
      "metadata": {
        "id": "ndThAdz3nzYV"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "You can label the axes and give a title to the plot as follows."
      ],
      "metadata": {
        "id": "los8m4Es9GA6"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plot(x, f(x))\n",
        "xlabel(\"x\")\n",
        "ylabel(\"f(x)\")\n",
        "title(\"Plot of $f(x) = \\cos(x)^2 - \\cos(x^2)$\")"
      ],
      "metadata": {
        "id": "3Bti5ooH8P4i"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "You can graph two or more functions on the same graph. For example\n",
        "\n",
        "$$f(x)=x^2,\\ \\ g(x)=1-x$$"
      ],
      "metadata": {
        "id": "bZzavLFJEyhG"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "plot(x, x**2, label=\"graph of $f(x)=x^2$\")\n",
        "plot(x, 1-x, label=\"graph of $g(x)=1-x$\")\n",
        "xlabel(\"x\")\n",
        "legend()\n",
        "title(\"two functions\")"
      ],
      "metadata": {
        "id": "rp9vTVKqBVC_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###<a name=\"graph\"></a>**Exercise 1**\n",
        "Graph the function $g(x)=\\cos(x)-\\cos(\\sqrt{2}x)+\\frac{1}{3}\\sin(3x)$ on the interval $[0,12]$. In your **plot** command, use 100 equal *subintervals* (therefore, the number of equally-spaced points is 101). Give your plot a title."
      ],
      "metadata": {
        "id": "DIQOwXFQHl-b"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 2**\n",
        "To find the derivative of f, use the command **diff**. For example,"
      ],
      "metadata": {
        "id": "Go2IWVSOKeW_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "df = diff(f)\n",
        "x=linspace(0,6,16)\n",
        "plot(x,df(x))"
      ],
      "metadata": {
        "id": "TVtBstKMCF-c"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "To evaluate the derivative of $f$ at $x=1$, namely $f'(1)$, just run the command **df(1.0)** or **df(1.)**. If you run **df(1)**, Python will give an error message because it wants a float-point number for its argument instead of an integer."
      ],
      "metadata": {
        "id": "jz1WNU_pCyuf"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "* With the function $g$ in [Exercise 1](#graph), evaluate $g'(0)$, $g'(\\pi)$, and $g''(1)$. *Hint: $g''$ is the derivative of $g'$.*\n",
        "* Graph $g$, $g'$, $g''$ together on the same plot on the interval $[0,12]$. In your **plot** command, use 100 equal subintervals."
      ],
      "metadata": {
        "id": "9OfTiqaLF-bl"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##**III. Bisection method**\n",
        "Suppose a continuous function $f$ has different signs at $x=a$ and at $x=b$. *Intermediate Value Theorem* guarantees that $f$ must have a root between $a$ and $b$. You can find approximately this root using the Bisection method, i.e. by the process of halfing the search interval until a certain number of divisions is reached or until a desirable precision is reached. If the original search interval $[a,b]$ is divided $n$ times, the Python code is as follows."
      ],
      "metadata": {
        "id": "Z_m2zyE8PBDJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def bisection(f,a,b,n):\n",
        "  for i in range(n):\n",
        "    c = (a+b)/2\n",
        "    if f(c)*f(a)<0: b = c\n",
        "    else: a = c\n",
        "    print(i, c)\n",
        "  return (a+b)/2"
      ],
      "metadata": {
        "id": "0x6z6B0_XHkG"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 3**\n",
        "* Explain what each line of the code above does.\n",
        "* Use this code, with $a$ and $b$ of your choice, to find approximately $\\sqrt{2}$ (root of $f(x)=x^2-2$). How large must $n$ be to get 4 correct decimal places?"
      ],
      "metadata": {
        "id": "F_ezhh0EZjX_"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 4**\n",
        "* Graph two curves $y=xe^{-x}$ and $y=\\frac{1}{x^2+1}$ on the same plot.\n",
        "* From the graph, how many intersection points do you see? Give rough estimates for the coordinates of those intersection points.\n",
        "* Use the Bisection method to find the coordinates of those points. The $x$-coordinates must be correct with an allowable error of $0.001$."
      ],
      "metadata": {
        "id": "J9xAgbPfOvHc"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##**IV. Newton-Raphson method**\n",
        "In the Newton-Raphson method to solve $f(x)=0$, you start with an initial guess $x_0$. Then successively update this guess by $x_1,x_2,x_3,...$ using the recursive formula\n",
        "\n",
        "$$x_{n+1}=x_n-\\frac{f(x_n)}{f'(x_n)}$$\n",
        "\n",
        "You stop when a certain number of iterations is reached or when a certain precision is reached, i.e. $|x_{n+1}-x_n|$ is less than some prescribed $\\epsilon$. The Python code for doing $n$ iterations is as follows."
      ],
      "metadata": {
        "id": "yTs3lEbvc53r"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def newton(f,x0,n):\n",
        "  df = diff(f)\n",
        "  x = float(x0)\n",
        "  for i in range(n):\n",
        "    x = x - f(x)/df(x)\n",
        "    print(i, x)\n",
        "  return x"
      ],
      "metadata": {
        "id": "ZXYoR04oWWHT"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 5**\n",
        "* Explain what each line of the code does.\n",
        "* Find $\\sqrt{2}$ by finding the root of $f(x)=x^2-2$. Use initial guess $x_0=2$. How large is $n$ in order to obtain a precision of $0.00001$ ? In other words, how large is $n$ such that $|x_{n+1}-x_n|<0.00001$ ?"
      ],
      "metadata": {
        "id": "PoZ9vnh11OPp"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 6**\n",
        "In projectile motion, the horizontal $x$ and the vertical $y$ positions are functions of time. From Newton's Second Law of Motion:\n",
        "\n",
        "$$x=v_0\\cos(\\theta)t,\\ \\ \\ y=h+v_0\\sin(\\theta) t-\\frac{1}{2}gt^2$$\n",
        "\n",
        "where:\n",
        "\n",
        "$x$ = horizontal displacement\n",
        "\n",
        "$y$ = vertical displacement\n",
        "\n",
        "$h$ = initial height\n",
        "\n",
        "$v_0$  = initial velocity\n",
        "\n",
        "$\\theta$ = angle of projection\n",
        "\n",
        "$g$ = $9.81$ m/s$^2$ = acceleration due to gravity\n",
        "\n",
        "From the first equation, one can extract $t$ in terms of $x$ and then substitute $t$ into the second equation. One arrives at the formula:\n",
        "\n",
        "$$y=h+x\\tan(\\theta)-\\frac{gx^2}{2v_0^2\\cos^2(\\theta)}$$\n",
        "\n",
        "The *horizontal range* is the value of $x$ such that $y=0$."
      ],
      "metadata": {
        "id": "tUzvkc6u3MWQ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "![shooting.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAggAAADmCAYAAACwLX24AAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAAFiUAABYlAUlSJPAAADb+SURBVHhe7Z0HmFNl2obtq7u2Xbuuil3RtWFBFJGiSBMQpEmRDkpHacICAgKCCtJEAelioalUBQGpUgQGqYJ0BQQEpSr7/f/9bg4bZoZhMpNkTpLnvq73ysxJ5mSSyZzv+d56mhNCCCGESIYEghBCCCFSIIEghBBCiBRIIAghhBAiBRIIQgghhEiBBIIQQgghUiCBIIQQQogUSCAIIYQQIgUSCEIIIYRIgQSCEEIIIVIggSCEEEKIFEggCOFz/vjjD7dv3z63ZcsWt3TpUjdjxgw3c+ZMX9iSJUvs9+L34/cUQsQPEghC+ISDBw+6rVu3uqSkJDd9+nQ3ZMgQ9/rrr7sWLVq4evXquXLlyrlixYq5MmXKmD333HNZZt7zlypVyj3//PPupZdest+zS5cubtiwYe7LL7808bBx40b322+/BV6hECKWkEAQIoocPnzY7dixw61cudJEwPjx421B7d69u3vllVdsoa1du7YtukWLFnXly5d3DRo0cG3btnVvvPGG++CDD9xXX31lC/DUqVOzzHh+bOTIkfa7//vf/7bfs0qVKiZk+P1r1Kjh6tat65o3b378dx8zZoybPHmyW7RokXke9u/f744ePer+85//BN4hIYRfkEAQIoLgdv/111/dsmXL3Pvvv2+LZc2aNW0HziKKGGjYsKFr2rSpeQsQCxMnTnSzZs2yHfj27dvdgQMH3LFjxwJn9Ccs8L///rv7+eef3YoVK9zs2bNNCCAgevTo4Vq2bGmvE8FQoUIFV7JkSRMTvO42bdrYax86dKi9bjwoiIe9e/eaoBJCZA0SCEKEGUTBrl273IIFC1yfPn1c9erVXZEiRSw8gJeAnfTnn39uu+h169ZZWOGXX36J28UQDwGL/bZt29z69estjwIvxPDhw90777xjAqFq1armMXnmmWdc6dKlLXSBJ6Jjx46Wc7F79+7A2YQQ0UICQYgwgKucxW/SpEkWh/d2yix8LIDjxo1zq1atcnv27PG9NyCaHDp0yDwsiCTeH0TVlClTTDzwPtaqVcu8Dbyf7733nvvuu+/s8UKIyCOBIEQGYVfLbvijjz4yFzohA0RBnTp1bGeMu5xdM8mHInSOHDlinhjeR0IQeBQQCy+//LL79NNP3Zo1a9yff/4ZeLQQItxIIAgRAsTYia8PGDDAYuq4xMuWLWtfk2OwePFi2+Eq6S68ELbZsGGDibEmTZrY+04eR+fOnS1pkxCNECK8SCAIkQaEA3B/EzPv2bOnJRUSI69cubJr3bq1Gzx4sFUkaCcbPXivCTX07t3b/g6IBbw2/fv3t7wOwhZCiMwjgSBEEOz8WWBIHiRvANc25XqIghdffNFK+sgz+OGHHyz5TmQtiDdyFjp16mRigTBEq1atLASxadMm/Y2EyAQSCEL8P8S7iWnTnIh6fjLq8+bNayWJHFu4cKHbuXNn4NHCb+DpQRBQIkqlSIECBVzBggUtN+Trr7+2EkwhRGhIIIiEBE8BCXC4qgcNGmSJbyQZEtcmxs2xuXPnWpKhiC2oFJk3b571X6CKhJJJOlESDqJHAxUnQohTI4EgEgp6DVB58O6771oJHaKARaR9+/Zu7Nix7vvvv1cpYpxAYuNPP/1kHoSuXbu6ihUrWkIpzZk+/PBDK6tU7ogQJ0cCQSQEzAPAI0CsukSJEq5w4cLmiiZ+TWWCBg3FNwi+zZs3W38Fwkb58uUzsfDWW2+5b775RlUQQqSCBIKIaxgW9PHHH1sIgeFCdDUcOHCgW7t2rbwECQqtqylH7devn3mPCCuRgEquyY8//hh4lBBCAkHEHXgDiDXTrAi3MmVwlCROmDDBvAVCeJBjwtApqlVowkTIiQoIxlcLkehIIMQolG9Riuf1qKc0j45zJGLRyOdUI3b5WXbWlPJRKubBAkrL4FiMzfKaeQ9ee+01K0vEiD0TWlA3Q5EWfN6pVGEq5dNPP22VLHPmzFGOgkhoJBBiFBr3MEIXFzolen379nW5c+d29957r3vsscesvCu1sjyy90naevbZZ12OHDnM2DXRAx9I0iNOjws2FkAoeR32uKizC2RkMrFmJaGJUEFI0pmxUaNGFpJCMJDUqn4KIhGRQIhByMxGAOAt4MKFx6BQoUKuQ4cOtgt6++233SOPPGKu0uQgGkjOYiGdNm2a1Y1TL05mN7XixGdJ3EIk+HnXTTUC45AJI/BamJT46quvuunTp7sdO3YEHiVExsAzx/8PIoGeGL169TIhKkQiIYEQg1COR1IViXZ4BLp16+YqVap0PNzAQs/9LVq0SLHIf/vttyYIEAYe1PzjjmfHDbhW6R6Ia95v0LsADwhhBEYoc/GmuyHtjlWJIMINnzdaOj/55JNW/UJYjjJYIRIBCYQYgx0+bk92+HgPGAyEax0j1AAcJ+mK/vTJL2aIC0IR8+fPDxxx7vPPP3f58+e3+D3goSDrH0+EX4YOkRcxYsQIa3iDMKCHAd8TYlE1gogkfL6SkpJc27Zt7bNXv359E6ma+SDiHQmEGIN2sgwMGjp0qH1PyIC6bmr6vTgpt+yqX3jhhRPc7Sz2eAvuu+8+C0V4TJ482XIRvvjiC/seoYEAYTHOyq5z/B7kRhDyoBqBjnhcpMk697wlQkQLPFQIA8Q4Ib1mzZqZR06IeEUCIcagfI+BNHgCABcogqFx48YnCAQSGKn55/5gSN576KGHThAIDB/KlSuXJT56kPRIrgLNZaINXhFCIIRIKFFEHNA2l7bISjoUWQ3VMiTFkp+AUOCzqf4JIh6RQIgxcHWWL1/+uEDgYkWNPyKB8AMgENhpk4mdvJ6bxETCCTNnzgwccW706NGuePHiNirXgyYyPE+0BAJJhwxLQsCwQ8Nb0LBhQ4v5El7wS6hDCA9EAf8nVAQhYhENe/fuDdwrROwjgRBjEHMnpMBCCiyc9JWndTAeAJISyS9g541rHrGAIR54LGOK2fW0a9fO2ssSgsD7wKhcbzARj+/YsaMlOka6YQy/F27aLl26mEjhYkv+BAmSCiMIv8P/1OrVq63fxlNPPWWiHE+XEPGABEKMQYVCmzZtXOfOnW3XDezyaRmbM2dOW2Dz5MljAmH58uV2P3HTN9980/IXiOuTn8BjWZBJuqIkktwEz31PsySSFBEYkdq5kzuBN4OES7LDqcLo37+/hVC81yVErMD/Fbkx3vTIkSNHypsgYh4JhBhkzJgx1gyIpkYehB6I2dMPgOQpduDe4j5+/HjL+qdvAJB4iCCgQRKeA9z4wZ0XGV5DmSP9FcLNli1b3CeffGLeCcQJ/Rfo5kjlhKoRRKyD8O3Zs6eV3xIikzdBxDISCDEIiyylVgMGDDhe2ggIAsIDyXf95BawOyc8EQxZ2cmT/nD5k+BI06FTtWsOBeK1iBIqKxAGiBh2XJrNL+IN/gcJ95HDg2Cn4ij4/1SIWEECIUahdwF5AiTwpQVigccQbkiP637ZsmWWnxCO8i0ulDwv3Q65WNKMqX379taISTXkIt5hxgn/owUKFHDNmze3ls1CxBISCDEKuQjE6081xx6BEEpp4Pbt2zPdlRDPAyEOKimIx5JfgNtVPe1FooEQpmQXgUzOz6hRo+x/V4hYQAJBhA36FxA2IJObUsoqVaq4YcOGpQhtCJFokEhMpUO+fPks/8ZLIBbCz0ggiEyBd4IEQxIhSZwkv4AOjPRW8MomhRD/rXSgW2mZMmWOVzogqoXwKxIIIkMQuli3bp0lSlIhwUWPKgpKF08V9hAikSFhl8mreNlI1qUCSQg/IoEgQoLchMWLF1ulA6KAdrP0ZKBSQv0LhEgfeBOodKADIyWRhOKUmyD8hgRCgsDFJzOVA5QjUn1Al0OaMeEi9cYsCyEyBvk5VPYwTlreBOE3JBDSAe70WDRg3PNnn31miVH0IfBEQvBj0oKfp6SSzop0PGSENLFTDacRIjzwPzllyhTzJjAgbcaMGYF7hMhaJBBOAd0EaRpEHTNDkWLF+J0ZAc3kx8cff9z961//slvuY8fCLa+Jr1OrzyaPgA6LlCjSY555DQgFEhKFEOEH0e1NMKU9uhBZjQRCGrDDpn7/wgsvdGeeeab7y1/+4ns7++yz7Xc999xz3cUXX+yyZ89urktex2OPPeb++te/2n08lsddcsklVpsNxEVJPKTrIm2YCSPQy4AdDZ4EIURk4f+MRmUlS5Y0z5/aj4usRAIhDRAIDCw677zz3GmnnRYTdsYZZ5igyZEjh+36u3XrZtMeGQ9NeACBEPx4HosgYK4DrxWPA4mHTFckGVGJU0JEF4Y80YGxYMGCNrckPaFAISKBBEIa8I/59ttvu7/97W8nLKp+NX5Pr3SK9sYMdaIfAQLh008/db1797ZQQ/KfoXcB3gKSD/E0rFmzRhclIbIQxqwTBmTqKiIhlG6oQoQLCYQ0iBWBQCiBLOjatWu7Xr16mWsSjwGhAxIKMb7GGJYU7BEhzMBAmY8++kiNjYTwEbt377aqITwJ/H8q3CCijQRCGvhdIFx66aXmMSB08MEHH5gwwFNA8yK8BSNGjLDdhycU+PrNN990d911l+UhcI7zzz/fDR48OPCKhRB+ghBfnz59LHGR/2HNMhHRRAIhDfwqEK699lrzGNDa+L333jNxwC15A8w/yJUrl7vvvvtc3bp17TjjZoO9CNWqVbPcA8510UUXmZAIBXYy2s0IER0QBQj+Rx55xEKAjGQXIhpIIKSBnwQCyYfZsmVzDz300HGPAW5HRMGDDz7obr31VnfnnXda1QKPIUnx+uuvd3fffbf1MCAfgVwEvAgkLt50003u9NNPN6FAF7dQ+O677+xCtXbt2sARIUQkwZOA9+/RRx+1MOLBgwcD9wgROSQQ0sAPAoEcgX/+85+26NPbAG/A8OHDbaEnz4DF3yvBpIyR8sSaNWu6GjVquPvvv9/yDfAqjBs3zrwIiAo8BlQrUNFAiGHIkCGBV5w+CElcc801FsrgPRJCRB4aKvXt29fCilyXGKsuRCSRQEiDrBQIZ511lvUowDvQqVMnCxWwwOMJYJEnVIBHgO5rPI6LBkmKhBW8W2bQ33LLLa5q1arHBQI/wzkYFnP11Ve7yy+/3L5PL8xiIHGK35EmSjt37gzcI4SINFQz8H+cJ08eq1RSuEFEEgmENMgKgYDL/8orr3RPPPGEa9WqlfUm8MoU8RrQ3dDzBLDzp0SR0sXq1aubMMB7gPE14oGExGCBgHEOQhTkMTz88MPu22+/DbziU/Pzzz9bJQS/K+dWxzchogsinfLl4sWLW4iRagchIoEEQhpESyDQ/ZAQAeIA1z/TEfEY0Np4woQJtrgTFqCrIWEGvAfjx483gVCoUCHzMjBrIVgg4EFgh3/bbbfZgh7sQeDiwvc8/vbbb3cTJ04MvOJTQ0MlQhkXXHCBPS9THenAKISILmwY8CTQMl0iQUQCCYQ0iIZA4NzkF9CLgGFIr732msX4EQiUN5G9jEuf3IKbb77Z3XHHHa5NmzbmASAXARFwww03WKMjRALCgNtatWrZ7AUSFRs0aGCCgAsKNnDgQJvTwPPiBZg8eXLgFZ+amTNnWkIkSY5NmjSxiojVq1cH7hVCRBOuAVQt0XmR+SlChBMJhDSIhkDAc0B1AmECqhDY0bPrZxH27MYbb7RqBMIOVClQ5oiYoPqAbmv8fuQTMFQJtyM102Q745GgVwKJjHgjEB2NGjVy+fLls/tJYsydO7f76quvAq/41KxYscK8FvzOeCSYEIlXQQgRfQg34GUk34i8ol27dgXuESLzSCCkQVYmKXrGwv/AAw9Y6KF+/fqucOHC7pxzznH33nuvhRjYObBY4ynglp09HgWEBWKCrwkxkM+AMKB1K2IDYcDtFVdcEVIfBC5IlFnxXFyYyKw+fPhw4F4hRFZALhCt0gk3aOKqCBcSCGmQVQKB57vuuussgbBChQoWLiBfgNABCYdUJtDngCFLuBh79OhhfQnIB+jatav1OWAR55ZERbwFiAWMUMY999xjz0EfBMocKZ0MBTwXeDXwSggh/MH06dNtA8GQti1btgSOCpFxJBDSINoCgee56qqrbEFHCJBs6IkDmiNxiyuRMASdEhEFVDdQ0UBlAuWK5BqQw0DogZGxhBF4LBULeBPofUBSpPecdFIMtVES+Q8SCEL4C65X06ZNs2mshBI3b94cuEeIjCGBkAbREgjMRSCUQAiAZESaHCEIMK+vQdmyZS3pEM8CuQXkI+BB8EofEQfkA7Ro0eJ46SNJiGQ5k7NAr4RgYeCZBIIQ8cWsWbNcmTJlbHPx448/Bo4KEToSCGkQaYFAl0OEAQs/qh9h4FUiYPyDIwzowY4wYJFHTLDQs/CzUFPuSIiA8senn37awg+UH5LZTB4CHQ9p05za82MSCELEH4sXL7bwIuGG9evXB44KERoSCGmQXoFA10OqEdJaiION1siEEljkS5cufbxvAZ4CjKZHlC3iAaC0kcciDBAUPA/nQARQBknOAXkFhBIQEpRBkqjIzyR/3tRMAkGI+GTZsmUWqiRJOSkpKXBUiPQjgZAG6REIdD0kFEBfARbm1B4TbCys5APgGUAMIAz4eTwH3JKJTBIhCzyLNxULNCXiNvg85BMQZiCR0RvMhJj4xz/+ccLjTmUSCELEL5QlFy1a1IyvhQgFCYQ0SE0g4CXAxY8rn8WcJkTeACQW6uDF1zMWd0IE7PrxGCAKPEGA0WwILwDljJQP/v3vf7c8A27xNqR2Tvob4C2g0RHeAx6f2uNOZRIIQsQ3S5cutYZqeBOWL18eOCrEqZFASAMEAiNWce17CyoVBFQHNGvWzEoMSRCkcoDJhizUwYsvQoJJjOzy6WPghRAQBYgEpiySUEiVAYmEeArIM/CEwcnEAYbnAoHAAp9a8mF6TQJBiPhn3bp1rl69epa8uGTJksBRIdJGAiENEAjdu3e3iYe4/ckZIJTAPAPmIXhTFRm7/O6771rjIRZdchIIEeARoCqBZEMvjPDSSy+ZkqdeGY8DCz29CHg8XgA8FPQnSL6QJzcel96ch7RMAkGIxGDDhg127aHb6ty5cwNHhTg5EginoF+/fpZb0LRpUwsj0F7YG7XsDT+ixJCOhuQAICZILqRfgScMPMNjQJtiShDpYMhjCT3goUAUpEcYhNskEIRIHDZt2mTeS0KaM2bMCBwVInUkEE4Bw5JYRNnxIwjoO8CtZ55AaNeuneUl4GXwhIGXgEj8j8REZi3gMaD0ECOkQOOicHgCMmoSCEIkFjRQ4rpEefWUKVMCR4VIiQRCGhBieOuttyzJkIUdLwFiIFggYIQY3nnnHUsYZBRyw4YNTRwgDAoUKGD9CMgroNIAz4EnDFJbsKNtEghCJB7btm2z6xS9U7788svAUSFORAIhDbwqBnoQ4P6nNwGCgNyD5AKBSYnkKZCHQHiBYUgsoiQpEqKg/JDF2Otj4BeTQBAiMdm7d6+Nl+e6xuA1BrEJEYwEQhokL3NkoX/99ddt/kFygfD+++/bhEUWXEoVyUdgEeVnL774YhMZJC96t94CnRV5B8EmgSBE4rJnzx7zjObKlctyrCQSRDASCGmQXCBQTkiZkFfF4AkE8hIYvUxIgZJGHsPYVRIbSQYif4HSSPomEIIgOZGFmdvMlCiGwyQQhEhs9u3bZzNcyJ8ihPrnn38G7hGJjgRCGiQXCNhtt93munTpYv0PPIGA8kY00BeBxZOvgycs4nGgDJIJa0899ZQlM3qdDyUQhBBZzW+//Wbj4imBZMNz9OjRwD0ikZFASIPUBAIJiyVKlDjeJMkTCRheBY4hDhAFo0ePNmFA74O8efNaSSOGOKC0MavFASaBIISA/fv3W9UWHk82OfIkCAmENEhNIGDkInTr1i1FRYMnEFDgJC3Wr1/fBihRucCCSiUESYp+EAaeSSAIITwOHz7sevToYZ4ENjnHjh0L3CMSEQmENDiZQKBEkUoFvAiIAk8YcIvHgBGrVDEwUAmjrJEuiX6rYMAkEIQQwRw8eNCuY3hKuaYdOnQocI9INCQQ0uBkAgHzvAjjx483T8KgQYPcyy+/7B599FETBbRO5ucQBsFVC34zCQQhRHIQCYyTL1iwoIkEhRsSEwmENEAg0CgpeFiTZ5QnMoCJfyI8BlQvkFtASSMeAyz5z/jRmAMxdOjQwCtOHxIIQsQ/iAQawFF9xSZIJZCJhwRCGiAQevbsaYsogoDcAc/wCpBbkDNnTputgDDgcV4ogWTG4Mdn1rzzeX0TvFvaNHNf8HHva6ZBet6L5Me9c9HAiSqMUJBAECIxILzQv39/CzeQZ3XkyJHAPSIRkEA4BQsXLrTyHxokEVLwjFJHkhCZwUDGL42SunbtasbjseDHp9cYL41Xom/fvicYQ6NoaJI9e3Zb1JnrgJcCcUINMyWTtHPmOJUShA7oy1CjRg17HOOkOc7XFSpUsOfiNeEh+f777wOvNn1IIAiROPz++++2USLcwLWOjZNIDCQQMkEkaoVR6PRGpxKCLGLPcPF5I6URAggDchzuv/9+EyP0Z0A40IiJcAdiYNy4cW7atGmWD8ExKiquvvpqO54ZJBCESCy4LrFRobqB8m31SUgMJBB8xuLFi61lM10WWYQ9IymS5EcWeEIb3rhovuc4JZRUVzD74YEHHnCXXXaZy5cvnw1jIfyBQHjwwQftcSQdZQYJBCESjwMHDpgnk+6w9HqRSIh/JBB8xqxZs2z3T2JQr169XPfu3a2SgmQhvqe/AuEMbgcMGGDxQf5puWWk9CWXXGJ5EYgG8iEwwguMdkU4SCAIITKKl7iISMDLqeqG+EYCwWcgEMgnoGRy/vz5bsaMGbYQs6ifzAhBEBssUqSICQTEQP78+e0We+KJJ0w85MiRw8INPDYzSCAIkbjQTIk8KTYxassc30gg+AxPIJAASf7A1KlT7WsWZDwLTIpMbngLuMVbQCXFrbfeasmKJDQGG56ESy+91NyDmUECQYjEhsRFLyeBDYeqG+ITCQSfESwQpk+f7r766iurRqCckRLFUxmli5Q2BpdIesY5EBF0fcwMEghCCHISCHsWLVrUrimqbog/JBB8RmoCoUqVKulu00xPBMIMeAqSG6WPJDeShZwZJBCEEEC4wWvLjEhQM6X4QgLBZ8yZMydTAoGuj+QZXHPNNVbhEGxUNhBmGDt2bODZMoYEghDCwws3kAOlcEN8kSGBQOYqiSl8EGThM95T8g4oacyoQDj33HPNS4BIoGIh2PAsUAYpgSCECCd0XGQKJKOiFW6IH0IWCLiQ6MJHmQsJKriWZOGxUqVK2bAndvpUMWREIPA4fp4eCcmNRkp4FiQQhBDhBpGAJ6FkyZLW7E3EPiELBHa6zzzzjPXzZ6dLxvwtt9wiC4PRDZGdP50SX3nllQwJBMQBHRMpb6Sromd58uSxBkwIBCUpCiEiwb59+9yrr77qSpcu7dasWRM4KmKVkAUCSSmFChWyBaJRo0auWbNmrmnTprIwWMuWLd1zzz1ni3yTJk0yJBDuuOMO17BhQ5vPgBfCM/5pK1eubH+3UIczJUcCQQhxMjZs2ODKli1r1509e/YEjopYJEMCgTgT3gPGHDdv3vyEhShcltoCGu8WDoFw++23W/4CAiH5uStWrGgLO1PZMoMEghAiLb7++mvzNFPhoMqG2CUiAoHjLHCZMc6RaMZ7iWsusx6EBg0apPAgtGrVyvopSCAIISINomDQoEEmEriGidgkQwKBAUBkw9epU8cWo5deesnsxRdftN1rp06dbGwxUwa9EcihGKOU+VmacPTu3dt6fyeCMV+B95EqBARDRgQCeQz8XQgzeH8XjPMhPihzVIhBCBFp9u/fbyForjurVq0KHBWxRIaSFPmDM2qYxYgdK25tjEQ7Jgbi0qb3P8OBHnroITsWqjF90JslEJxsF89GIiHvJ+WI7PozIhBotczfIfjvgvE9rZpJgqR/emaQQBBCpIfVq1e78uXL2wZl586dgaMiVghZINADYfjw4fYHD45xe9amTRuLd7PIsdATU6fsJRRDgDCqmLbACAUWSOLn8Wzea6TKgI6HKG9PIDz//POpioGTGe8dP08lRPDfhrAFyYrLli0L/DUzhgSCECK9cB2jjJt8BE1/jC1CFgg0wCC+RKghudHoh3GgjB4mBNG+fXtrHTx58uSQjAmGnIOdb7t27axcZsWKFW758uVxa95rLFiwoAkEr8yRYU1UizB8iZLSYMODwy39DYIFAg1LvEZWqf2dMvtPKoEghAiFoUOHWvUb13c1UYodQhYI6WHAgAEmENq2bWvZrBMnTgzJWBhRmwgEzkH8ikV06dKlcWsrV6601/jUU08dFwh0Vfziiy9snDN5A8mNqYx4c1DnDGnyBALNSiKJBIIQIhR+/fVXu6ZxrSLsIGKDsAsEdqeeQCDcgIuchSQUowsX88bZIXOOpKQkt2TJErdo0aK4NcQBr5HcDQQCYZrZs2ebB4H3ENGU3PDOICIIQQRPeySxM5JIIAghQmXt2rWWj0D4E8Eg/I8Egk8sWCBccMEFVhHC1EU8BXgJgm3YsGFmDEbhllIiTxxIIAgh/Mo333xjoQauUQo1+B8JBJ+YJxAIMZx11lnWZyJnzpwpqjuoCuE4Rktlvue9Pvvss639tQSCEMKvIAoo56ZUnrVB+BsJBJ8Yr3Hx4sUuf/78Jg6qVq1qjY2CrUKFCq5mzZqW/Nm5c2e7rVWrllVA1K1b1ypGmONAkmIkkUAQQmSUvXv3WiO3Z5991n3//feBo8KPSCD4xHiN3D7xxBNWyUCOwZQpU06o7uC9mTdvns1fP3bsmE1PW7BggSUyzpw50w0ePNjeM5pURRIJBCFEZti4caMrV66cq1evnvvll18CR4XfkEDwiXkeBAQCYQbeg0mTJrkJEyYct3HjxplAoHzRe6/nz5/vxo8fb4KC1qaUPUogCCH8DuXs5E+xXnAtE/5DAsEnFiwQSFT0PAZjx441YYDRAZHKBXpNAEKBpB/KIHk8vSPoovjWW2/Z/ZFCAkEIkVnIRxgyZIjN9qG8XUmL/kMCwSd2MoFAHwn6QKxbt87qhzdv3nx8Ohrv9ZYtW6x8iEZLhBhuvvlmeRCEEDEB5Y40gqODrvIR/IcEgk8suUDwwgs//vhj4J09Nbx3LNwSCEKIWIHND4nWTLPdvXt34KjwAxIIPrHUBALvBXMTaI+M14CQQmrGfQcOHHAff/yxVUAoxCCEiCXIR8ibN68bOHCgQg0+QgLBJ5ZcIHghBrwIVCiQa0D+QWrGfXPnznXvv/++chCEEDEHGx1G3lP6iFgQ/sD3AqF169bWRIjFc+HChXFrvEaEQrBAoHzxs88+syoFL1HxZEZJJA1IJBCEELEI5Y5MnC1Tpoz74YcfAkdFVuJrgcCwpo4dO1rNrJekF69GrgGJhgUKFDhBIKT2/qRmvGfyIAghYhkSFfEidOjQwfq8iKzFtwKBiYR4EOgUOGbMGIuvjxo1Km6N18jtAw88YI2S8AhIIAghEg2ufUWLFrW5M6wnIuvwrUDo16+flexdeeWV7t5773X33HOPu/vuu+PWeI133XWXO//88+2fQwJBCJGI4DmgEot5DZR5i6zD9zkI1apVs6mFqMmhQ4fGpXmTGfmnYOHlH0MCQQiRqJCPUL16dZs/s3Xr1sBREW18KxCYSHjHHXdYp61E4aeffnK5cuWyPAQJBCFEIkPiNvMaGEzndY8V0cX3AoFzJQokYzLGWQJBCCGcVXDRinnkyJHqj5AF+F4gsOglClQyxIJAYChUtmzZrARTCCEiBf0RXn/9dRtgR9M4EV0kEHxErAgEup3x90XdCyFEJNm2bZurW7euVbTxtYgeEgg+IlYEAg2diAuuXLkycEQIISIH15wSJUq4bt26mVdBRAcJBB/hd4HA5LUVK1ZooIoQIupQ6UUJOKFN5SNEBwkEH+F3gTB//nxXtmxZ1SYLIaIOA+nwXHINUj5CdJBA8BF+FwjTp0+3To8TJ060748dO2a3QggRDTZt2uQqVqzoGjRo4Hbs2BE4KiKFBIKP8LtAYKpk7ty53YgRI9zo0aMtuxhvgmKCQohoMXv2bFesWDGb/sh6IyKHBIKP8LtAYKT0/fffb9PWXnjhBffYY4+5PHny2PPrH1UIES3oQFu6dGk3bdo05SNEEAkEH+F3gTBv3jyXPXt216hRIxvHunTpUhss9dxzz7k9e/YEHiWEEJHl999/d23btrW29EyAFJFBAsFH+F0gzJo1yz388MNu3LhxgSPONW7c2DwJtIkWQohosWrVKlekSBHXvHlzVVZFCAkEHxELOQjMiuD54OjRo65Vq1b2+/788892TAghogXXpFKlStm1T2HO8COB4CNiRSB4LZYRCC1atJBAEEJkGf379z9+zVQ+QniRQPARfhcIVCxQ5jh27Fj7HoHQtGlTq2xQiEEIkRXs3LnT1a9f33Kh1q5dGzgqwoEEgo/wu0BISkpyL7/8sjVMAv7WgwcPdq+++qrbu3evHRNCiGhDPsLzzz9vHk0lTIcPCQQf4XeBwN+WuezcetAD4dChQ3LtCSGyFK6BVFVRAqnrUXiQQPARfhcIQgjhV1h7WDdKlixpPVtE5pFA8BGxIhCoQWbsKt4EIYTwC+QjNGzY0Jq5rV+/PnBUZBQJBB8RKwKB5ytevLj1RRBCCD+xevVqCzW0b9/eHT58OHBUZAQJBB8RKwLhvffec1dcccXxagYhhPATXJsKFy5sI6KVj5BxJBB8RKwIhA8++MDdeOON9rxCCOE3CH927NjRhjotXrw4cFSEigSCj5BAEEKI8MA46Nq1a9tgOfVpyRgSCD5CAkEIIcLHokWLXNGiRV2XLl2sHFuEhgSCj5BAEEKI8DJq1CgTCbSIVz5CaEgg+AgJBCGECC9UMnTu3NlEgvIRQkMCwUdIIAghRPjZsmWLK1u2rKtZs6bbvn174Kg4FRIIPkICQQghIsPs2bPNi9CjRw8bNCdOjQSCj5BAEEKIyDF8+HBXpEgRN378eOUjpAMJBB8hgSCEEJGD/ght27a1TrBMgBRpI4HgIyQQhBAisqxbt85VrlzZNW7c2GY3iJMjgeAjJBCEECLyzJgxwxUqVMgNGjRIoYY0kEDwERIIQggRef744w9bYypUqOCWLVsWOCqSI4HgIyQQhBAiOmzdutVaMb/88stuz549gaMiGAkEHyGBIIQQ0YPrZp48edyQIUMUakgFCQQfIYEghBDRgy6LTH2kiVJSUlLgqPCQQPAREghCCBFdmPRYrVo117JlSw10SoYEgo+QQBBCiOgzadIkV7hwYffJJ58o1BCEBIKPkEAQQojoc+DAAdepUydXunRpt3bt2sBRIYHgIyQQhBAia1i/fr3lInTo0EGhhgASCD5CAkEIIbKOsWPHuvz587sxY8Yo1PD/SCD4CAkEIYTIOvbv3+9effVVV6lSJbdx48bA0cRFAsFHSCAIIUTWsmbNGlemTBkLNRw5ciRwNDGRQPAREghCCJH1fPjhh65gwYJ2DU7kUIMEgo+QQBBCiKxn7969rkmTJq5mzZpu27ZtgaOJhwSCj5BAEEIIf7BixQr3/PPP21rEcKdERALBR0ggCCGEfxg5cqT1Rpg1a1ZChhokEHyEBIIQQviH3bt3uwYNGpgngemPiYYEgo+QQBBCCH8xe/Zs9+ijj7p3333XHTt2LHA0MZBA8BESCEII4S+OHj3qBg8ebF6EBQsWBI4mBhIIPkICQQgh/MeePXtcvXr1zPg6UZBA8BESCEII4U9mzpzpnnrqKQs1ZCRhkXwGRkuHoyICr8auXbvc4cOHA0cigwSCj5BAEEIIf8Ki/Oabb7oSJUq4pUuXBo6mDxZywhSMkz548GDgaMbhfJ9++qmbOnVqREWCBIKPkEAQQgj/8vPPP7vq1auHHGqYM2eOa9++vVuyZEngSErwSgR7Jk7lpaD0skWLFi4pKSlwJPxIIPgICQQhhPA306ZNs1ADPRLSU9Wwb98+17FjR/fee+/Z+rh9+3a7dm7atMnuxzPxzTffuEWLFp0QfkCMENbYsWNH4Mh/cyHmzp1r52BOxGuvveb69etn540EEgg+QgJBCCH8DQtz9+7dXalSpdK1e58/f76rUqWKCQu8AgyDqlixomvXrp07cOCAiYDnnnvOjRo16oSFfu3ate6FF15wgwYNOi5EJk+ebD87b948O9dHH33kXnzxRbdlyxa7P9z4XiAMHDgwcOb4B0UpgSCEEP5m8+bNrnz58u6VV16xuQ1pQe5BtWrV3Pr16+171shJkyaZaBg9erRr1aqVrZW//PKL3e9BbgEeAkIaeBHwLnTu3NnVrl3bEhTh22+/NXGBcDhVSCIj+F4gcK5EgQ8dDTmefPJJ+4NLIAghhD9hkSfUQOLhyRZnjr/99tuufv36xxd1wCPA9bpw4cLmAVi5cmXgnhPBU4AQId9gw4YNrkaNGhZS8J4PLwOjqfv3759YAqFv377u1ltvdUWLFrW53Cgpkjzi1YhRNW7c2F1zzTWuUKFC8iAIIYSPIdTQpUsXV7lyZVuoUwMhQOUD7ZqDBQIwUvqyyy5z5cqVs5yC1Pjtt9/MS9GtWzcLJ+CJWLhwYeBe53744QcTEH369EksgUCt6c033+wuvPBCd8MNN7hs2bLFtbHgXnfdde6cc85xxYoVs/IVCQQhhPAvGzdudFWrVjXXP2tfavTq1cvVqVPHeiB4rFu3zoQB+QRFihSxBf5kCY/jxo2zx9HJkXDE77//HrjHueXLl1vZ5YgRIxJLIPCG3XbbbeaaYbGcOHGiLZjxap7HgBwE3FYKMQghhP8ZO3asLfLjx49PdZHmfhb4ZcuW2fckFDZq1MhVqlTJ8hJIQnzmmWcsjPDrr79afwNaOnvn2rlzp/38XXfdZZUTwc/BdR8PRqh9GdKLr0MMCARe/JAhQyxZkfPGq/Eaub3vvvtcwYIFFWIQQogYgDLGZs2aWS4ALv/k4C2oW7fu8VyF6dOn27rG2sj3dFhs27at5d2tWrXK1axZ067hwUKAEHvx4sVTnJ91snXr1haKiAQxEWLABX/99dfHtfEayT84++yzFWIQQogYYsWKFZYL0LVrV+trEAzVB/RAYCGnt8H+/fvdtm3bTghJ4DmgUgGxwNrneQoIJ+AdwNvQqVOnE/okEN5o0qSJrRPBYiKc+N6DQB3osGHDzNhlx6sNHz7cFnZEkZIUhRAitqBkkc2d1+8gGCoQmjdvbo8JFgbJQSiwOeTxQG5B3rx57bzfffedHQOEAmvGG2+8YSGISOFbgeCVORKfSRTIZKXMkT4IykEQQojYgX4IJBGy2/cW+GDogEiewaFDhwJHUoL3gZAFIDIQC1S4ITqomvCgRwLtmwk5RMp7AL4XCJwrUcBlpEZJQggRmxBqyJMnj5U/Jg81AOtjOBZ0zhFJYeDhe4HAopcoqNWyEELELqx/JNSTUIi3IBqLeCSRQPAREghCCBHbkIRIcyMqFyKZHxANJBB8hASCEELEPsxIILGQZPtY9iJIIPgICQQhhIh9vNJGRgXQ9ChWkUDwERIIQggRH9DXoFatWmZ8HYtIIPiIWBEIQ4cOlUAQQohTwDTGkiVLmjeBpkcnm7fgVyQQfEQsCARqcRk+wgCtzz77zBp7BA8PEUII8V+4PjJ34cEHH7QWykzupR8C62QsIIHgI2JBICxatMg98cQT5kGgAyTtPzkmhBDif7Bxom0yg5wYOkgbZiYy4lGgEV4siAQJBB8RCwKB2l6eA4FQr149V6NGjVQHlAghRCLDNTl37txWyUC5Ix0SGdzUuHFja6e/Zs2awCP9iwSCj4gFgcC0MXqDn3feeS5Hjhz24Q8eICKEEInOwYMHzbuKx4Dr+p49eyxREZHAtZqR/qNGjfJ9CaQEgo+IBYGwa9cuy8o97bTTbD45/cCFEEL8D3IPmLTYoEEDEwdMcWSCI18vW7bMlStXzvXr18/3SYthFwgoIlpNhksgcK5EYdOmTe6RRx5xTz75pL0HEydOdBMmTEiXMcyD9wqB8PbbbwfOGBkQIKeffrqrVq2afeCFEEL8D7yqXIefeeYZm89Ad0WGOXGL5wAv7KRJk+LTg4DqwVOQ3DjOgAoWqmuvvda1a9fOzZw50xIyQrEZM2ZYckf27Nld//797bwne854MT4oeBAQCLifWPB5L/gQpcd4zwYPHuxuvfVW171794i8Z5wP8BBdfvnldsvvndpjZTKZLFGN6+KSJUtc+fLlXdOmTS13i/As454RDQ0bNjSvQrTWtYwKkZAFAk/GgvTOO++4Pn36WAw62HBzt2zZ8rgHgZ3w+PHjQzLc65yb3TBvMCIBd0zy54ono062Q4cO7vrrr3d33nmna9asmX2wcFOlx1q0aOGqV6/urrzySleqVCl7zxBZqT1XZozx25UrV3bnn3++uc9ompTa42QymSxRjWsv13S8rBdddJFt3HLlyuUuvfRSd8UVV9hYaNbKSK9rPXv2dCNHjrTwRkYIWSDgIeDFUdf56KOPuscee+y48X2+fPlc/vz53bnnnutuu+02G33J8VDs8ccfd/fee6+dg2x5zhH8PPFovGZe+4UXXmjx/bPPPtudc845dpse47Fnnnmm/Wy2bNnsfGTQpvZcmTHO+fDDD7u7777bvB2ReA6ZTCaLdWPd4vrIhg+BcPvtt9uml7WNUvFoXDsfeugh29AtXbo0sIKHRsgCAVfF5s2b3fLlyy22ktxwo0ydOtXVrl3bVahQIcOG56BKlSru448/dqtXr3bff/99qs8XL8b7tnLlSpsjTuYrf1RefyjGz1SsWNENGzbMzhep94zzJsLfRCaTyTJj3rWS6zHmfZ3aYyNhrNNr167NcDO7sCcpeoQjNu3FvBMNXjfvX0YsUd8zIYQQ4SViAkEIIYQQsYsEghBCCCFSIIEghBBCiBRIIAghhBAiBRIIQgghhEiBBIIQQgghUiCBIIQQQogUSCAIIYQQIgUSCAkKTZUOHz5sDamEyAq8QTXc0unt0KFDgXucfc3kO+bqe/CZpdU7t0eOHLH7+T453Mf5vOZhyR9z4MAB99tvv9nEPSEiBZ9trrGxjARCgsGFkXafDBLB9u3bF7hHiOjC57Bz585mJUqUsOE1v/zyi/voo4+sZTi97IsXL24Dwpivv3XrVpsDQztypuEx/ObFF1+0lrKAEGA4HC3HCxcubOd7/fXXbbCcJxqGDh3qnn32WZuY2rp1a/sdEBFChBuutUzY5fPLKP9Y7HIrgRDn8KHctWuXW7dunV08X3nlFRvkVKhQIZtXzvyHH374QSaLqK1Zs8bGmQd7CRYuXGhD3xhc06hRIxtx/tlnn7lKlSqZAGABZy5Lzpw53eeff+42bNhgA2huvvlmEwZME73//vtdrVq1TFjMnz/fPf30065YsWJ2PuaTMHa+bt26Nov/ww8/tPM1b97cRtFXrVrV7uO8wfA/s337duthz/9Naq9HJjuVMYehV69e9vkuWbKkiVOOMeY5VkSpBEKcw67pjTfecHfccYeN4D799NNt6mP27Nntg8sujOmMMlmkjAWeiXZFihRxy5YtC3wynVuwYIEJhDfffNN2+MCCjEeAi+iWLVvcJ598YiKAsbgIDD6zDILbvXu3XWT5bD/55JNu8eLFrmvXrnYhZmFnkef8CAqEBMcYglanTh03ceJEN3nyZPfWW2/Z78UuL/iCjYhh1DpT95iGl9prksnSMj7zXFv5vN50003ujDPOsEm9XHcZ5R8roQcJhDiHi93w4cPNa3DZZZe5Cy64wMaMskNjJ8UtF1CZLFLGZ4xdfps2bU7Yrc+dO9cVKFDAvAMeP/30kxs4cKBr3LixLeaMjr/66qvdkCFD7GcJDfTp0yfwaOfGjBnjypUr52bNmmU/07JlyxPyaurXr2+/A94KhAYeBUaVc6FmBC/CmVBEcD4CYqV///7mXUj+WmSy9Bqfe8QAn30+b+edd54JYrwKsZL/IoGQAPBhZGfGGOiXXnrJLqh8eEePHp3hMaBCZBYEAiGBSZMm2ffkEHDxJH+A3AFcsngF2I2RL7N+/XoTCDzGEwF8hhkNP3v2bLsYs/P3Lr48pmbNmnax/vbbb90jjzxiX0+YMMFCGePGjbOvER7KQxDhBqE5ZcoUE7p8RsmdWbRoUaqJtX5FAiHBIB+BXAR2WlyEd+zYEbhHiOgyZ84cCw+wSAOJiOy2yB0gB4C8Aj6jV111leXLIBDwKPTo0eO4QPj4449d6dKlTQCQDIZ3DK8CIQpExY033mhimPAE58WSkpJMGONVIIGRPBwhwg0J4H379rV8F8JdSlIUMQMZtlyEY70MR8Qu8+bNs7wETyBwAcXdT94AngUqGPLly+euueYa16FDB7dq1Sp7fO/evY8LhE8//dQ8YkuXLjXxS3IiHgcex3kIq+FVoFwSL0OpUqVMROClQJxw38aNG+1cQoQTPAgIVcpqYxUJBCFElrBt2zZb4Mn49tizZ48bO3as7brwGlCZQI4CXi8utjyehERvN0aZIuEC7mPHtmTJEjdq1CjbuXFL+WT79u2PCwqSJAlRvPbaa27EiBEmkoUQqSOBIITIEoj7e5ac5O5Yvk/t8d7X3E+ODT0QyEsgaZFQBCWPfJ+c5OcXQqREAkEIERfgPaBqgYodmiyRm0AJIz0QhBChI4EghIgbCFGQazB9+nQLXXihBSFE6EggCCGEECIFEghCCCGESIEEghBCCCFSIIEghBBCiBRIIAghhBAiBRIIQgghhEiBBIIQQgghkuHc/wESb5ZgfTXJtgAAAABJRU5ErkJggg==)"
      ],
      "metadata": {
        "id": "tIyVetMUZ8XX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "With an initial velocity of $50$ m/s and an initial height of $10$ meters, what is the angle $\\theta$ to obtain a horizontal range of $225$ meters? Give your answer in degree (not radian), correct to 4 decimal places."
      ],
      "metadata": {
        "id": "BrKsivKMeLGg"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "##**V. Gradient Descent/Ascent methods**\n",
        "The *Gradient Descent method* is a simple algorithm to find a local minimum of a function $f$. You start with an initial guess $x_0$ and successively update it using the recursion\n",
        "\n",
        "$$x_{n+1}=x_n-\\alpha f'(x_n)$$\n",
        "\n",
        "The *Gradient Ascent method* is to find a local maximum. It uses the same recursion formula above, except that the minus sign is changed to plus sign. The Python code to compute $x_1$, $x_2$,..., $x_n$ in the Gradient Descent method is as follows."
      ],
      "metadata": {
        "id": "RUrzT4qp7Ytx"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "def gradd(f,x0,alpha,n):\n",
        "  df = diff(f)\n",
        "  x = float(x0)\n",
        "  for i in range(n):\n",
        "    x = x - alpha*df(x)\n",
        "    print(i, x, df(x))\n",
        "  return x"
      ],
      "metadata": {
        "id": "UUfoQ-PVjoW_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 7**\n",
        "* Explain what each line in the code above does.\n",
        "* Use the code to find the minimum value of $f(x)=x^2$ with initial guess $x_0=1$. Stop when $|f'(x_n)|<0.01$. Give a value of the learning rate $\\alpha$ for which the sequence $(x_n)$ converges quickly to the minimum. Give a value of $\\alpha$ for which the sequence $(x_n)$ fails to converge."
      ],
      "metadata": {
        "id": "lv3_kr_4jmgu"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 8**\n",
        "* Graph the function $2x^4-9x^3+18x-4$ on the interval $[-2,4]$.\n",
        "* Use the Gradient Ascent method to find its local maximum. Choose your own initial guess and learning rate. Stop when $|f'(x_n)|<0.005$."
      ],
      "metadata": {
        "id": "4SX6-2mUntYO"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "###**Exercise 9**\n",
        "If you take a loan of $P$ (dollars) at a mortgage rate APR$=r$ and if the loan term is $n$ months, then your monthly payment is\n",
        "\n",
        "$$q=\\frac{Pr/12}{1-(1+r/12)^{-n}}$$\n",
        "\n",
        "You might have seen this formula if you took Math 107. Imagine that you have found your dream house that is listed for \\\\$400,000. You can only afford \\\\$2,000 per month. Suppose the loan term is 30 years. How low must the mortgage rate be (in percentage) for you to afford this house?\n",
        "\n",
        "*Hint: you can model this problem as an optimization problem or as a root-finding problem.*"
      ],
      "metadata": {
        "id": "WlszgpYxpTQH"
      }
    }
  ]
}