Lesson 3 - Numerical calculus.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lesson 3 - Numerical calculus"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Authors:**\n",
    "* Yilber Fabian Bautista\n",
    "* Keiwan Jamaly\n",
    "* Julia Lienert \n",
    "* Sean Tulin"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lecture, we will cover methods for calculating derivatives and integrals numerically. While we will cover some of the theory behind these methods, the main practical takeaway is to gain familiarity with existing functions within the `numpy` and `scipy` libraries that implement these tasks. \n",
    "In Python, most basic tasks have been coded up already in a very efficient way, so knowing your way around Python's libraries saves a lot of work. \n",
    "We have already introduced `numpy`. `scipy` is another important and vast library with many algorithms for scientific computing.\n",
    "\n",
    "By the end of this lecture you will be able to:\n",
    "* Do **numerical differentiation** and understand the principles behind it.\n",
    "* Perform basic **numeric integration** using python integration libraries\n",
    "* Use several integration techniques including: **the rectangle, trapezoidal, and Simpson's rules**. \n",
    "* Do numerical integration in 1 dimensions using quadratures using `scipy` (**quad**)\n",
    "* Perform higher dimensional integrals using `scipy` **dblquad, tplquad** "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Differentiation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Differentiation is a fundamental tool in science. Naively one would expect the usual derivative definition will  be used for numerics, that is\n",
    "\n",
    "\\begin{equation*}\n",
    "\\frac{f(x+h) - f(x)}{h}\n",
    "\\end{equation*}\n",
    "\n",
    "Expanding  $f(x+h)$ around $h=0$ we verify, that it is indeed the differential of $f(x)$.\n",
    "\n",
    "\\begin{equation*}\n",
    "\\frac{f(x+h) - f(x)}{h} = \\frac{f(x) + hf'(x) + \\mathcal{O}(h^2) - f(x)}{h} = f'(x) + \\mathcal{O}(h)\n",
    "\\end{equation*}\n",
    "\n",
    "As we can see,  the numerical  errors propagates with $\\mathcal{O}(h)$. This is a good approximation to start with, however  we can do better by simply changing our definition of the differential  and expanding around $h=0$\n",
    "\\begin{equation*}\n",
    "\\frac{f(x+h) - f(x-h)}{2h} = \\frac{f(x) + hf'(x) + h^2/2 f''(x) - (f(x) - hf'(x) + h^2/2 f''(x)) + \\mathcal{O}(h^3)}{2h} = f'(x) + \\mathcal{O}(h^2)\n",
    "\\end{equation*}\n",
    "This is called the central difference formula and is primarily used for any numerical differentiation, where now the numeric errors propagate as $\\mathcal{O}(h^2)$\n",
    "\n",
    "## Exercise 1\n",
    "Using the techniques from above, prove that \n",
    "\\begin{equation*}\n",
    "\\frac{f(x+h) -2f(x) + f(x-h)}{h^2} = f''(x) + \\mathcal{O}(h^2)\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Differentiation with scipy and numpy\n",
    "\n",
    "The above numeric differentiation techniques are already implementations in `scipy` and `numpy` libraries. The two functions have different applications.\n",
    "\n",
    "* [scipy.misc.derivative](https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html) for calculating derivatives of *known* functions. That is, you can know what $f(x)$ is and can calculate it at any point $x$, and you want to calculate its derivative at a specific point.\n",
    "* [numpy.gradient](https://docs.scipy.org/doc/numpy/reference/generated/numpy.gradient.html) for calculating derivatives of *unknown* functions, where you know $f(x)$ only on a grid of specific points $x$. \n",
    "\n",
    "Both functions use the central difference formula.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example with `scipy.misc.derivative`\n",
    "\n",
    "SciPy has the function [scipy.misc.derivative](https://docs.scipy.org/doc/scipy/reference/generated/scipy.misc.derivative.html), which allows us to  directly compute the derivative of an analytic function. The syntax for using `scipy.misc.derivative` function is the following:\n",
    "\n",
    "```py\n",
    "scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)\n",
    "```\n",
    "\n",
    "where `func` is the given function, `x0` the point  the derivative is calculated at, and `dx` is a float for spacing, in our previous formulas `dx=h`.\n",
    "\n",
    "Let us see an specific example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4.999999999921734\n"
     ]
    }
   ],
   "source": [
    "# include the function 'derivative' from SciPy\n",
    "from scipy.misc import derivative\n",
    "\n",
    "# define the function of which the derivative should be taken\n",
    "def test(x):\n",
    "\n",
    "    return x**3 + x**2\n",
    "\n",
    "# calculate the derivative, store the result of the calculation in the variable 'result' and print it\n",
    "result = derivative(test, 1.0, dx = 1e-6)\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From analytic differentiation, we  would have expected $test'(x)=3x^2+2x$, and evaluating at $x=1$ will produce $test'(1)=5$. The difference with our numeric differentiation is then"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-7.826628234397504e-11"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result-5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see this is a number approximately of order $\\mathcal{O}(dx^2)\\approx10^{-12}$ as expected from the central difference formula."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example with `numpy.gradient`\n",
    "\n",
    "Now, suppose we do not having a known function $f(x)$, but only data points are known on a grid of points in $x$. \n",
    "The derivative at each point is calculated using the formula of the centered derivative. Since both the leftmost and rightmost point do not have two neighbouring points, their derivatives are calculated using the forward (for the leftmost point) and the backward (for the rightmost point) derivative.\n",
    "\n",
    "`numpy` includes the function [np.gradient](https://numpy.org/doc/stable/reference/generated/numpy.gradient.html) that does the job for us, provided an array with data points is given. The syntax is\n",
    "\n",
    "```py\n",
    "np.gradient(f,dx)\n",
    "```\n",
    "\n",
    "where `f` is the array with the data points and `dx` is the spacing between the data points.\n",
    "\n",
    "Notice that the `numpy` functions for taking derivatives only work for data points given as arrays but not for specific functions. For the latter case, one has to use `scipy.misc.derivative`.\n",
    "\n",
    "Also note that `f` need not be a one-dimensional array or list. If $f(x)$ is a vector function, evaluating $f(x)$ on a grid of $x$ points will yield a 2D array. Then `np.gradient` will return an object of the same shape, representing the *gradient* $\\nabla f(x)$ at each grid point $x$.\n",
    "\n",
    "Let us see an specific example of how to use the `np.gradient` function in practice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAEKCAYAAADTgGjXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA0EElEQVR4nO2de7gcVZXofysnJ48DMSEPQgjkJIyIPJRHMkEBkyAaIMhDxe+CIQqJZPCKg4KOkVxExFwZVGRwHL0B5HWOIDMiIk8RCTIyjBycAEGEREhCIJDwioS8T9b9o6qSOt1V3VXVXV3V1ev3ffvrrl27undX79pr77XWXltUFcMwDMNIQr+sK2AYhmE0LyZEDMMwjMSYEDEMwzASY0LEMAzDSIwJEcMwDCMxJkQMwzCMxPTPugKNZuTIkTp+/Pisq2EYhtFUPP7446+p6qjS/JYTIuPHj6enpyfrahiGYTQVIrIiKN/UWYZhGEZiTIgYhmEYiTEhYhiGYSTGhIhhGIaRmMyFiIj8VETWiMiSkPMiIleJyDIReVJEDvOdO05EnnXPzUurjt3dMH489OvnvHZ3p/VNBcJummHkgtQfRVXNNAFTgMOAJSHnZwD3AAJ8APhvN78N+CuwDzAAeAI4oNr3TZw4UePQ1aXa0aEKO5OI8zpihJNEVDs7nbItR1eX8+NFdt4Q/02ym2YYDSHqo9jRkeyxA3o0qI8Oymx0AsZXECL/Dzjdd/wsMAb4IHCfL//rwNerfVdcIdLZ2fcPqJS8P6vwfaPXWoNaaNyUtEUbhpH4UezsjP9dYUIkc3VWBMYCL/qOV7l5YflliMhcEekRkZ61a9fG+vKVK6OXVXdrlhUrYO7cgmlwvDmxCMya5fxI2Pmjk7JhA5xxhqm8DCMm3d1OP5PkUYzTr1WjGYSIBORphfzyTNWFqjpJVSeNGlW24LIi48bFKr6DQvWNtbTWqKxY4QgnkYLcNMNIB288d8YZTj+ThKT9WhDNIERWAXv7jvcCXq6QX1cWLICOjuTXN3XfWI/WGodCT+UMIzlhioAkdHQ4/Vq9aAYhcgfwGddL6wPAOlVdDTwG7CsiE0RkAHCaW7auzJwJCxdCZ6dzLEHznyo0Zd9YOvtIgnezkty0Qk3lDCM59VAEeI9gZ6fTn82cWb/65cGofjOwGtiKM7uYA5wDnOOeF+BHOJ5YTwGTfNfOAJ5zz82P8n1xDeulRPWAqLdRqyH4f1xbW3yPgjDPq7CbZsZ3w6hKEueeNJwgybN3ViNTrUIkDL+XRNQ/O1deXEG+zGm6osX9vpZwezOMncTtU9J+PEyIpCxEPOL2jZkPsuO21Hq21iT+iZnfMMNInzj9SKMeCRMiDRIiqsn6xkwG2XmSeHGEWW71gYZRG1EfgyzWpJkQaaAQ8eM3B+RlRLGDKK21ra2xureogs1UW0bByHvTDxMi4pxrHSZNmqRZbUo1fnx1Z6fOTli+POWKdHfD/PnVK9PRkYIrRwTyXj/DSIHc9A8hiMjjqjqpNL8ZXHwLQ5Q1JytWpOzVGtV1NxVfwIjMnOk8KV1dlW/Yhg2OsDGMJsZbAxJlzFTP9R31woRIAyldcxJGqutJ5s+vvHCwo8PpvJcvz36EH+WGrVhhkYKNpqUZxnRVCdJxFTk12iYSRhT9Z13tx1Esdnm2M0Sx35jnltFkVGvWeWrSNHEAxkISdZBdlwF2lOGOp2zN5VCHaLpAU28ZTUIUFVauZx8+TIhkiKf6ryZIalZtRVFh5VHZ6scvdSuFUUndqGQYtVGEMZ0fEyI5oNogu+YBdqW4z80y3IGdUnf79gZIXsNIhyKM6fyYEMkBUVRbieL/e3PmMDfuZhrulJK65DWMdCjKmM7DhEhOqKbaUo2ppak2Z2624U4pqUlew0iHoo7pTIjkjEoD7Fhamkpz5mYc7gRRd8lrGOlQ5DGdCZGcUW2AHVlLEzYKF2nO4U4l6iZ5DSMdijymMyGSQ7wBdpgTUkUtTbU5cz33xcwLdZO8hpEORR7TmRDJMWH9faiWpshz5mrUJHkNIx26u52xTZHHdJkLERE5TkSeFZFlIjIv4PxXRWSxm5aISK+IDHfPLReRp9xz2URVTJHYWpoiz5mjElvyGkY6eGO6sPFLYcZ0QcvYG5WANpytbfcBBgBPAAdUKH8i8Dvf8XJgZJzvzEvYk6hUi1bSJzRKWLx5kYxqnwHV4snkKY6EUWiqPbfN1gzJadiTycAyVX1eVbcAtwAnVyh/Os6e7C1DJC1NK9pBwjD7iJETimwH8ZO1EBkLvOg7XuXmlSEiHcBxwC982Qr8RkQeF5G5YV8iInNFpEdEetauXVuHajeecC2NMn7Wh+hecURwgcLMmWNg9hEjQ1rBDuInayES9JSH7ZJ1IvAHVX3Dl3ekqh4GHA98QUSmBF2oqgtVdZKqTho1alRtNc6IcPuIsELHMZer6eb0vqdaxQ4ShtlHjAbTMnYQH1kLkVXA3r7jvYCXQ8qeRokqS1Vfdl/XAL/EUY8VkqpaGnZhPv93Z0bR5sxJsPUjRoNpRd+WTLfHFZH+wHPAMcBLwGPAp1X16ZJyQ4EXgL1V9R03bxegn6q+7b6/H/iWqt5b6Tuz3B63XvTrFzxVFraznTbnIMt9NPNEta127T4ZdST02RQnbmgzk8vtcVV1G3AucB/wDHCrqj4tIueIyDm+oh8HfuMJEJfRwH+KyBPAH4G7qgmQojBu+PrAfEUYzwt0t59ZvDlzUsw+YqSM59fSr0JvWjQ7iJ/+WVdAVe8G7i7J+0nJ8fXA9SV5zwMHp1y9XLKAC5nLd9jALiVnhBWMZ65cDfSnYLPm2hg3Lng24tlHFiwonp7BSB3PBlKk0O5xydomYiRg5hv/ykLOppPlBPkhbNjS37xYSzH7iJECYTaQtjZn8ltUO4ifTG0iWVAEmwjjxsGLjmd0P3oJWu5TBB1s3TH7iFFnimwDKSWXNhEjJp7y9cWdS2vGEazTL7IONjFmHzHqTNhz1krPnwmRZiEouKIIC7iQDimfT9tW4xWwJ9+oEW88FzSpLboNpBQTIs1CkPJVlZmdj7Dwpo7A9SOm6g8hzD5ikteIQKVg2a1gAynFbCLNQgTla9jIyFT9AVSyj3R0tF5PYESmVZ8zs4k0O2PGBOf7VDBhKn1T9QdQaWtdC9JoVMCes76YEMk7XjS3lwOiwZQoX03VnwDrEYyY2HPWFxMieSYompvnWRSgfDVVfwKsRzAiYsb0YEyI5JkQY/oO5WuJzr5SkEYzsodgkteIgBnTwzHDep6pYSVTqxr/EmFGdqMK9jyZYb052WOP4PwIqhZT9cfAjOxGFex5CseESB7xjOmrV5efi6h8NVV/AqynMEKw5ykcEyJ5I6YxPQxT9SfAegqjBDOmV8eESN6IaUwPw4zsCTDJa/gwY3o0MhciInKciDwrIstEZF7A+Wkisk5EFrvpG1GvbUrqqFIxVX9MTPIaPsLCvMcczxWerLfHbcPZHvejOPutPwacrqp/9pWZBnxFVT8W99ogcu+dtdde8NJL5fk1uIG0UrjqumHuOC2PPTd9yat31mRgmao+r6pbgFuAkxtwbf7wlK9BAqRG5aup+hNgRvaWx56baGQtRMYCL/qOV7l5pXxQRJ4QkXtE5MCY1+afkDDvQF2Ur0Gq/vZ2MwpWxHqQlufii8u3njFjejlZC5Gg3YFKJ5B/AjpV9WDgh8DtMa51CorMFZEeEelZu3Zt0rqmR52M6WH4Vf0i0L8/bNsGs2aZvTiUIMk7aJD1IC2ApxSYPdt5DIcObZ2tbpOQtRBZBeztO94L6BNpUFX/pqrr3fd3A+0iMjLKtb7PWKiqk1R10qhRo+pZ//rQANWJZ2S/6SZHiKg6yezFIZRKXnAU4SZ5C02QUmDrVue5MWN6MFkLkceAfUVkgogMAE4D7vAXEJE9RJynWEQm49T59SjXNg177RWcn4LqZP582LSpb555aoXgl7zt7bBli0neghOkFLDnozKZChFV3QacC9wHPAPcqqpPi8g5InKOW+xUYImIPAFcBZymDoHXNv5X1IGjjirPS0n5avbiBMyf7wxH/VjPUkjs+YiPBWDMku5umDcPVq1ydExDh8IbbzgzkAULUpk7m+dqAszXs2Ww5yOcvLr4ti6e8nXVKud42zbYuDF15astyk6AeWoVHgtvkhwTIlmRkfLVFmUnwHykC03KHvaFx9RZWZEDFYlN3WPg7TmycqXz33n/UYqqR6Mx2HMQDVNn5Y0a9gqpF2ZEjEGpp5b5SBcGew5qw4RIVgStV2mw8tVU/QkwH+nCYc9BbZgQaTTd3c4s5MknYfBgGDEis+WwpupPgA1bC8dFF5XnmTE9OiZEGolnwXv1Ved448aGeGSFUboou73dwqFUxYathcHzyPrc55xjC2+SDBMijSSHy2H9qv62NlP1V8V8pAuBhTepH+ad1Uhy4JEVhnmoxMDz1ApbVGDD2Nxj7T0+Yd5ZJkQayciR8Prr5fk5aLk5lm/5xXqipsXae3zMxTcrPMVrv35OSJOcblBgqv4EmJG9abH2Xj9MiKSJX/HqGRv69cvUIyuMIFV/W1su5Ft+sZ6oaTnvvPK8nIznmg4TImkSZEjv7YVdd3XmzDmy4JV6ag0Z4lT1q1915J7ZjAMIkrwDB1pPlGM8xcD55zvHu+2Wu/Fc02FCJE2aTN3heWpt3w6XX+7krV5t3lqhBG1c1dtrPtI5Jcgja/Nm88iqFTOsp0kTG16buOrZ0N0NZ53Vd98R89TKFdama8MM61mwYIGzgs9Pkyhem2wSlT22cVXusTadDpkLERE5TkSeFZFlIjIv4PxMEXnSTY+IyMG+c8tF5CkRWSwiOQjN6+IpXmfNcjqW9vamU7yazTgm1kPlHmvT6ZCpEBGRNuBHwPHAAcDpInJASbEXgKmq+n7gUmBhyfmjVfWQoGlWJpR6ZIHj5tRkitcgm3GTTKKywXqo3OOFN/Fjbbp2sp6JTAaWqerzqroFuAU42V9AVR9R1Tfdw0eBvRpcx3gEeWRt2tR0ao2gzava281mHIp5auUWTzHgBVocPrzpFAO5JmshMhZ40Xe8ys0LYw5wj+9Ygd+IyOMiMjeF+sWnQGoNz1vr2mud43XrzFMrlFJPLRGLZpkDgjyyNm1qOsVArslaiEhAXqC7mIgcjSNEvubLPlJVD8NRh31BRKaEXDtXRHpEpGft2rW11rkyBVRrfOtb5XlmMw6gdOOq3l6TuhmTw5inhSNrIbIK2Nt3vBfwcmkhEXk/cA1wsqruCD6lqi+7r2uAX+Kox8pQ1YWqOklVJ40K2gyqnjSxR1YYBZpcNYb582HLlr551nNlgrXd9KkoRETkU+7rhJS+/zFgXxGZICIDgNOAO0rqMA64DZilqs/58ncRkSHee2A6sCSlekZnxgzHkD54cGEUrwWcXKWL9Vy5wdpu+lSbiXzdff1FGl+uqtuAc4H7gGeAW1X1aRE5R0TOcYt9AxgB/FuJK+9o4D9F5Angj8BdqnpvGvWMhGe9Gz7cUbpedFHuQpskxTy1YmI9V2747GfL86zt1hlVDU3A/cCDwJs4M4Q+qdK1eU0TJ07UutPVpdrR4YVYdFJHh5NfELq6VDs7d/68oUNVRZy8Av3M+hDUHgYMsBvVQErb64gR1l5rBejRgD61YtgTV8V0GHATUOZlraoP1VuopU0qYU9aKJ7CddfB7Nl98yy6RwDexlUrVzrqzdGjnff9sjZDFh/PI8tvULc2Wjs1bUolIqNUNdStSUR+qKpfrLGODSEVIdJCO9y0kLysH7fcAqefvnNTsnHjHH2K9WipYG00HWqKnVVJgLgcmahWRaGFdOBmM07A1q3OgOK118zltwFYG20sNreuBxdeWJ5XUOtdC8nL+nHRReUzVXP5TQ1ro43FhEgteB5Z//APzvGwYYVx6w3DPLUSYEPjhnLmmeV51kbTo15CJGjlebEJiqewZUvh4ylYTK0E2NC4IXhjuksucY5zuAt1IYkkRETk0CpF/qUOdWkuWjiegsXUiolN31InaEy3cWPhx3S5IKp31oPAGODfgVtU9em0K5YWdfPOaiGPrDDMCyYGnsuvd8OGDIH1681Tq05YW0yfWr2zjgamAWuBhe5GUP+nvlVsMkxFYar+OHjTtxtvdAYab79t07c6Ym0xOyLbRFT1FVW9CjgHWIwTjqR1+UbAz28xFYXJ0QSYp1YqWFvMjqg2kf1F5JsisgT4V+AR8r45VFp41rs5c5zjoUNb1npnqv4E2JA5Fc4+uzzP2mJj6B+x3HXAzcB0dcOvtyRB8RS2bnWsdy0kPDy8n+xX9V94YUveiuiMGxesvLchc0386U9O4OwRI+Cll8zU1Eii2kQ+oKr/0tICBFraIysMT9X/5pvOQ7xggeNzYO6+Idj0ra50d8Oee8Jtt8GAAXDZZYUJnt00VNtP5Fb39SkRedKXnhKRJxtTxRxhqohQ7rrL2Q1240azF1ckaKFN//620CYBnmJg9WrneN06a3NZUC2K7xhVXS0inUHnVTVgXp5vanLxNT/CUOzWJOCGG+Css/oa2i3cbGSszTWWRC6+rgBpA65V1RWlKbXa5pVLLnGM6H5MFQHYJC0RF19snlo1YG0uH1S1iahqL7BBRIamUQEROU5EnhWRZSIyL+C8iMhV7vknReSwqNfWDc8j68wznYf+Xe9qWY+sMMzFMgHWC9bEmDHB+dbmSvD6r5SMlVG9szYBT4nI/cA7Xqaq/mMtX+7Ocn4EfBRYBTwmIneo6p99xY4H9nXT4cCPgcMjXls7QR5Z27a1rEdWGAsWlN+mQYNsklYR89Sqib/7O3i5xNXHFAMllPZfnrES6tZ/RV1seBdwEfB74HFfqpXJwDJVfV5VtwC3ACeXlDkZuNHdofFRYJiIjIl4be2YR1Yk/PZiT+OnavbiipinViK6u2HsWHj4YRg40AItVqQB/VdUF98bgFuBR1X1Bi/V4fvHAi/6jle5eVHKRLkWABGZKyI9ItKzdm21/bVKMJVDZDx335tuchyONm82T62KWEjk2HgDa28GsnmzBVqsSAP6r6gr1k/ECXVyr3t8iIjcUYfvDwohX+ouFlYmyrVOpupCVZ2kqpNGjRoVr4am7I/N/PmOxs+PTd5C8CTvNdc4xxYSuSKmGIhJA/qvqOqsb+Koj94CUNXFwIQ6fP8qYG/f8V5A6YLGsDJRrq0dUznExiZvCbj00vI86x3LsLYVkwULnNmtnzr3X1GFyDZVXVeSVz2GfHUeA/YVkQkiMgA4DSid4dwBfMb10voAsE5VV0e8tnZKlf2meK2KTd4SYL1jJPbcMzjf2lYIJ57o6JYHD06t/4oqRJaIyKeBNhHZV0R+iBOEsSZUdRtwLnAf8Axwq6o+LSLniMg5brG7geeBZcDVwP+udG2tdQrEUzlYPIVIBE3evJAoRggmeSNx0EHleaYYCMBz6x061DEaXXhhev2XqlZNQAewAGf03+O+HxTl2ryliRMnqpE+XV2qnZ2qIqqgOmiQ876z0zlnlNDVpdrR4dwsL3V02M1y6epSHTvWuS0DBqiOGGHtKZSU2hLQo0HyISizyMmESGPp6lJta7O+MRKe5PVu1OWXZ12jXGDyNSb+NuRPnZ01fWyYEKkWO+vXVLB9qOpJdZwUNYS6bY9rRMLiGyXgtdechRDt7Y5xvcXjmlsbiklKW3eHxc6qtmL9e+7rJ4A9gC73+HRgeeLaGC2D2YsTcN99zsP+jhscIoVVxs2EtaGYNDgSQrUAjA+p6kPAoar6v1T11276NHBUKjUyCoXZixNgC236MDZwCbG1oVAuvrg8L0Xvg6jeWaNEZB/vQEQmADFX7RmtiHlqJcCG3n047LDyPPPICsDzyJo92zluUKDYqELky8AiEVkkIouAB4HzUqmRUSiCYmqBRfaoiE3f+gSe/fWvnWCLtlSrAl48GL8aywsUm/ayhCBre1ACBgIHu2lg1Ovylsw7KzvMUysiQe5Igwe3zI0K+vmDBrXMz09GSh5ZfkjinVVEzDsrO8zLJgbd3Y4NZOVKpzsYOBC2bGkJTy1rJwlIySOr70cl2NnQMOqJqfpj0MIhka2dJCBDFWhFISIiUTetMoyqmKo/AS3oqWXtJAEZbt1dbSbyqIjc7sayGp96bYxCY55aCWjBYXkDAs8Whxxs3V1tncgkdnphXSkij4nID0RkuogMTL12RqEwT60EtNCw3OsPZ82CrVsdQWLeWBXI0iPLT5C1PSwB7cCHgcuBPwJ3xbk+D8m8s/KBeWpFpEU8tcwjKwEN8MjyQxreWSIyVlVfqocwaxTmnZUPzAMnBqWeWoMGOYb2AnlqWXtIQAM8svp+bAreWc0mQIz80IKq/uSUempt2lQ4Ty1rDwnIiarTXHyNTMhJ+28uCuypZe0hARl6ZPmp5uJ7k/ta9xAnIjJcRO4XkaXu624BZfYWkQdF5BkRedpfDxH5poi8JCKL3TSj3nU00sM8tRJQ4OG6eWTFIAceWX6qzUQmikgnMFtEdnM7/h2pxu+eBzygqvsCD7jHpWwDLlDV/YEPAF8QkQN853+gqoe46e4a62M0kCBPrWOPLYR6Pz0KPFyfMQPa2lLdCrwY5MUjy0c1IfIT4F7gvcDjJalW6/TJwA3u+xuAU0oLqOpqVf2T+/5tnL3UQwJDG82Gf+v6yZPhV79ybIXm7htCAadv3qB6+HDH1HPRReltBV4I5s93VJh+slZpBrlslSbgx1HKxUnAWyXHb1YpPx5YCbzLPf4mzsZYTwI/BXarcO1cHKHXM27cuJpd3Yz60tXluHOau28ESjevHziwaTcbt21vE+D976VJJPWvJosAjCLyW5wdEUuZD9ygqsN8Zd9U1TK7iHtuV+AhYIGq3ubmjQZew9m+91JgjKrOrlYnc/HNH+bemYDubkcn7je0d3Q0lQ7I/vcEZHjTwlx8M4viKyLPAtNUdbWIjAEWqep+AeXagTuB+1T1ipDPGg/cqaoHVfteEyL5o8Hu7sWgAD2w/e8JuO66nZtOeTRo8JDHKL53AJ91338W+FVpARER4FrgmVIB4goej48DS1Kqp5EyBbYXp0cBPLXsf49BRrsWRiFLIXIZ8FERWQp81D1GRPYUEc/T6khgFvDhAFfey0XkKRF5EjgaZ/dFowkJshf379/U9uL0KUAP/NWvlueZW28AOfTI6kOQoaTIyWJn5RO/vXjgwJ22wia0FzeGJg425f3XXrV3283+64o0OEZWGIQY1m3FupEL/JE9+rmtUosV2aO+BC202b499yGRgwbVmzfnZ1CdS3KuurTtcY1cUQB7cePp7oazznLip3vk1FPL/t8E5OSm5dGwbhhl5HzQlU/mz+8rQCD7BWgh2P+bgHkBwTxyZDwyIWLkigLYixtPE/XM9v/GwPPI+vznneNhw3LjkeXHhIiRK4I8tcCZzedY1Z8tTdAze/1hkFYmR4Pq/BBkPNqyJZfGIxMiRq7w24tLMSN7CEGSd+DA3PTMQf2h5wuQs0F1fshjjKwQzLBu5Jac2BObA//uh+C4uG3fnovdD+1/TEAOl/ObYd1oOppI1Z89fh/p9nbo7c2Nj7T9jwloAhWlhwkRI7c00XOUH+bPd3TnfjJWg9j/mIA5c8rzcmo8MiFi5JYgVf+AAbl8jvJDDof9l16ai11cmwPPA+Eb33COhw/PpUeWHxMiRm4pXZTdr58TMijni7KzJUfDfq8//MxnHM3akCG57w+zJcgDYdOmXHpk+TEhYuSaUlX/9u25UfXnk6DpW3t7w4f9Qf1hb2/u+8NsaSKPLD/mnWU0BebhEwO/p1b//s5qdpGGemrZ/5WAHHpk9a2GeWcZTUwOVf35xT99a2tz8ho8fbP/KwF77RWcn3MPBBMiRlOQI1V/8zB/vqNT99Mg9Yj9XzHwjEcvvlh+rgk8EDITIiIyXETuF5Gl7mvY/urL3c2nFotIT9zrjWIQpOpva8v985UtGU4HvvjF8rwm6A8bTwGW82c5E5kHPKCq+wIPuMdhHK2qh5To4+JcbzQ5pZ5a73qXY6j9ylccVbJ5awWQwXTAG1R/5SvO8W67mUdWRYKM6ao7jUfNcMOCdqpqRAKeBca478cAz4aUWw6MTHp9abKdDYvB1VeXb/TW0WE74/UhaPdD/654db5ZQV9n/0kVRIL/H5Gsa1YGOdzZcLSqrgZwX3cPKafAb0TkcRGZm+B6o4B8+9vleU3gDdlYGhzNskk9VLOlAMajVIWIiPxWRJYEpJNjfMyRqnoYcDzwBRGZkqAec0WkR0R61q5dG/dyI4eY909EPE+tIEFS5x7e/pMEnHZaeV6TGY9SFSKq+hFVPSgg/Qp4VUTGALiva0I+42X3dQ3wS2CyeyrS9e61C1V1kqpOGjVqVP1+oJEZBRjANZYG9PBN6qGaDd3djmD/5392jEYjRjSt8ShLddYdwGfd958FflVaQER2EZEh3ntgOrAk6vVGcQny1hKxzatCSVHqNrmHauPxPLI8Aa4KGzc273L+IENJIxIwAseraqn7OtzN3xO4232/D/CEm54G5le7vloyw3px6Opy7MNBdkkz6JaQkpE96GM9W3EKtvtiENZoOzuzrllFCDGsW9gTo+mxEBsR8cKhhO1Rm0CNYvc+ATkPbxKGhT0xCosZdCOSgpHd7n0C9tgjOL9JjUcmRIymx4zsMaljz2/G9Bh4xvTVq8vPNbHxyISI0fSYkT0mdZC6ZkyPSakxHZouvEkYJkSMpidoTZ2ncrZ9RwIIk7rr10eKIVOAcE+NpwjhTUIww7pRKMzQGxH/niODBjkupn4qGNrtHiegSY3pfsywbrQEZuiNiGdk374dghbgVjC02z1OwOjRwfkFMB6ZEDEKhRnZExBk2IBQqWDG9Bh4xvRXXik/VxDjkQkRo1CYkT0BESWvGdNjUmBjuh8TIkahMCN7AiJIXjOmJ6DAxnQ/Zlg3CosZgGNQZTX7+MGvsuL1XctO2b2sQAGM6X7MsG60HGYAjoFnaA9SbW3YwMrXO8rzsXsZiKf3CxugF8x4ZELEKCxmZE9AiJF9NAGGYexelhGk9/NTQOORCRGjsASp+sGM7BUpNaZzOp0s5xXGlBUtYH9YO0F2EI+CGo9MiBiFpcG7wxYDn+Tt5nTmcjUr6QQcK7rgqGgK2h/WTph+T6RQxnQ/Zlg3WgIzssfANbKPX7GIFYwvO905Yj3LXys3shvA2LHw8svl+QVoaGZYN1oaM7LHwDWyryTY4BFmZG9ZPEN6v37w6qvl5wuu98tMiIjIcBG5X0SWuq+7BZTZT0QW+9LfRORL7rlvishLvnMzGv4jjKYhzACsavaRUnY4F7kqrFLGYZJ3B35Duir09kJbW1PvmR6XLGci84AHVHVfnO1t55UWUNVnVfUQVT0EmAhsAH7pK/ID77yq3t2IShvNSZiRHcw+4qevc1G5EOngHRZwoUlejyBDem8v7LqrsxakoHYQP1kKkZOBG9z3NwCnVCl/DPBXVQ3xnTOMcCoZ2SHxxn6FI9y5SOlkOQs5m5ncbJLXw/SkmQqR0aq6GsB93b1K+dOAm0vyzhWRJ0Xkp0HqMA8RmSsiPSLSs3bt2tpqbTQt3no6CdbStNJzH0qocxHKciY4AsSjlSVviy0orESqQkREfisiSwLSyTE/ZwBwEvDvvuwfA38HHAKsBr4fdr2qLlTVSao6aVRQ2GujpTD7SDlV+8QwO0grSt4WXFBYiVSFiKp+RFUPCki/Al4VkTEA7uuaCh91PPAnVd3h+qCqr6pqr6puB64GJqf5W4ziYPaRvkTqE0dcEXyyFSVvCy4orET/DL/7DuCzwGXu668qlD2dElWWiIzx1GHAx4ElSSrx6quv8sYbbyS51IjI8OHDGR22KU8GeM93WLxBT0vTKv1AtT5xwQKYyeEw99rggp7khda4adUWFLYYmS02FJERwK3AOGAl8ClVfUNE9gSuUdUZbrkO4EVgH1Vd57v+JhxVlgLLgX/wCZVQShcbPvPMM7znPe+hra2tXj/N8NHb28tzzz3H/vvvn3VVAilYoNVERL4HlSL9QiEW1FWkxX9/7hYbqurrqnqMqu7rvr7h5r/sCRD3eIOqjvALEDd/lqq+T1Xfr6onRREgYZgASY+839tWtY/418dFtg23smeC2UFCsRXrDUBVOeGEE5gyZQq9vb078l966SW+9rWvBV7zpS99qU/ZUmbPns2mTZsCz82dO5d33nkHgHnz5rFq1aoaal9sWtE+Uro+LoiKfWIrSl6zg4RiQiQK/mFbgodk9erVDBkyhN///vd9RuY//vGPOeOMMwKvufLKK0NH8UuXLmX48OEMGjQo8PxJJ53EzTc7JqQzzjiDn/zkJ7Hq20q04vqRsP6wrS3iIutWkrzesx82AylwYMXIqGpLpYkTJ6qfP//5z1qRri7Vjg5VZ5zlpI4OJz8iM2fO1N13313nzJnTJ//II49UVdXNmzfrCSecoNOmTdNPfepTqqo6depU3bp1q1588cU6Z84cPeaYY3Zcf8UVV+h//Md/qKrqBRdcoHfeeaeuXr1aP/KRj+i2bdv07bff1lNOOWXH9xxzzDGR65oGVe9xThDp+zf7U2dnrL88l3R1Ob8j7DeK1PHDOjvT+RGNJOjZL+LvjAjQowF9apbeWfniS1+CxYvL8x99FDZv7pu3YQPMmQNXX903/5BD4Moryz7i29/+NgDXXHNNn/wtW7YAsHLlSkaOHMmdd96JBugXDjzwQK655hqmT5/OW2+9xdKlS5kyZQoAl156KccffzzDhg3ju9/9Lm1tbey66668/vrrZd9jVGbcuPABZ7M7IHkqrDCNDMRcHzdzppPCjCrepi0LFjTnDYPKKixoaTuIH1NnVaNUgFTLj4G4Bsp3v/vdvO9972PmzJn84Ac/KCt30EEHAbDnnnuybt26PoJm8ODBfPSjH2Xjxo0ccsghNdeplamkpYHmVm2l1h9WkjzNqtqqpsKClreD+LGZiEfADAKovBHFokU1fWV7ezsAmzdv5stf/jL9+vVj+vTpzCxpmOLzhlFV3vOe97B8+XImTpzI6tWreeihh5gwYQKLFi1i2rRprF+/nuHDh++4ZsCAATXVs1Wotn4Emm+AXc0rFXxrQZL8ngULKk9xmm3RTZQpW8FdeeNiM5FqBA1P6zSNnTJlCkuWLGHFihVMnTqVD37wg4waNYrdd68cRuyEE07gv/7rvwA477zz+P73v893v/tdLrnkEjZu3MiDDz7IjBmOl/RTTz3F5Mm2mD8qnhdrmKEdmmeAXc0rFXb2h4n7+GqeCdBc+xGbCis+QYaSIqfYhnXVnUZEkbpaWF988UX9p3/6p0TXzp49Wzdu3Bh47uyzz9b169erqurXvvY1XblyZeI61oNmMaz7KYJNtZLdO4F/SA6/sI5UcxQoindFDRBiWM+8U290SiREjJpo1nvcrP1KZvWOInnzeNOKMGJoAGFCxNRZhhFCM6q2GqLCCiOKagvyd9NMhVUTJkQMowpRvLbOOCNbtb/nUHTGGRn3h1EkL+TD1c28sOqCCRHDqELeB9hRZh/Q4P6wmuSFbA3umU7ZioUJEcOIQJwBdqNmJVFnH5BBf5hXyZubKVtxMCFiGDGIMsCG9PvGqLMPyLA/9CRvV1e2+kBPcIjArFk5m7I1PyZEIlBj/MU+XH/99Vx//fWRylaL5Ovn7rvv5tBDD+Xaa6/tk3/BBRfwyiuvxK3mDjZs2MBcL95HCZs3b2b27Nk7jmfNmuW4/BWYqANsSKdvjDP7gJz0h3FmJbNmOZ19vW5aqbSt1j5NhRUbEyJVKA2b3ajZ9/bt2ytG8i3l9ttv5+abb2bOnDk78v72t7+xZs0a9thjj8T1+NnPfsaJJ54YeG7gwIEMHz6cpUuXAnD44YfzwAMPJP6uZiHqANtjxQo46ywYOTLZQCTuQBqcenV15ag/jKoP9Dr5WgVKXGkLpsJKSpDfbyMS8CngaWA7MKlCueOAZ4FlwDxf/nDgfmCp+7pblO8NWydy3nmqU6eWp4EDg93GBw4sL3veecH+1Zs3b9aTTjpJjz32WD3ttNP0uuuuU1XVSy65RKdOnapHH320vvDCC/rCCy/otGnT9JOf/KRee+21OyL5TpkyRXt7e1VV9dRTT9XVq1f3+fyHH35YR48erZMnT9aHH354R/5tt92m3/ve91RVdf369frJT35Sp0yZomeeeaaqqt500016+OGH6xFHHKGLFy8OjCZ80kkn6Wuvvabbtm3T6dOn66pVq/Tee+/V89wfe9ddd+n3v/99VVV9+umn9Ytf/GLZ72/WdSJRiLImIyxiLqiOGOEkkeD3/rJRU96WYfQh6lqS0tTeXvlGFfqm5QNyuE5kCfAJ4PdhBUSkDfgRcDxwAHC6iBzgnp4HPKCq+wIPuMd1px7xF2+//XYmT57Mvffey5AhQwAnHMlLL73EokWL+NGPfsR3vvMdANasWcPPf/7zPmqiI444gj/84Q9s2LCBd955p2xmcdRRR3HcccfR3d3NUUcdtSN/6dKljB8/HoCFCxcyffp0HnroIa699lp6e3u56qqrePjhh+nu7mb+/Pk7ogk/+OCD/PznPwecvVBGjBhBW1sbV1xxBZ///Of5zne+wwJ3xLbPPvvwl7/8pex9qxB3VuLhDbhff91JqsHv/WWrkbvZRxBx9IF+tm6tfKMKfdPyTWYBGFX1GegbXDCAycAyVX3eLXsLcDLwZ/d1mlvuBmARELxNYATSjL/4/PPPc+ihhwIwceJEwNnb3QuYCDBmzBgADj744DIV1qmnnsqNN97ImjVrOP7446N9KfSxTzz33HN84QtfAKBfv3688sordHZ20t7ezvjx41m3bl2faMITJ07k/PPP7/N5Bx54IBs3buTYY49ll1122fEd/v+wyv9ZWKIEb0yTmoIoNhovjHyUYIdpIOIImaa6afkl7zaRscCLvuNVbh7AaHX3VXdfQ6MWishcEekRkZ61a9fGqkA94i9OmDCBJ554AoD/+Z//AWC//fZj+vTpLFq0iEWLFnHjjTcCTgdfysSJE1m8eDG/+MUv+MQnPhH5e71ov973Pfroo4Bjbxk1ahTLly9n69atLF++nKFDh+6IJtzd3c29997Lq6++ypgxY3bsTXLPPfew//7789vf/pbXXnsNgBdeeIH99tsPcISl974VSTorqYWmHkiXzkoaMQDp7ISbbnKESFPetPyRqhARkd+KyJKAdHLUjwjIi+3+o6oLVXWSqk4aNWpUrGv97TzS1qEBnHLKKTzyyCMce+yxvPXWW4Az49hjjz2YNm0aRx99NNddd13Fzzj88MNZtWoVY8eOrVjOz4c//GEWuxttnX322dxzzz1MnTqVz33uc7S1tXHuuefyoQ99iE9/+tNceumlgdGEvYjBb7/9NpdffjmXXXYZl112GRdccAEAv/vd7/jYxz4GwAMPPLDjfStT2mZGjIB6RuP3+tpceF7Viid5VZ3OPS2B0tTSNucEGUoamXDUUIGGdeCDwH2+468DX3ffPwuMcd+PAZ6N8n2tFoDx/PPPLzPEx2H9+vV69tlnB57bvHnzDiO9qrMNsOcA4Kfo9zgKfgN8XJuv/5qWsf/6I2ePGKE6YIDdtIwhr1F8qwiR/sDzwARgAPAEcKB77ru43lo4RvXLo3xfEYTIlVdeqVOnTt2RrrzyyqyrVJFmvMdpUto/VnM6sj5Q7ablgDAhIs65xiMiHwd+CIwC3gIWq+qxIrIncI2qznDLzQCuBNqAn6rqAjd/BHArMA5YCXxKVd+o9r2TJk3Snp6eHcfPPPMM733ve1vWIJw2qspf/vIX9t9//6yrYhhGDYjI46o6qTQ/S++sXwK/DMh/GZjhO74buDug3OvAMbXWo729nU2bNjF48OBaP8oIYNOmTTu2ATYMo3i0/B7rI0eO3OHBZKSD575sGEbxaHkhMmzYMIYNG5Z1NQzDMJqSvK8TMQzDMHKMCRHDMAwjMZl5Z2WFiKwFkgamGAm8Vsfq1AurVzysXvGwesWjqPXqVNWy1dotJ0RqQUR6glzcssbqFQ+rVzysXvFotXqZOsswDMNIjAkRwzAMIzEmROKxMOsKhGD1iofVKx5Wr3i0VL3MJmIYhmEkxmYihmEYRmJMiBiGYRiJMSHiIiLHicizIrJMRMr2axeHq9zzT4rIYVGvTbleM936PCkij4jIwb5zy0XkKRFZLCI9pdemXK9pIrLO/e7FIvKNqNemXK+v+uq0RER6RWS4ey6V+yUiPxWRNSKyJOR8Vm2rWr2yalvV6pVV26pWr4a3Lfez9xaRB0XkGRF5WkTOCyiTXhsLig/fagknzPxfgX3YuW/JASVlZgD34Oy2+AHgv6Nem3K9jgB2c98f79XLPV4OjMzofk0D7kxybZr1Kil/IvC7BtyvKcBhwJKQ8w1vWxHr1fC2FbFeDW9bUeqVRdtyP3sMcJj7fgjwXCP7L5uJOEwGlqnq86q6BbgFKN3C92TgRnV4FBgmImMiXptavVT1EVV90z18FNirTt9dU71Surben306cHOdvjsUVf09UGmvmyzaVtV6ZdS2otyvMDK9XyU0pG0BqOpqVf2T+/5t4BmgdB/t1NqYCRGHscCLvuNVlP8JYWWiXJtmvfzMwRlteCjwGxF5XETm1qlOcer1QRF5QkTuEZEDY16bZr0QkQ7gOOAXvuy07lc1smhbcWlU24pKo9tWZLJsWyIyHjgU+O+SU6m1sZYPBe8StK1hqe9zWJko1yYl8meLyNE4D/pRvuwjVfVlEdkduF9E/uKOphpRrz/hxNpZL87ulLcD+0a8Ns16eZwI/EH77oaZ1v2qRhZtKzINbltRyKJtxSGTtiUiu+IIri+p6t9KTwdcUpc2ZjMRh1XA3r7jvYCXI5aJcm2a9UJE3g9cA5yszo6PwI5dIlHVNTi7SE5uVL1U9W+qut59fzfQLiIjo1ybZr18nEaJuiHF+1WNLNpWJDJoW1XJqG3FoeFtS0TacQRIt6reFlAkvTaWhqGn2RLOjOx5YAI7jUsHlpQ5gb6GqT9GvTbleo0DlgFHlOTvAgzxvX8EOK6B9dqDnYtZJwMr3XuX6f1yyw3F0W3v0oj75X7meMINxQ1vWxHr1fC2FbFeDW9bUeqVYdsS4EbgygplUmtjps4CVHWbiJwL3IfjrfBTVX1aRM5xz/8EZ5/3GTgP1QbgrErXNrBe3wBGAP8mIgDb1InUORr4pZvXH/iZqt7bwHqdCnxeRLYBG4HT1Gm1Wd8vgI8Dv1HVd3yXp3a/RORmHI+ikSKyCrgYaPfVqeFtK2K9Gt62Itar4W0rYr2gwW3L5UhgFvCUiCx28y7EGQSk3sYs7IlhGIaRGLOJGIZhGIkxIWIYhmEkxoSIYRiGkRgTIoZhGEZiTIgYhmEYiTEhYhiGYSTGhIhhZIyIjBeRjT4f/7Byg91Q4lvcFdqGkTkmRAwjH/xVVQ+pVEBVN7plsgjlYRiBmBAxjBQRkb93NwEaJCK7uJsGHVSh/LtFZK27idFiEXlDRP4qIu9qZL0NIyoW9sQwUkRVHxORO4BvA4OBLlUN3BnPLb9MRP4TuEJVHxaRRcAXtTwqq2HkAhMihpE+3wIeAzYB/xih/IGAJ2jeCzybUr0Mo2ZMnWUY6TMc2BVn69JBlQqKyGBgkKq+KSJ7A6+rs+OcYeQSEyKGkT4LgYuAbuCfq5Q9AGd7U4D9fe8NI5eYEDGMFBGRz+CEUP8ZcBnw9yLy4QqX+FVZG4HDROS9KVfTMBJjoeANI2PcfbHvVNVQr62S8suBSar6Wpr1Mowo2EzEMLKnFxgadbEhzkZI2xtQL8Oois1EDMMwjMTYTMQwDMNIjAkRwzAMIzEmRAzDMIzEmBAxDMMwEmNCxDAMw0iMCRHDMAwjMSZEDMMwjMSYEDEMwzAS8/8BG/hQunj9IoAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# note: numpy and matplotlib are already included\n",
    "# create list of evenly spaced numbers\n",
    "x = np.linspace(0.0, 2.0*np.pi, 100)\n",
    "\n",
    "# save the spacing in the variable 'h'; since the numbers in 'x' are evenly spaced it doesn't \n",
    "#matter which difference we use\n",
    "h = x[10] - x[9]\n",
    "\n",
    "# create array with the sine of that numbers\n",
    "f = np.sin(x)\n",
    "\n",
    "# calculate the derivatives of the points in 'f' using the spacing 'h' and store the results in the variable 'deriv_bf\n",
    "deriv_f = np.gradient(f, h)\n",
    "\n",
    "# plot the array 'b' and its derivative 'deriv_b'\n",
    "plt.plot(x/np.pi, f,       marker = 'o', color = 'red',  label = 'f (sin(x))')\n",
    "plt.plot(x/np.pi, deriv_f, marker = 'o', color = 'blue', label = 'deriv_f (cos(x))')\n",
    "plt.xlabel(r'x [$\\pi$]')\n",
    "plt.ylabel('f / deriv_f')\n",
    "plt.legend(loc = 'lower left', fontsize = 8)\n",
    "plt.show()\n",
    "\n",
    "# Note: As expected, 'f' shows the sine function and its derivative 'deriv_f' is the cosine!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "\n",
    "Compute analytically derivative of $\\sin{x}$ and evaluate it at the given `x`-array. Convince yourself that the difference of this result and the numeric derivative `deriv_f` is approximately of  order $\\mathcal{O}(h^2)$. In addition, evaluate the derivative using the  **scipy.misc.derivative**, and compute the difference of this result and the previous two methods of computation. Plot the analytic, and the two numeric derivatives in a scale that makes visible their differences. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3 (optional)\n",
    "\n",
    "To get more practice at the computation of numerical derivatives with Python:\n",
    "1. Calculate the derivative of the function $f(x) = x^4 - 4x^2 + 3x - 1$ at the point $x = 1.5$ both analytically (by hand + using your calculator) and numerically (using $h = dx = 0.1$ and $h = dx = 0.01$). Compare the results. What do you expect to change when using different values for $dx$?\n",
    "2. Calculate the derivative of the function $g(x) = 1/x$ at the point $x = 0.1$ both analytically (by hand + using your calculator) and numerically (using $h = dx = 0.1$). What happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution to task 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution to task 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical Integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a 1-dimensional integral of the form\n",
    "$$\\int_{a}^{b} f(x) dx,$$\n",
    "exploiting the  linearity of the integral, this integral  can be rewrite it as:\n",
    "$$\\int_{a}^{b} f(x) dx = \\sum_{i=0}^{N} \\int_{x_i}^{x_{i+1}} f(x) dx$$\n",
    "where $x_i = a + i h$ and $h = (b-a)/N$. i.e. the  $x_i$ points are  equally spaced.\n",
    "\n",
    "If $x_i$ and $x_{i+1}$ are close to each other, the integral $\\int_{x_i}^{x_{i+1}} f(x) dx$ can be approximated with different forms. \n",
    "\n",
    "### Rectangle rule\n",
    "\n",
    "The rectangle rule is the simplest form of the integral. Each integral in the previous sum is approximated by the area of a rectangle of base $(x_i - x_{i+1})$ and height $ f(\\frac{x_i + x_{i+1}}{2})$, i.e. the curve passes through the central points in one of the rectangle's edges. \n",
    "\n",
    "$$\\int_{x_i}^{x_{i+1}} f(x) dx = (x_i - x_{i+1}) f\\left(\\frac{x_i + x_{i+1}}{2}\\right) + \\mathcal{O}(1/N^2)$$\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/2/26/Integration_rectangle.svg\" alt=\"drawing\" width=\"60%\"/>\n",
    "\n",
    "### Trapezoidal rule\n",
    "\n",
    "A little bit more complex is the trapezoidal rule. Each integral in the sum is approximated by the area of a trapezoid. Notice that the trapezoid has two different heights, and to compute its area we have have to evaluate the function at two different $x_i$ values. The area of the trapezoid is then computed by\n",
    "\n",
    "$\\int_{x_i}^{x_{i+1}} f(x) dx = (x_i - x_{i+1}) \\frac{f(x_i) + f(x_{i+1})}{2} + \\mathcal{O}(1/N^3)$\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/0/03/Integration_trapezoid.svg\" alt=\"drawing\" width=\"60%\"/>\n",
    "\n",
    "### Simpson's rule\n",
    "\n",
    "The most complicated rule out of the three is the Simpson's rule. The function of each interval is approximated  by a polynomial of order 2 and integrating that polynomial instead of the original function. It's error propagation is of order $\\mathcal{O}(1/N^5)$\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/1/13/Integration_simpson.svg\" alt=\"drawing\" width=\"60%\"/>\n",
    "\n",
    "<style>\n",
    "img {\n",
    "  background-color: white;\n",
    "}\n",
    "</style>\n",
    "\n",
    "As you can see, the error for these methods evolve differently for different values of $N$. It's important to keep in mind, that different methods also require different numbers of function calls. So if f(x) is a very expensive function without a lot of fluctuations, it's better to use the rectangle rule, than the trapezoid rule.\n",
    "\n",
    "There is also a way to vary the size of the intervals depending on the function, but this is a more advanced method. We will just mention it here, for you to know that this exists. \n",
    "\n",
    "Also keep in mind, that the methods presented here are just the common three out of many which are optimized for different problems. \n",
    "\n",
    "[Source of images](https://en.wikipedia.org/wiki/Numerical_integration)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n",
    "\n",
    "Using `numpy`, program the rectangle and trapezoid rules in the next code block to integrate the function $f(x) = e^{-x}$ from $x=0$ to $x=1$. Compare it with the analytic result varying the the number of intervals and see what the difference for the mentioned methods of computation. Plot all your findings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your solutions here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Simpson's rule is quite complicated to implement and other integration methods can be even more complicated. For that reason, clever people have already done this for us and bundled them in a library. The`scipy` library contains many methods for numerical integration. We will now take a look at some of them, but a comprehensive list is [here](https://docs.scipy.org/doc/scipy/reference/integrate.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `scipy.integrate.quad`\n",
    "\n",
    "Scipy's `quad` function is the main workhorse that you should use for numerical one-dimensional integrals. \n",
    "It is based on the FORTRAN library **QUADPACK**. Scipy has generally a really good [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) that we encourage you read.\n",
    "\n",
    "Given a 1-dimensional analytic integrand (or interpolated) function  `func`, the 1-dimensional quad integration has the following syntax:\n",
    "```py\n",
    "result, error = integrate.quad(func, x_min, x_max, args=())\n",
    "```\n",
    "where `args` are possible parameters entering in your function. \n",
    "Run the commands `import scipy.integrate as integrate` and  `integrate.quad?` for additional documentation.\n",
    "\n",
    "Let us see how to use `quad` integration with an specific example.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999546000702375 2.8326146575791917e-14\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate as integrate\n",
    "\n",
    "def f(x):\n",
    "    return np.exp(-x)\n",
    "\n",
    "# calculate the integral of 'f' from 0 to 10\n",
    "result, error = integrate.quad(f, 0, 10)\n",
    "\n",
    "print(result, error)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first number gives you the result, the second one the numerical error of the integration. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While `quad` is useful for known integrand functions, other methods implemented in the `scipy` library are useful when having integrand functions that are `np.arrays`. One of such a methods is the already mentioned **trapezoid** rule.\n",
    "The syntax for implementation of this integration method is the following:\n",
    "\n",
    "```py\n",
    "scipy.integrate.trapezoid(y, x=None, dx=1.0, axis=-1)\n",
    "```\n",
    "or \n",
    "```py\n",
    "scipy.integrate.trapz(y, x=None, dx=1.0, axis=-1)\n",
    "```\n",
    "\n",
    "Let us see an specific example of this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.000834863090761"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(0, 100, 1000)\n",
    "y = f(x)\n",
    "\n",
    "integrate.trapz(y, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the result is fairly similar to the quad method, but you don't have an error estimate.\n",
    "\n",
    "The last method we will be presenting here, is **Simpson's method**, which works similar as the trapezoid method. The syntax for the implementation is:\n",
    "\n",
    "```py\n",
    "scipy.integrate.simpson(y, x=None, dx=1.0, axis=-1)\n",
    "```\n",
    "or \n",
    "```py\n",
    "scipy.integrate.simps(y, x=None, dx=1.0, axis=-1)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000402922553666"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linspace(0, 100, 1000)\n",
    "y = f(x)\n",
    "integrate.simps(y, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Depending on the problem we are trying to solve, we can use either of the different integration methods in `scipy` library.  While `quad` works just with functions and the step size is estimated on the fly, the implementation of the `trapezoid` and `simpson` methods work on arrays. \n",
    "\n",
    "Keep in mind that one can convert an array into a function and a function into an array, for that, you we have prepared an  optional section on interpolating functions, which we will extend in Lecture 7. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 5\n",
    "\n",
    "In this exercise we will compute numerically the integral of the function \n",
    "\n",
    "$$y(x,a) = \\int_0^{x} e^{-a\\phi^2} d\\phi,$$\n",
    "using the `quad` integration method described above. Here $a$ is a positive real number. For $x\\to \\infty$, this integral has an analytic result \n",
    "\n",
    "$$\\int_0^{\\infty} e^{-a \\phi^2} d\\phi = \\sqrt{\\frac{\\pi}{4a}}$$\n",
    "\n",
    "which we will use to compare to our numeric integration. \n",
    "\n",
    "1. Write a function, which represents the integrand and a function `y(x,a)` which calculates the integral. \n",
    "Recall the `quad` function allows for integration when parameters are included inside the integrand, as well as infinite integration interval.\n",
    "\n",
    "2. Check that your function is well defined by calculating `y(np.inf,10)` and compared to the analytic result using `a = 10`.\n",
    "\n",
    "3. Use the interval `x = np.linspace(0, 100, 100)` to evaluate `y(x,10)`\n",
    "\n",
    "5. Plot your result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your results here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two-dimensional integrals (`scipy.integrate.dblquad`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Higher-dimensional integration is also possible to perform. In particular, for two dimensions we have `scipy.integrate.dblquad`. \n",
    "\n",
    "For instance, the area of a circle in Cartesian coordinates can be computed from the double integral.\n",
    "\n",
    "$$A = \\int_{x=-1}^{1} \\int_{y=-\\sqrt{1-x^2}}^{\\sqrt{1-x^2}} dy dx$$\n",
    "\n",
    "To understand this further, lets rewrite the integral, to make the variables match with the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.dblquad.html):\n",
    "\n",
    "$$\\int_{x=a}^{b} \\int_{y=\\text{gfun}(x)}^{\\text{hfun}(x)} \\text{fun}(y,x) dy dx$$\n",
    "\n",
    "where:\n",
    "\n",
    "$$a = -1, \\quad\n",
    "b = 1, \\quad\n",
    "\\text{gfun}(x) = -\\sqrt{1-x^2}, \\quad\n",
    "\\text{hfun}(x) = \\sqrt{1-x^2}, \\quad \n",
    "\\text{fun}(y,x) = 1$$\n",
    "\n",
    "Let's program this integral using our `dblquad` documentation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3.141592653589797, 2.000471344132393e-09)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = -1\n",
    "b = 1\n",
    "\n",
    "def gfun(x):\n",
    "    return -np.sqrt(1-x**2)\n",
    "\n",
    "def hfun(x):\n",
    "    return np.sqrt(1-x**2)\n",
    "\n",
    "def fun(y, x): # the order of the arguments is important\n",
    "    return 1\n",
    "\n",
    "\n",
    "integrate.dblquad(fun, a, b, gfun, hfun)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And it is as easy as that! "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4\n",
    "\n",
    "Using the same approach, calculate the volume of a sphere using the function `scipy.integrate.tplquad`, which is documented [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.tplquad.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation (optional)\n",
    "\n",
    "Given two arrays with data points (y, x), we have precise values where the data points are recorded, but sometimes one wants an estimate of the values in between these data points. A process called interpolation helps us to do this. We will roughly outline interpolation and discuss it in detail in  Lecture 7. \n",
    "\n",
    "Let's say that we have the following data points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1ce3844a520>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAT0UlEQVR4nO3df6hc553f8fdnZdW5JDFy8LWRr5TKBK2onbBSPagGQ0nTdKVNl5USCCilsf8wKBgHErq4tfafTf4ICXiTtIbaoGyM7TaNEcSxRWqvVusmhIBjZWQrkWVHjVh7Y/3AutsgYoNwbeXbP+Yoncg30r3SvTNX93m/YJgz33nOOc8Z44+OnnlGT6oKSVIb/mDcHZAkjY6hL0kNMfQlqSGGviQ1xNCXpIZcMe4OXMg111xTa9asGXc3JOmysn///n+sqslz64s+9NesWUO/3x93NyTpspLkH2aqO7wjSQ0x9CWpIYa+JDXE0Jekhhj6ktSQRT97R5Ja8vjzx7h3z2GOnzrN9SsmuHvTOrZumJq34xv6krRIPP78MXY8dpDTb50B4Nip0+x47CDAvAW/wzuStEjcu+fwbwP/rNNvneHePYfn7RyGviQtEsdPnZ5T/WIY+pK0SFy/YmJO9Yth6EvSInH3pnVMLF/2O7WJ5cu4e9O6eTuHX+RK0iJx9staZ+9IUiO2bpia15A/l8M7ktQQQ1+SGmLoS1JDLhj6Sd6VZF+SnyY5lOSLXf0LSY4lOdA9Pja0z44kR5IcTrJpqH5zkoPde/clycJcliRpJrP5IvdN4CNV9UaS5cCPkjzVvff1qvqr4cZJbgS2ATcB1wN/l+QPq+oM8ACwHfgx8CSwGXgKSdJIXPBOvwbe6F4u7x51nl22AI9W1ZtV9TJwBNiYZCVwVVU9U1UFPAJsvaTeS5LmZFZj+kmWJTkAnAT2VtWz3VufTfKzJA8mubqrTQGvDu1+tKtNddvn1mc63/Yk/ST96enp2V+NJOm8ZhX6VXWmqtYDqxjctX+QwVDNB4D1wAngq13zmcbp6zz1mc63s6p6VdWbnHzHYu6SpIs0p9k7VXUK+AGwuape6/4w+A3wDWBj1+wosHpot1XA8a6+aoa6JGlEZjN7ZzLJim57Avgo8PNujP6sjwMvdNu7gW1JrkxyA7AW2FdVJ4DXk9zSzdq5DXhi/i5FknQhs5m9sxJ4OMkyBn9I7Kqq7yX5b0nWMxiieQX4DEBVHUqyC3gReBu4q5u5A3An8BAwwWDWjjN3JGmEMphIs3j1er3q9/vj7oYkXVaS7K+q3rl1f5ErSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhs1kY/V1J9iX5aZJDSb7Y1d+XZG+SX3TPVw/tsyPJkSSHk2waqt+c5GD33n3dAumSpBGZzZ3+m8BHquqPgPXA5iS3APcAT1fVWuDp7jVJbgS2ATcBm4H7u0XVAR4AtgNru8fm+bsUSdKFXDD0a+CN7uXy7lHAFuDhrv4wsLXb3gI8WlVvVtXLwBFgY5KVwFVV9UwNVmN/ZGgfSdIIzGpMP8myJAeAk8DeqnoWuK6qTgB0z9d2zaeAV4d2P9rVprrtc+sznW97kn6S/vT09BwuR5J0PrMK/ao6U1XrgVUM7to/eJ7mM43T13nqM51vZ1X1qqo3OTk5my5KkmZhTrN3quoU8AMGY/GvdUM2dM8nu2ZHgdVDu60Cjnf1VTPUJUkjMpvZO5NJVnTbE8BHgZ8Du4Hbu2a3A09027uBbUmuTHIDgy9s93VDQK8nuaWbtXPb0D6SpBG4YhZtVgIPdzNw/gDYVVXfS/IMsCvJHcAvgU8CVNWhJLuAF4G3gbuq6kx3rDuBh4AJ4KnuIUkakQwm0ixevV6v+v3+uLshSZeVJPurqndu3V/kSlJDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqyGwWRl+d5PtJXkpyKMnnuvoXkhxLcqB7fGxonx1JjiQ5nGTTUP3mJAe79+7rFkiXJI3IbBZGfxv486p6Lsl7gf1J9nbvfb2q/mq4cZIbgW3ATcD1wN8l+cNucfQHgO3Aj4Engc24OLokjcwF7/Sr6kRVPddtvw68BEydZ5ctwKNV9WZVvQwcATYmWQlcVVXP1GA19keArZd6AZKk2ZvTmH6SNcAG4Nmu9NkkP0vyYJKru9oU8OrQbke72lS3fW59pvNsT9JP0p+enp5LFyVJ5zHr0E/yHuA7wOer6tcMhmo+AKwHTgBfPdt0ht3rPPV3Fqt2VlWvqnqTk5Oz7aIk6QJmFfpJljMI/G9V1WMAVfVaVZ2pqt8A3wA2ds2PAquHdl8FHO/qq2aoS5JGZDazdwJ8E3ipqr42VF851OzjwAvd9m5gW5Irk9wArAX2VdUJ4PUkt3THvA14Yp6uQ5I0C7OZvXMr8GngYJIDXe0vgE8lWc9giOYV4DMAVXUoyS7gRQYzf+7qZu4A3Ak8BEwwmLXjzB1JGqEMJtIsXr1er/r9/ri7IUmXlST7q6p3bt1f5EpSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDZrNG7uok30/yUpJDST7X1d+XZG+SX3TPVw/tsyPJkSSHk2waqt+c5GD33n3dWrmSpBGZzZ3+28CfV9U/A24B7kpyI3AP8HRVrQWe7l7TvbcNuAnYDNyfZFl3rAeA7QwWS1/bvS9JGpELhn5Vnaiq57rt14GXgClgC/Bw1+xhYGu3vQV4tKrerKqXgSPAxiQrgauq6pkaLMz7yNA+kqQRmNOYfpI1wAbgWeC6qjoBgz8YgGu7ZlPAq0O7He1qU932ufWZzrM9ST9Jf3p6ei5dlCSdx6xDP8l7gO8An6+qX5+v6Qy1Ok/9ncWqnVXVq6re5OTkbLsoSbqAWYV+kuUMAv9bVfVYV36tG7Khez7Z1Y8Cq4d2XwUc7+qrZqhLkkZkNrN3AnwTeKmqvjb01m7g9m77duCJofq2JFcmuYHBF7b7uiGg15Pc0h3ztqF9JEkjcMUs2twKfBo4mORAV/sL4CvAriR3AL8EPglQVYeS7AJeZDDz566qOtPtdyfwEDABPNU9JEkjksFEmsWr1+tVv98fdzck6bKSZH9V9c6t+4tcSWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNmc3C6A8mOZnkhaHaF5IcS3Kge3xs6L0dSY4kOZxk01D95iQHu/fu6xZHlySN0Gzu9B8CNs9Q/3pVre8eTwIkuRHYBtzU7XN/kmVd+weA7cDa7jHTMSVJC+iCoV9VPwR+NcvjbQEerao3q+pl4AiwMclK4KqqeqYGK7E/Amy9yD5Lki7SpYzpfzbJz7rhn6u72hTw6lCbo11tqts+tz6jJNuT9JP0p6enL6GLkqRhFxv6DwAfANYDJ4CvdvWZxunrPPUZVdXOqupVVW9ycvIiuyhJOtdFhX5VvVZVZ6rqN8A3gI3dW0eB1UNNVwHHu/qqGeqSpBG6qNDvxujP+jhwdmbPbmBbkiuT3MDgC9t9VXUCeD3JLd2snduAJy6h35Kki3DFhRok+TbwYeCaJEeBvwQ+nGQ9gyGaV4DPAFTVoSS7gBeBt4G7qupMd6g7GcwEmgCe6h6SpBHKYDLN4tXr9arf74+7G5J0WUmyv6p659b9Ra4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ15IKhn+TBJCeTvDBUe1+SvUl+0T1fPfTejiRHkhxOsmmofnOSg91793Vr5UqSRmg2d/oPAZvPqd0DPF1Va4Gnu9ckuRHYBtzU7XN/kmXdPg8A2xkslr52hmNKkhbYBUO/qn4I/Oqc8hbg4W77YWDrUP3Rqnqzql4GjgAbk6wErqqqZ2qwKO8jQ/tIkkbkYsf0r6uqEwDd87VdfQp4dajd0a421W2fW5ckjdB8f5E70zh9nac+80GS7Un6SfrT09Pz1jlJat3Fhv5r3ZAN3fPJrn4UWD3UbhVwvKuvmqE+o6raWVW9qupNTk5eZBclSee62NDfDdzebd8OPDFU35bkyiQ3MPjCdl83BPR6klu6WTu3De0jSRqRKy7UIMm3gQ8D1yQ5Cvwl8BVgV5I7gF8CnwSoqkNJdgEvAm8Dd1XVme5QdzKYCTQBPNU9JEkjlMFkmsWr1+tVv98fdzck6bKSZH9V9c6t+4tcSWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhlzwn1aWpEvx+PPHuHfPYY6fOs31Kya4e9M6tm5wtdRxMfQlLZjHnz/GjscOcvqtwbIax06dZsdjBwEM/jFxeEfSgrl3z+HfBv5Zp986w717Do+pRzL0JS2Y46dOz6muhWfoS1ow16+YmFNdC++SQj/JK0kOJjmQpN/V3pdkb5JfdM9XD7XfkeRIksNJNl1q5yUtbndvWsfE8mW/U5tYvoy7N60bU480H3f6/6qq1g+txXgP8HRVrQWe7l6T5EZgG3ATsBm4P8mymQ4oaWnYumGKL3/iQ0ytmCDA1IoJvvyJD/kl7hgtxOydLcCHu+2HgR8A/6mrP1pVbwIvJzkCbASeWYA+SFoktm6YMuQXkUu90y/gb5PsT7K9q11XVScAuudru/oU8OrQvke72jsk2Z6kn6Q/PT19iV2UJJ11qXf6t1bV8STXAnuT/Pw8bTNDrWZqWFU7gZ0AvV5vxjaSpLm7pDv9qjrePZ8EvstguOa1JCsBuueTXfOjwOqh3VcBxy/l/JKkubno0E/y7iTvPbsN/DHwArAbuL1rdjvwRLe9G9iW5MokNwBrgX0Xe35J0txdyvDOdcB3k5w9zv+oqr9J8hNgV5I7gF8CnwSoqkNJdgEvAm8Dd1XVmZkPLUlaCBcd+lX198AfzVD/P8C//j37fAn40sWeU5J0afxFriQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqyEIsjC4taY8/f4x79xzm+KnTXL9igrs3rXPhb102DH1pDh5//hg7HjvI6bcG6/8cO3WaHY8dBDD4dVlweEeag3v3HP5t4J91+q0z3Lvn8Jh6JM2NoS/NwfFTp+dUlxabkQ/vJNkM/BdgGfDXVfWV+T6HY65z4+c1e9evmODYDAF//YqJMfRGmruR3uknWQb8V+BPgBuBTyW5cT7PcXbM9dip0xT/f8z18eePzedplgw/r7m5e9M6JpYv+53axPJl3L1p3Zh6JM3NqId3NgJHqurvq+r/Ao8CW+bzBI65zo2f19xs3TDFlz/xIaZWTBBgasUEX/7Eh/ybkS4box7emQJeHXp9FPgX5zZKsh3YDvD+979/TidwzHVu/LzmbuuGKUNel61R3+lnhlq9o1C1s6p6VdWbnJyc0wl+39iqY64z8/OS2jLq0D8KrB56vQo4Pp8ncMx1bvy8pLaMenjnJ8DaJDcAx4BtwL+bzxOc/Wu3s1Fmx89Lakuq3jG6srAnTD4G/GcGUzYfrKovna99r9erfr8/iq5J0pKRZH9V9c6tj3yeflU9CTw56vNKkvxFriQ1xdCXpIYY+pLUEENfkhoy8tk7c5VkGviHi9z9GuAf57E7lwOvuQ2tXXNr1wuXfs3/tKre8evWRR/6lyJJf6YpS0uZ19yG1q65teuFhbtmh3ckqSGGviQ1ZKmH/s5xd2AMvOY2tHbNrV0vLNA1L+kxfUnS71rqd/qSpCGGviQ1ZEmGfpLNSQ4nOZLknnH3ZxSSPJjkZJIXxt2XUUiyOsn3k7yU5FCSz427TwstybuS7Evy0+6avzjuPo1KkmVJnk/yvXH3ZRSSvJLkYJIDSeb1nxlecmP63eLr/xv4NwwWbfkJ8KmqenGsHVtgSf4l8AbwSFV9cNz9WWhJVgIrq+q5JO8F9gNbl/J/5yQB3l1VbyRZDvwI+FxV/XjMXVtwSf4D0AOuqqo/HXd/FlqSV4BeVc37D9KW4p3+gi++vhhV1Q+BX427H6NSVSeq6rlu+3XgJQZrMC9ZNfBG93J591had20zSLIK+LfAX4+7L0vBUgz9mRZfX9Jh0Loka4ANwLNj7sqC64Y5DgAngb1VteSvmcGiS/8R+M2Y+zFKBfxtkv1Jts/ngZdi6M9q8XUtDUneA3wH+HxV/Xrc/VloVXWmqtYzWF96Y5IlPZSX5E+Bk1W1f9x9GbFbq+qfA38C3NUN386LpRj6C774uhaHblz7O8C3quqxcfdnlKrqFPADYPN4e7LgbgX+rBvjfhT4SJL/Pt4uLbyqOt49nwS+y2DYel4sxdD/7eLrSf4Jg8XXd4+5T5pn3Zea3wReqqqvjbs/o5BkMsmKbnsC+Cjw87F2aoFV1Y6qWlVVaxj8v/y/qurfj7lbCyrJu7vJCSR5N/DHwLzNyltyoV9VbwOfBfYw+HJvV1UdGm+vFl6SbwPPAOuSHE1yx7j7tMBuBT7N4M7vQPf42Lg7tcBWAt9P8jMGNzd7q6qJKYyNuQ74UZKfAvuA/1lVfzNfB19yUzYlSb/fkrvTlyT9foa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1Jasj/A9/dFnwtseHkAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(0.0, 5, 5)\n",
    "\n",
    "y = x**2*np.exp(x)\n",
    "\n",
    "plt.scatter(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want to estimate y(3), but as you can see, a value of y at x=3 doesn't exist. The simplest estimate is a linear interpolation, which means drawing a straight line between the points. This can be done with the `scipy.interpolate.interp1d` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(284.8654386)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import scipy.interpolate as interpolate\n",
    "\n",
    "y_linear = interpolate.interp1d(x, y)\n",
    "\n",
    "y_linear(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[scipy.interpolate.interp1d](https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpolate.interp1d) generates a function in the domain $(min(x), max(x))$ It easily can be understood by plotting it with a high point density, to see that it truly continuous."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1ce39499130>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX1ElEQVR4nO3db4xc1X3G8e+Ds4UNCVoQBtm7JrZSxyoG1ZSR68pSlZKk69IodqjSmKqBF0gbIaMSNXJr501JJQtLJNAiFSQnIEybxLUEMRZ/QhxMFBEZnDU2GGNc3ELCrld4o8SJUS0Xm19fzJkwrOfOzuzO//t8pNHcOXPv3DMYHh/OPfO7igjMzCwfzmt3B8zMrHUc+mZmOeLQNzPLEYe+mVmOOPTNzHLkQ+3uwHQuvfTSWLhwYbu7YWbWVfbt2/fLiJg7tb3jQ3/hwoWMjo62uxtmZl1F0s8rtXt6x8wsRxz6ZmY54tA3M8sRh76ZWY449M3McqTjV++YmeXJjv3j3PX0EY6dOMX8gX7WDy9hzTWDDft8h76ZWYfYsX+cjY8e5NS7ZwEYP3GKjY8eBGhY8Ht6x8ysQ9z19JHfBX7JqXfPctfTRxp2Doe+mVmHOHbiVF3tM+HQNzPrEPMH+utqnwmHvplZh1g/vIT+vjkfaOvvm8P64SUNO4cv5JqZdYjSxVqv3jEzy4k11ww2NOSn8vSOmVmOOPTNzHLEoW9mliPThr6kCyTtlfSSpEOSvp7a75A0LulAelxfdsxGSUclHZE0XNZ+raSD6b17Jak5X8vMzCqp5ULuaeC6iHhHUh/wnKSn0nv3RMQ3yneWdCWwFlgKzAd+JOkTEXEWuB8YAZ4HngRWAU9hZmYtMe1IP4reSS/70iOqHLIa2BYRpyPiDeAosFzSPOCiiNgTEQE8DKyZVe/NzKwuNc3pS5oj6QBwHNgVES+kt26T9LKkByVdnNoGgbfKDh9LbYNpe2p7pfONSBqVNDo5OVn7tzEz6xI79o+zcvNuFm14gpWbd7Nj/3hLzltT6EfE2YhYBgxRHLVfRXGq5uPAMmAC+GbavdI8fVRpr3S+LRFRiIjC3Lnn3MzdzKyrlappjp84RfB+Nc1WBH9dq3ci4gTwY2BVRLyd/jJ4D/gWsDztNgYsKDtsCDiW2ocqtJuZ5UorqmlmqWX1zlxJA2m7H/g08Fqaoy/5PPBK2t4JrJV0vqRFwGJgb0RMACclrUirdm4CHmvcVzEz6w6tqKaZpZbVO/OArZLmUPxLYntEPC7p3yUtozhF8ybwZYCIOCRpO/AqcAZYl1buANwKPAT0U1y145U7ZpY78wf6Ga8Q8I2spplFxYU0natQKMTo6Gi7u2Fm1jBT75AFxWqad95wdcPq7kjaFxGFqe0uuGZm1mKtqKaZxaFvZtYGza6mmcW1d8zMcsShb2aWIw59M7McceibmeWIQ9/MLEe8esfMrIl27B9vy9LMLA59M7MmmfojrFJhNaBtwe/pHTOzJmlnYbUsDn0zsyZpZ2G1LA59M7MmySqg1orCalkc+mZmTbJ+eAn9fXM+0NbfN4f1w0va1CNfyDUza5p2FlbL4tA3M2uidhVWy+LpHTOzHHHom5nliEPfzCxHarkx+gWS9kp6SdIhSV9P7ZdI2iXp9fR8cdkxGyUdlXRE0nBZ+7WSDqb37k03SDczsxapZaR/GrguIv4QWAaskrQC2AA8ExGLgWfSayRdCawFlgKrgPvSTdUB7gdGgMXpsapxX8XMzKYzbehH0TvpZV96BLAa2JratwJr0vZqYFtEnI6IN4CjwHJJ84CLImJPFO/G/nDZMWZm1gI1LdlMI/V9wO8D/xYRL0i6PCImACJiQtJlafdB4Pmyw8dS27tpe2p7pfONUPw/Aq644orav42ZWZt0WjXNLDVdyI2IsxGxDBiiOGq/qsrulebpo0p7pfNtiYhCRBTmzp1bSxfNzNqmVE1z/MQpgverae7YP97urp2jrtU7EXEC+DHFufi305QN6fl42m0MWFB22BBwLLUPVWg3M+tqnVhNM0stq3fmShpI2/3Ap4HXgJ3AzWm3m4HH0vZOYK2k8yUtonjBdm+aCjopaUVatXNT2TFmZl2rE6tpZqllTn8esDXN658HbI+IxyXtAbZLugX4BfAFgIg4JGk78CpwBlgXEaW/Am8FHgL6gafSw8ysq80f6Ge8QsC3s5pmFhUX0nSuQqEQo6Oj7e6GmVmmqXfIgmI1zTtvuLptF3Ml7YuIwtR2F1wzM5ulTqymmcWhb2bWAJ1WTTOLa++YmeWIQ9/MLEcc+mZmOeLQNzPLEYe+mVmOePWOmVmNuqWoWjUOfTOzGkz9AVapqBrQVcHv6R0zsxp0U1G1ahz6ZmY16KaiatU49M3MapBVPK0Ti6pV49A3M6vB+uEl9PfN+UBbf98c1g8vaVOPZsYXcs3MatBNRdWqceibmdWoW4qqVePpHTOzHHHom5nliEPfzCxHarkx+gJJz0o6LOmQpNtT+x2SxiUdSI/ry47ZKOmopCOShsvar5V0ML13b7pBupmZtUgtF3LPAF+NiBclfRTYJ2lXeu+eiPhG+c6SrgTWAkuB+cCPJH0i3Rz9fmAEeB54EliFb45uZtYy0470I2IiIl5M2yeBw0C1y9ergW0RcToi3gCOAsslzQMuiog9Ubwb+8PAmtl+ATMzq11dc/qSFgLXAC+kptskvSzpQUkXp7ZB4K2yw8ZS22Dantpe6TwjkkYljU5OTtbTRTOzWduxf5yVm3ezaMMTrNy8mx37x9vdpYapOfQlfQR4BPhKRPyW4lTNx4FlwATwzdKuFQ6PKu3nNkZsiYhCRBTmzp1baxfNzGatVE1z/MQpgverafZK8NcU+pL6KAb+dyLiUYCIeDsizkbEe8C3gOVp9zFgQdnhQ8Cx1D5Uod3MrGP0SjXNLLWs3hHwAHA4Iu4ua59XttvngVfS9k5graTzJS0CFgN7I2ICOClpRfrMm4DHGvQ9zMwaoleqaWapZfXOSuBLwEFJB1Lb14AbJS2jOEXzJvBlgIg4JGk78CrFlT/r0sodgFuBh4B+iqt2vHLHzDrK/IF+xisEfLdV08wybehHxHNUno9/ssoxm4BNFdpHgavq6aCZWSutH17ygTtkQXdW08zigmtmZmV6pZpmFoe+mdkUvVBNM4tr75iZ5YhD38wsRxz6ZmY54tA3M8sRh76ZWY549Y6Z5daO/eM9uzQzi0PfzHKpVFit9COsUmE1oKeD39M7ZpZLvV5YLYtD38xyqdcLq2Vx6JtZLmUVUOuVwmpZHPpmlkvrh5fQ3zfnA229VFgtiy/kmlku9XphtSwOfTPLrV4urJbF0ztmZjni0Dczy5Fa7pG7QNKzkg5LOiTp9tR+iaRdkl5PzxeXHbNR0lFJRyQNl7VfK+lgeu/edK9cMzNrkVpG+meAr0bEHwArgHWSrgQ2AM9ExGLgmfSa9N5aYCmwCrhPUukS+f3ACMWbpS9O75uZWYtMG/oRMRERL6btk8BhYBBYDWxNu20F1qTt1cC2iDgdEW8AR4HlkuYBF0XEnogI4OGyY8zMrAXqmtOXtBC4BngBuDwiJqD4FwNwWdptEHir7LCx1DaYtqe2VzrPiKRRSaOTk5P1dNHMzKqoecmmpI8AjwBfiYjfVpmOr/RGVGk/tzFiC7AFoFAoVNzHzKxWeaymmaWm0JfURzHwvxMRj6bmtyXNi4iJNHVzPLWPAQvKDh8CjqX2oQrtZmZNk9dqmllqWb0j4AHgcETcXfbWTuDmtH0z8FhZ+1pJ50taRPGC7d40BXRS0or0mTeVHWNm1hR5raaZpZaR/krgS8BBSQdS29eAzcB2SbcAvwC+ABARhyRtB16luPJnXUSU/onfCjwE9ANPpYeZWdPktZpmlmlDPyKeo/J8PMCnMo7ZBGyq0D4KXFVPB83MZmP+QD/jFQK+16tpZvEvcs2sp+W1mmYWF1wzs56W12qaWRz6Ztbz8lhNM4und8zMcsShb2aWIw59M7McceibmeWIQ9/MLEcc+mZmOeIlm2bWM1xNc3oOfTPrCa6mWRtP75hZT3A1zdo49M2sJ7iaZm0c+mbWE7KqZua1mmYWh76Z9QRX06yNL+SaWU9wNc3aOPTNrGe4mub0PL1jZpYjtdwY/UFJxyW9UtZ2h6RxSQfS4/qy9zZKOirpiKThsvZrJR1M792bbo5uZmYtVMtI/yFgVYX2eyJiWXo8CSDpSmAtsDQdc5+k0pWV+4ERYHF6VPpMMzNromlDPyJ+Avyqxs9bDWyLiNMR8QZwFFguaR5wUUTsiYgAHgbWzLDPZmY2Q7OZ079N0stp+ufi1DYIvFW2z1hqG0zbU9srkjQiaVTS6OTk5Cy6aGZm5WYa+vcDHweWARPAN1N7pXn6qNJeUURsiYhCRBTmzp07wy6aWS/asX+clZt3s2jDE6zcvJsd+8fb3aWuMqMlmxHxdmlb0reAx9PLMWBB2a5DwLHUPlSh3cysZi6qNnszGumnOfqSzwOllT07gbWSzpe0iOIF270RMQGclLQirdq5CXhsFv02sxxyUbXZm3akL+l7wCeBSyWNAf8EfFLSMopTNG8CXwaIiEOStgOvAmeAdRFR+hO6leJKoH7gqfQwM6uZi6rN3rShHxE3Vmh+oMr+m4BNFdpHgavq6p2ZWZn5A/2MVwh4F1WrnX+Ra2Zdw0XVZs+1d8ysa7io2uw59M2sq7io2ux4esfMLEcc+mZmOeLQNzPLEYe+mVmOOPTNzHLEoW9mliNesmlmHWnH/nGvx28Ch76ZdRxX02weT++YWcdxNc3mceibWcdxNc3mceibWcfJqprpapqz59A3s47japrN4wu5ZtZxXE2zeRz6ZtaRXE2zOaad3pH0oKTjkl4pa7tE0i5Jr6fni8ve2yjpqKQjkobL2q+VdDC9d2+6V66ZmbVQLXP6DwGrprRtAJ6JiMXAM+k1kq4E1gJL0zH3SSpNzN0PjFC8WfriCp9pZmZNNm3oR8RPgF9NaV4NbE3bW4E1Ze3bIuJ0RLwBHAWWS5oHXBQReyIigIfLjjEzsxaZ6eqdyyNiAiA9X5baB4G3yvYbS22DaXtqu5mZtVCjl2xWmqePKu2VP0QakTQqaXRycrJhnTMzy7uZrt55W9K8iJhIUzfHU/sYsKBsvyHgWGofqtBeUURsAbYAFAqFzL8czKz7ubBaa810pL8TuDlt3ww8Vta+VtL5khZRvGC7N00BnZS0Iq3auansGDPLqVJhtfETpwjeL6y2Y/94u7vWs2pZsvk9YA+wRNKYpFuAzcBnJL0OfCa9JiIOAduBV4EfAOsiolQ16Vbg2xQv7v438FSDv4uZdRkXVmu9aad3IuLGjLc+lbH/JmBThfZR4Kq6emdmPc2F1VrPtXfMrG1cWK31HPpm1jYurNZ6rr1jZm3jwmqt59A3s7ZyYbXW8vSOmVmOOPTNzHLEoW9mliMOfTOzHHHom5nliEPfzCxHvGTTzFrC1TQ7g0PfzJquVE2zVFytVE0TcPC3mKd3zKzpXE2zczj0zazpXE2zczj0zazpXE2zc3hO38waqtIF2/XDSz4wpw+uptkuHumbWcNk3f4Q4M4brmZwoB8BgwP93HnD1b6I2wYe6ZtZ3bKWX1a7YPvTDdc55DvArEJf0pvASeAscCYiCpIuAf4TWAi8Cfx1RPw67b8RuCXt/3cR8fRszm9mrVdt+aUv2Ha+Roz0/ywifln2egPwTERslrQhvf5HSVcCa4GlwHzgR5I+UXbjdDPrMJVG9NVG8/MH+hmvEPC+YNs5mjGnvxrYmra3AmvK2rdFxOmIeAM4CixvwvnNrAGy5ucrhToUR/O+/WHnm23oB/BDSfskjaS2yyNiAiA9X5baB4G3yo4dS23nkDQiaVTS6OTk5Cy7aGYzkTWinyNV3H/+QD9rrhn0BdsON9vpnZURcUzSZcAuSa9V2bfSvylRaceI2AJsASgUChX3MbPGqTSNkzUPfzaC/r45mcsvffvDzjarkX5EHEvPx4HvU5yueVvSPID0fDztPgYsKDt8CDg2m/Ob2exlTeMMfLiv4v6l0btH891pxiN9SRcC50XEybT958A/AzuBm4HN6fmxdMhO4LuS7qZ4IXcxsHcWfTezOtVzYfb8D52XOaL3aL57zWZ653Lg+yrO730I+G5E/EDSz4Dtkm4BfgF8ASAiDknaDrwKnAHWeeWOWetkLbWcGvglvzn1Lvd8cZnLIfcYRXT2lHmhUIjR0dF2d8Osa2T9cGrl5t0VV97MkThbIQcGB/r56YbrWtFlawJJ+yKiMLXdv8g16yEz+eHUdBdmrbc49M26VKN+ODVYdqyncXqfp3fMutDUET1wzmi9nIB7vris4jFeedObPL1j1kOq/XCq0vx86YdTpWM9os8vh75Zh/MPp6yRXE/frIP5h1PWaB7pm3UI/3DKWsEjfbMOUG9Fy9+cetcjepsRj/TNWqyeEf10F2Yd8lYvh75ZC9VbCsE/nLJG8/SOWRPs2D/Oys27WbThCVZu3s2O/eNA/TXqfWHWGs0jfbMGa3QpBE/jWCN5pG/WYNOVQqjEI3prFY/0zWahnh9OHTtxKrMUgkf01ioOfbMyWWWJK7UDFadxBj7cx6//991zPtulEKwTuOCadb16groUrrWEOBRH4X917SCP7Bs/p/2CvvMqhvtAfx+nz7znwmbWVlkF1xz61lT1Bm+1sG5EUN95w9WZx2SFeNZa+SylipYezVs7OfStLs0M5GrBmxXW9Y62q90NCsj8pWsj+I5T1gk6JvQlrQL+FZgDfDsiNlfbfyah34jAauRotNvOnVWrvVGBXC14s8K63tF2ltJq+Ho+KevcnsaxTtYRoS9pDvBfwGeAMeBnwI0R8WrWMfWGfqMCq5Gj0W47911PH2lqIM8keOs1k5F+VohP9+fkaRzrRJ0S+n8C3BERw+n1RoCIuDPrmHpDv96bP88kHOr9rG4797FU9KtZGvn9ZhrU1aad6v0/MrNO1Cl3zhoE3ip7PQb88dSdJI0AIwBXXHFFXSeo9ovHetqzPmcmn9Vt5866l2qjAnkmF2Cz2u/43FKgclAXPnZJ1aDOeq9SmHsNvfWKVod+pQIj56RFRGwBtkBxpF/PCRoVWPPbONpu97nXDy9peiCX1BPW1UK83qB2iFte9dz0juf0Z3/uRl+QNrPW65TpnZ8BiyUtAsaBtcDfNPIE1X7xOJMRZKM+q9vOnTUSrrfdzDpLO5ZsXg/8C8Ulmw9GxKZq+3udvplZ/TplpE9EPAk82erzmpmZSyubmeWKQ9/MLEcc+mZmOeLQNzPLkY6vsilpEvj5DA+/FPhlA7vTDfyd8yFv3zlv3xdm/50/FhFzpzZ2fOjPhqTRSkuWepm/cz7k7Tvn7ftC876zp3fMzHLEoW9mliO9Hvpb2t2BNvB3zoe8fee8fV9o0nfu6Tl9MzP7oF4f6ZuZWRmHvplZjvRk6EtaJemIpKOSNrS7P60g6UFJxyW90u6+tIKkBZKelXRY0iFJt7e7T80m6QJJeyW9lL7z19vdp1aRNEfSfkmPt7svrSDpTUkHJR2Q1NAywz03pz+Tm6/3Akl/CrwDPBwRV7W7P80maR4wLyJelPRRYB+wppf/nCUJuDAi3pHUBzwH3B4Rz7e5a00n6e+BAnBRRHy23f1pNklvAoWIaPgP0npxpL8cOBoR/xMR/wdsA1a3uU9NFxE/AX7V7n60SkRMRMSLafskcJjiPZh7VhS9k172pUdvjdoqkDQE/CXw7Xb3pRf0YuhXuvl6T4dB3klaCFwDvNDmrjRdmuY4ABwHdkVEz39nijdd+gfgvTb3o5UC+KGkfZJGGvnBvRj6Nd183XqDpI8AjwBfiYjftrs/zRYRZyNiGTAELJfU01N5kj4LHI+Ife3uS4utjIg/Av4CWJembxuiF0N/DFhQ9noIONamvlgTpXntR4DvRMSj7e5PK0XECeDHwKr29qTpVgKfS3Pc24DrJP1He7vUfBFxLD0fB75Pcdq6IXox9H9383VJv0fx5us729wna7B0UfMB4HBE3N3u/rSCpLmSBtJ2P/Bp4LW2dqrJImJjRAxFxEKK/y3vjoi/bXO3mkrShWlxApIuBP4caNiqvJ4L/Yg4A9wGPE3x4t72iDjU3l41n6TvAXuAJZLGJN3S7j412UrgSxRHfgfS4/p2d6rJ5gHPSnqZ4uBmV0TkYgljzlwOPCfpJWAv8ERE/KBRH95zSzbNzCxbz430zcwsm0PfzCxHHPpmZjni0DczyxGHvplZjjj0zcxyxKFvZpYj/w/FlujlbhUtZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x_high_density = np.linspace(0.0, 5, 50)\n",
    "\n",
    "plt.scatter(x_high_density, y_linear(x_high_density))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also the possibility to use a quadratic formula for interpolation, in which case the function *and* its first derivative are continuous at each point. The keyword `kind=2` specifies that we want a polynomial of second order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1ce39561f40>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD4CAYAAAAAczaOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAXf0lEQVR4nO3db4xc1X3G8e+D45otCdogNsjetWMrctwCVu0wcqksRflvl0ax4zaJqQq8oNoIGYWoERXOm5IXkZFIQssLkJwEYdoUxxIErCTEIUAUpTKYMTYxtuPGLSTs2sKOEjdGtVxsfn0xZ2FY5s7O7Py5M3OfjzTamTP3zpyb4Gfv/u655ygiMDOzYrgg7w6YmVn3OPTNzArEoW9mViAOfTOzAnHom5kVyDvy7sBMLr300li8eHHe3TAz6yt79+79bUSMTG/v+dBfvHgx5XI5726YmfUVSb+u1e7yjplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFUjPj94xMyuSR/ZNcueuIxw7dYYFw0PcumYZ61eOtu3zHfpmZj3ikX2TbH74AGdeOw/A5KkzbH74AEDbgt/lHTOzHnHnriNvBP6UM6+d585dR9r2HQ59M7MecezUmabaZ8Ohb2bWIxYMDzXVPhsOfTOzHnHrmmUMzZ3zlrahuXO4dc2ytn2HL+SamfWIqYu1Hr1jZlYQ61eOtjXkp3N5x8ysQBz6ZmYF4tA3MyuQGUNf0oWS9kh6XtJBSV9J7bdLmpS0Pz2uqdpns6Sjko5IWlPVfpWkA+m9uyWpM4dlZma1NHIh9yzwkYh4VdJc4OeSHkvv3RURX6veWNLlwEbgCmAB8BNJ74+I88C9wDjwNPBDYC3wGGZm1hUznulHxavp5dz0iDq7rAO2R8TZiHgROAqskjQfuDgidkdEAA8A61vqvZmZNaWhmr6kOZL2AyeAxyPimfTWzZJ+Iek+Se9ObaPAy1W7T6S20fR8enut7xuXVJZUPnnyZONHY2ZmdTUU+hFxPiJWAGNUztqvpFKqeR+wAjgOfD1tXqtOH3Xaa33f1ogoRURpZORti7mbmdksNTV6JyJOAT8F1kbEK+mXwevAN4FVabMJYGHVbmPAsdQ+VqPdzMy6pJHROyOShtPzIeBjwC9TjX7Kp4EX0vOdwEZJ8yQtAZYCeyLiOHBa0tVp1M71wKPtOxQzM5tJI6N35gPbJM2h8ktiR0R8X9K/SlpBpUTzEvB5gIg4KGkHcAg4B2xKI3cAbgLuB4aojNrxyB0zsy5SZSBN7yqVSlEul/PuhplZX5G0NyJK09s94ZqZWQ46vRZuFoe+mVmXdWMt3Cyee8fMrMu6sRZuFoe+mVmXdWMt3CwOfTOzLuvGWrhZHPpmZl3WjbVws/hCrplZl3VjLdwsDn0zsxx0ei3cLC7vmJkViEPfzKxAHPpmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYhD38ysQBpZGP1CSXskPS/poKSvpPZLJD0u6Vfp57ur9tks6aikI5LWVLVfJelAeu/utEC6mdnAemTfJKvveJIlt/2A1Xc8ySP7JnPtTyNn+meBj0TEnwErgLWSrgZuA56IiKXAE+k1ki4HNgJXAGuBe9Ki6gD3AuPA0vRY275DMTPrLVMrZE2eOkPw5gpZeQb/jKEfFa+ml3PTI4B1wLbUvg1Yn56vA7ZHxNmIeBE4CqySNB+4OCJ2R2U19geq9jEzGzh5rpCVpaGavqQ5kvYDJ4DHI+IZ4LKIOA6Qfr4nbT4KvFy1+0RqG03Pp7fX+r5xSWVJ5ZMnTzZxOGZmvSPPFbKyNBT6EXE+IlYAY1TO2q+ss3mtOn3Uaa/1fVsjohQRpZGRkUa6aGbWc/JcIStLU6N3IuIU8FMqtfhXUsmG9PNE2mwCWFi12xhwLLWP1Wg3MxtIea6QlaWR0TsjkobT8yHgY8AvgZ3ADWmzG4BH0/OdwEZJ8yQtoXLBdk8qAZ2WdHUatXN91T5mZgNn/cpRtmxYzujwEAJGh4fYsmF5LounTGlk5az5wLY0AucCYEdEfF/SbmCHpBuB3wCfAYiIg5J2AIeAc8CmiJi6knETcD8wBDyWHmZmAyuvFbKyqDKQpneVSqUol8t5d8PMrK9I2hsRpentviPXzKxAHPpmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFUgjd+SamdkMHtk3yZ27jnDs1BkWDA9x65plPXUn7hSHvplZi6YWS5maO39qsRSg54Lf5R0zsxb14mIpWRz6ZmYt6sXFUrI49M3MWtSLi6VkceibmbWoFxdLyeILuWZmLZq6WOvRO2ZmBdFri6VkcXnHzKxAGlkjd6GkpyQdlnRQ0i2p/XZJk5L2p8c1VftslnRU0hFJa6rar5J0IL13d1or18zMuqSR8s454EsR8ZykdwF7JT2e3rsrIr5WvbGky4GNwBXAAuAnkt6f1sm9FxgHngZ+CKzF6+SamXXNjGf6EXE8Ip5Lz08Dh4F6hat1wPaIOBsRLwJHgVWS5gMXR8TuqCzM+wCwvtUDMDOzxjVV05e0GFgJPJOabpb0C0n3SXp3ahsFXq7abSK1jabn09vNzKxLGg59Se8EHgK+GBF/oFKqeR+wAjgOfH1q0xq7R532Wt81LqksqXzy5MlGu2hmZjNoKPQlzaUS+N+JiIcBIuKViDgfEa8D3wRWpc0ngIVVu48Bx1L7WI32t4mIrRFRiojSyMhIM8djZmZ1NDJ6R8C3gcMR8Y2q9vlVm30aeCE93wlslDRP0hJgKbAnIo4DpyVdnT7zeuDRNh2HmZk1oJHRO6uB64ADkvanti8D10paQaVE8xLweYCIOChpB3CIysifTWnkDsBNwP3AEJVROx65Y2bWRaoMpOldpVIpyuVy3t0wM+ubhVIAJO2NiNL0dk/DYGbWgH5aKKUeT8NgZtaAfloopR6HvplZA/ppoZR6HPpmZg3op4VS6nHom5k1oJ8WSqnHF3LNzBrQTwul1OPQNzNrUL8slFKPyztmZgXi0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYjH6ZuZTdNPUyg3y6FvZlZlUKZQzuLyjplZlUGZQjmLQ9/MrMqgTKGcpZGF0RdKekrSYUkHJd2S2i+R9LikX6Wf767aZ7Oko5KOSFpT1X6VpAPpvbvTAulmZj1jUKZQztLImf454EsR8afA1cAmSZcDtwFPRMRS4In0mvTeRuAKYC1wj6Sp+UjvBcaBpemxto3HYmbWskGZQjnLjKEfEccj4rn0/DRwGBgF1gHb0mbbgPXp+Tpge0ScjYgXgaPAKknzgYsjYndUVmN/oGofM7OesH7lKFs2LGd0eAgBo8NDbNmwfCAu4kKTo3ckLQZWAs8Al0XEcaj8YpD0nrTZKPB01W4Tqe219Hx6e63vGafyFwGLFi1qpotmZi0bhCmUszR8IVfSO4GHgC9GxB/qbVqjLeq0v70xYmtElCKiNDIy0mgXzcxsBg2FvqS5VAL/OxHxcGp+JZVsSD9PpPYJYGHV7mPAsdQ+VqPdzMy6pJHROwK+DRyOiG9UvbUTuCE9vwF4tKp9o6R5kpZQuWC7J5WCTku6On3m9VX7mJlZFzRS018NXAcckLQ/tX0ZuAPYIelG4DfAZwAi4qCkHcAhKiN/NkXE1J0ONwH3A0PAY+lhZmZdospAmt5VKpWiXC7n3Q0zs74iaW9ElKa3+45cM7MC8YRrZlZYgzybZhaHvpkV0qDPppnF5R0zK6RBn00zi0PfzApp0GfTzOLQN7NCGvTZNLM49M2skAZ9Ns0svpBrZoU0dbHWo3fMzApikGfTzOLyjplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFYhH75jZwCvixGpZHPpmNtCKOrFaFpd3zGygFXVitSwOfTMbaEWdWC1LIwuj3yfphKQXqtpulzQpaX96XFP13mZJRyUdkbSmqv0qSQfSe3enxdHNzDqqqBOrZWnkTP9+YG2N9rsiYkV6/BBA0uXARuCKtM89kqZmNLoXGAeWpketzzQza6uiTqyWZcbQj4ifAb9r8PPWAdsj4mxEvAgcBVZJmg9cHBG7o7IS+wPA+ln22cysYetXjrJlw3JGh4cQMDo8xJYNywt5ERdaG71zs6TrgTLwpYj4PTAKPF21zURqey09n95ek6RxKn8VsGjRoha6aGZWzInVssz2Qu69wPuAFcBx4OupvVadPuq01xQRWyOiFBGlkZGRWXbRzMymm1XoR8QrEXE+Il4HvgmsSm9NAAurNh0DjqX2sRrtZmbWRbMK/VSjn/JpYGpkz05go6R5kpZQuWC7JyKOA6clXZ1G7VwPPNpCv83MbBZmrOlLehD4EHCppAngn4APSVpBpUTzEvB5gIg4KGkHcAg4B2yKiKm7Im6iMhJoCHgsPczMrItUGUzTu0qlUpTL5by7YWZ9wHPsvEnS3ogoTW/33DtmNhA8x05jPA2DmQ0Ez7HTGIe+mQ0Ez7HTGIe+mQ0Ez7HTGIe+mQ0Ez7HTGF/INbOBMHWx1qN36nPom9nA8Bw7M3N5x8ysQBz6ZmYF4tA3MysQ1/TNrK94qoXWOPTNrG94qoXWubxjZn3DUy20zqFvZn3DUy20zqFvZn3DUy20zqFvZn3DUy20zhdyzaxveKqF1jWyXOJ9wCeBExFxZWq7BPgusJjKcomfjYjfp/c2AzcC54EvRMSu1H4Vby6X+EPgluj1ZbvMrOd4qoXWNFLeuR9YO63tNuCJiFgKPJFeI+lyYCNwRdrnHklTf4vdC4xTWSx9aY3PNDOzDpsx9CPiZ8DvpjWvA7al59uA9VXt2yPibES8CBwFVkmaD1wcEbvT2f0DVfuYmVmXzLamf1lEHAeIiOOS3pPaR4Gnq7abSG2vpefT22uSNE7lrwIWLVo0yy6aWT/znbed0e7RO6rRFnXaa4qIrRFRiojSyMhI2zpnZv1h6s7byVNnCN688/aRfZN5d63vzTb0X0klG9LPE6l9AlhYtd0YcCy1j9VoNzN7G9952zmzDf2dwA3p+Q3Ao1XtGyXNk7SEygXbPakUdFrS1ZIEXF+1j5nZW/jO286ZMfQlPQjsBpZJmpB0I3AH8HFJvwI+nl4TEQeBHcAh4EfApoiY+nV9E/AtKhd3/wt4rM3HYmYDwnfeds6MF3Ij4tqMtz6asf1Xga/WaC8DVzbVOzMrpFvXLHvLbJrgO2/bxXfkmlnP8Z23nePQN7Oe5DtvO8MTrpmZFYjP9M0sV74Jq7sc+maWGy9/2H0u75hZbnwTVvc59M0sN74Jq/sc+maWG9+E1X0OfTPLjZc/7D5fyDWz3PgmrO5z6JtZrnwTVnc59M2sKzwevzc49M2s4zwev3f4Qq6ZdZzH4/cOh76ZdZzH4/cOh76ZdZzH4/cOh76ZdZzH4/cOX8g1s47zePze0VLoS3oJOA2cB85FREnSJcB3gcXAS8BnI+L3afvNwI1p+y9ExK5Wvt/Mek/W0EyPx+8N7SjvfDgiVkREKb2+DXgiIpYCT6TXSLoc2AhcAawF7pE0p9YHmll/mhqaOXnqDMGbQzMf2TeZd9cs6URNfx2wLT3fBqyvat8eEWcj4kXgKLCqA99vZjnx0Mze12roB/BjSXsljae2yyLiOED6+Z7UPgq8XLXvRGp7G0njksqSyidPnmyxi2bWLR6a2ftaDf3VEfEB4C+BTZI+WGdb1WiLWhtGxNaIKEVEaWRkpMUumlm3eGhm72sp9CPiWPp5AvgelXLNK5LmA6SfJ9LmE8DCqt3HgGOtfL+Z9RYPzex9sw59SRdJetfUc+ATwAvATuCGtNkNwKPp+U5go6R5kpYAS4E9s/1+M8vPI/smWX3Hkyy57QesvuPJNy7Url85ypYNyxkdHkLA6PAQWzYs96idHtLKkM3LgO9Jmvqcf4+IH0l6Ftgh6UbgN8BnACLioKQdwCHgHLApIs7X/mgz61UzTZ7moZm9TRE1y+o9o1QqRblczrsbZpasvuNJJmtcmB0dHuI/bvtIDj2yWiTtrRpK/wZPw2BmTfEInf7maRjMLFOtu2sXDA/VPNP3CJ3+4DN9M6sp6+7aD//JiEfo9DGHvpnVlHV37VO/POkROn3M5R0zq6le7d4jdPqXQ9/MXLsvEJd3zArOtfticeibFZxr98Xi8o5ZgdQq47h2XywOfbOCyJo+YfiP5/L7/33tbdu7dj+YHPpmfSprWcKs9qwyzrx3XMDQ3Dlvec+1+8Hl0DfrcbVCHKh51l7+9e94aO9kzcnQsso4/3PmNe763AovWl4QnnDNrMuaOUOHt4Y7VM7CL5x7Qc2SzByJ8zX+TY+mUo0nSiuOrAnXfKZv1oJmSyxZdfWsM/QL515QsyQzvW1KrcCHyln+XZ9bUfMXiMs4xeLQN6vSyll4IyWWrLr6g8+8/LbArhfuWbLO9BcMD71RrnEZp9hc3rGBlRXgWe9B7VLKX181+pYQn2qfTYnlWLoBqlXDQ3M5e+71hvvq8fXFk1XecehbbmZTGmm1Fr5lw/LM95oN8WYJMqc2yPqOrHCfOo5m/neyYumZ0Je0FvgXYA7wrYi4o972Dv3OafZMuF2BXKu2DfXPVJttzwrwehc026Xemf6ta5Y1ddz1wt2snp4IfUlzgP8EPg5MAM8C10bEoax9ZhP67Qqmdn5Wr303NH8m3K5A3rJhOXfuOtLUGW+z7VmUfjbzX32zZ+EzlVhm89+CWbN6JfT/Arg9Itak15sBImJL1j7Nhn67ziC7EX55fvdszoTbFcjtrG03q97xzSbEwSUW6029Evp/A6yNiL9Pr68D/jwibs7ap9nQz1q0eTbBBJ0Nvzy/O8tszoSbNZvadrPtM9XC6/2V4xC3QdAr4/RVo+1t/2IljQPjAIsWLWrqC7LuOswKvnrjmrM0+1m9+N1ZFnThF86CWdS2m22//VNXAPVr4Vnv1QpzTzxmg6LboT8BLKx6PQYcm75RRGwFtkLlTL+ZL2jXGWQ3wi/P7846E57N0MVm26sDtlbwlt57SVvaoXaAT7U7xK2Iul3eeQeVC7kfBSapXMj924g4mLWPa/qd/Y68LmCbWWf1RE0/deQa4J+pDNm8LyK+Wm97j95x8JpZ83om9JvlcfpmZs3LCn0vl2hmViAOfTOzAnHom5kViEPfzKxAHPpmZgXS86N3JJ0Efj3L3S8FftvG7vQDH3MxFO2Yi3a80PoxvzciRqY39nzot0JSudaQpUHmYy6Goh1z0Y4XOnfMLu+YmRWIQ9/MrEAGPfS35t2BHPiYi6Fox1y044UOHfNA1/TNzOytBv1M38zMqjj0zcwKZCBDX9JaSUckHZV0W9796QZJ90k6IemFvPvSDZIWSnpK0mFJByXdknefOk3ShZL2SHo+HfNX8u5Tt0iaI2mfpO/n3ZdukPSSpAOS9ktq6zTDA1fTlzSHykItH6eyUtezwLURcSjXjnWYpA8CrwIPRMSVefen0yTNB+ZHxHOS3gXsBdYP8v/PkgRcFBGvSpoL/By4JSKezrlrHSfpH4AScHFEfDLv/nSapJeAUkS0/Ya0QTzTXwUcjYj/joj/A7YD63LuU8dFxM+A3+Xdj26JiOMR8Vx6fho4DAz0yjBR8Wp6OTc9BuusrQZJY8BfAd/Kuy+DYBBDfxR4uer1BAMeBkUnaTGwEngm5650XCpz7AdOAI9HxMAfM5WV9v4ReD3nfnRTAD+WtFfSeDs/eBBDXzXaBv5sqKgkvRN4CPhiRPwh7/50WkScj4gVwBiwStJAl/IkfRI4ERF78+5Ll62OiA8AfwlsSuXbthjE0J8AFla9HgOO5dQX66BU134I+E5EPJx3f7opIk4BPwXW5tuTjlsNfCrVuLcDH5H0b/l2qfMi4lj6eQL4HpWydVsMYug/CyyVtETSHwEbgZ0598naLF3U/DZwOCK+kXd/ukHSiKTh9HwI+Bjwy1w71WERsTkixiJiMZV/y09GxN/l3K2OknRRGpyApIuATwBtG5U3cKEfEeeAm4FdVC7u7YiIg/n2qvMkPQjsBpZJmpB0Y9596rDVwHVUzvz2p8c1eXeqw+YDT0n6BZWTm8cjohBDGAvmMuDnkp4H9gA/iIgftevDB27IppmZZRu4M30zM8vm0DczKxCHvplZgTj0zcwKxKFvZlYgDn0zswJx6JuZFcj/A5y91hp5NBxHAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_quadric = interpolate.interp1d(x, y, kind=2)\n",
    "\n",
    "x_high_density = np.linspace(0.0, 5, 50)\n",
    "\n",
    "plt.scatter(x_high_density, y_quadric(x_high_density))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A nice little hack to know is, that if you function is bijective, you can calculate the inverse with the inverse interpolation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x1ce395cdc10>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAD4CAYAAADFAawfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAARNElEQVR4nO3dcYwU53nH8d/DceALBhHKYdnYlAQhVDclxj05RFRRGoXguFFCVVVNa7f5Iyr/JJJdaiojW7Vd1VJTWpRWrapCY9VtnMSOgq9WlZigKshKRHCO3OEzxVfbqdN6scxZFjVuTw7GT//Ywb47Zm5n72Z2nt35fqTVDe/O7T6MuB9z8z6zr7m7AABxLaq6AADA3AhqAAiOoAaA4AhqAAiOoAaA4BaX8aKrV6/29evXl/HSANCTTpw48aq7D6Y9V0pQr1+/XiMjI2W8NAD0JDP7adZzXPoAgOAIagAIjqAGgOAIagAIjqAGgOBydX2Y2YuSzku6KOktdx8qsygA6Bb3DI/rqz/8rxlj2zas0sO//+HC3qOd9rxfdfdXC3tnAOgytx48ph+88FrL/X7wwmu69eCxwsK6lD5qAOh22/cf1XNn/3fe358n0PPKG9Qu6btm5pL+3t0PzN7BzHZJ2iVJ69atK6xAACjT8GhDe745pgtvV11JtrxBvc3dz5jZGklHzOxZd39y+g5JeB+QpKGhIVYjABDO8GhDux8ZU+BMTpUrqN39TPL1rJk9JukmSU/O/V0AUJ20Sb5O2rZhVWGv1TKozWyZpEXufj7Z/oSkPymsAgBYoLyTfJ2ycc2yjnd9XCXpMTO7tP/X3P2JwioAgDZEC+Xpli5epC/9xmbt3LK20NdtGdTu/hNJHyz0XQEgh4V2XpRp45plOrL7ox15L9rzAIQQ+Uy56BtY2kVQA+io6O1wt21dpz/d+UtVlzEDQQ2gNNHb4SKGchqCGkAhqm6Hm8viRaa/+M0PFj7J1ykENYC2RQ7lsjovqkRQA5gTk3zVI6gBvCNyO1xdQjkNQQ3UUPTOizqHchqCGuhxka8nS4RyHgQ10EOih3K3tMNFQ1ADXSryJF+3t8NFQ1ADXSByKPdiO1w0BDUQSPRJvk5+EBHeRVADFYkeykzyxUFQAx3AJB8WgqAGCkYoo2gENbAATPKhEwhqIKfIt1cTyr2NoAZmYZIP0RDUqLXoH2xPKEMiqFEjTPKhWxHU6EmRJ/kkQhntIajR9SKHMpN8KAJBja4SOZS5vRplIagRVuRQZpIPnURQo3LR2+G4noyqEdToqOjtcIQyIiKoUZrI7XB8sD26CUGNQkQOZTov0O0IarSNST6gswhqzIkPIgKqlzuozaxP0oikhrt/qrySUIXonRecKaPO2jmjvl3SaUkrSqoFHRL5erJEKAOz5QpqM7tW0q9JekDS7lIrQqGihzLtcEBrec+ovyzpjyQtz9rBzHZJ2iVJ69atW3BhaF/kST7a4YD5axnUZvYpSWfd/YSZfTRrP3c/IOmAJA0NDXlRBSJd5FBmkg8oVp4z6m2SPm1mt0i6QtIKM/uqu99WbmmQ4k/y8UFEQPlaBrW775W0V5KSM+o7CelyRA9lJvmAatBHXREm+QDk1VZQu/tRSUdLqaSHEcoAFoIz6oIxyQegaAR1QYZHG7rjkbGqy3gHoQz0DoK6AFVf2mCSD+htBPUCdTqkCWWgfgjqBRgebZQa0kzyAZAI6gXZ/ehYYa9FKAPIQlDP04ceOKK353GjPJN8ANpFUM/D9v1H9cr5n7Xcb8XSPj19/80dqAhAL1tUdQHdJu+KJxvXLCOkARSCoG7DrQeP5Q5pPqgIQFEI6pyGRxu57zgkpAEUiaDOKW+Hx5d/64ZS6wBQPwR1DpvvfSJXh8dtW9fRzQGgcAR1C5vvfUKvv3mx5X70QQMoC0E9hw89cCRXSG/bsIqQBlAagjpD3l7pjWuW8dkbAEpFUKfI2yt91fIldHgAKB1BPcs9w+O5QnrF0j4dv3t7ByoCUHcE9Sx5Pw2Puw4BdApBPc3me5/ItR+90gA6iaBOtNOGR680gE4iqEWvNIDYah/U9EoDiK7WQX3rwWO5eqWvWr6EXmkAlaltUN8zPJ7r0/BowwNQtVoGdd5FaVmhBUAEtQzqP3hkLNd+hDSACGoX1JvvfUJ51qSlVxpAFLUKanqlAXSj2gQ1bXgAulXLoDazK8zsKTM7aWanzOz+ThRWJNrwAHSzxTn2eVPSx9z9DTPrl/R9M/uOu/+w5NoKQRsegG7XMqjd3SW9kfyxP3nkmY+r3D3D47ThAeh6ua5Rm1mfmY1JOivpiLsfT9lnl5mNmNnI5ORkwWW2j15pAL0iV1C7+0V3v0HStZJuMrMPpOxzwN2H3H1ocHCw4DLbR680gF7RVteHu5+TdFRS6HTL2yt929Z1pdcCAAuVp+tj0MxWJtsDkj4u6dmS65q3vG14G9csow0PQFfI0/VxtaSHzKxPzWB/1N3/tdyy5qedlcNZlBZAt8jT9fG0pC0dqGVBbj14jJXDAfSknrgzkV5pAL2s64OaXmkAva6rgzpvr7REGx6A7tXVQb370bFc+/GRpQC6WdcG9fb9R/V2jmbpbRtW8ZGlALpaVwb19v1Hc3V4bFyzjE/DA9D1ui6o2wlp2vAA9IKuCmp6pQHUUdcEdd5eaUn0SgPoKV0R1Hl7pSU6PAD0nvBB3U6vNIvSAuhF4YP6zm+ezLXfbVvX8Wl4AHpS6KC+Z3hcb+VolmblcAC9LHRQf+1460se9EoD6HWhg7rVyTS90gDqIHRQz4WQBlAXoYO6f47qCGkAdRE6qJct7U8dXzmQPg4AvSh0UJ+bupA6/j8Z4wDQi8IG9fBoQ5bx3DUrBzpaCwBUKWxQ7zs8obSmD5O0Z8emTpcDAJUJG9SNc1Op4y5xmziAWgkb1H2WfuEjaxwAelXYoL7o6Xe7ZI0DQK8KG9TvfU96C17WOAD0qrBB/eaFi6njnFADqJuQQT082tD/XXg79Tl6qAHUTcig3nd4IvM5eqgB1E3IoD6T0Zon0UMNoH5CBvXKjAnD9/QvoocaQO2EDOqsCcMli/s6WwgABNAyqM3sOjP7npmdNrNTZnZ72UVlTRgykQigjvKcUb8l6Q/d/RckbZX0BTO7vsyisi59ZI0DQC9rGdTu/rK7/zjZPi/ptKRSLxTTQw0A72rrGrWZrZe0RdLxlOd2mdmImY1MTk7OuyB6qAFgptxBbWZXSvqWpDvc/fXZz7v7AXcfcvehwcHBeRdEDzUAzJQrqM2sX82QftjdD5VZED3UADBTnq4Pk/QVSafdfX/ZBdFDDQAz5Tmj3ibpdyV9zMzGksctZRVEDzUAzLS41Q7u/n0pc/nCwmUtaJs1DgC9LtydiazsAgAzhQtqVnYBgJnCBTUruwDATOGCOuvEmRNqAHUVLqizJg25KxFAXYUK6uHRRmZ7CXclAqirUEG97/CE0q5wmLgrEUB9hQrqRsbt4y5xVyKA2goV1PRQA8DlQgU1PdQAcLlQQU0PNQBcLlRQs7ILAFwuTFCzsgsApAsT1KzsAgDpwgQ1K7sAQLowQc3KLgCQLkxQs7ILAKQLE9RZE4ZMJAKouzBBnXXpI2scAOoiTFDTQw0A6UIENT3UAJAtRFDTQw0A2UIENT3UAJAtRFBnnTWvHOinhxpA7YUI6j07Nmmgf2a/9EB/n+779C9WVBEAxLG46gKkd1dv2Xd4QmfOTemalQPas2MTZ9MAoCBBLTXDmmAGgMuFuPQBAMhGUANAcAQ1AATXMqjN7EEzO2tmz3SiIADATHnOqP9R0s0l1wEAyNAyqN39SUmvdaAWAECKwq5Rm9kuMxsxs5HJycmiXhYAaq+woHb3A+4+5O5Dg4ODRb0sANQeXR8AEFyYOxOHRxvcQg4AKfK0531d0jFJm8zsJTP7fNFFDI82tPfQuBrnpuSSGuemtPfQuIZHG0W/FQB0nZZn1O7+22UXse/whKZmLcU1deGi9h2e4KwaQO2FuEadtXDAXAsKAEBdhAjqrIUDWIYLAIIEddbCASzDBQBBuj5YOAAAsoUIaomFAwAgS4hLHwCAbAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AARHUANAcAQ1AAS3OM9OZnazpL+S1CfpH9z9z4ouZHi0oX2HJ3Tm3JSuWTmgPTs2aeeWtUW/DQB0nZZBbWZ9kv5W0nZJL0n6kZk97u7/XlQRw6MN7T00rqkLFyVJjXNT2ntoXJIIawC1l+fSx02Snnf3n7j7zyR9Q9Jniixi3+GJd0L6kqkLF7Xv8ESRbwMAXSlPUK+V9N/T/vxSMjaDme0ysxEzG5mcnGyriDPnptoaB4A6yRPUljLmlw24H3D3IXcfGhwcbKuIa1YOtDUOAHWSJ6hfknTdtD9fK+lMkUXs2bFJA/19M8YG+vu0Z8emIt8GALpSnq6PH0naaGbvk9SQ9FlJv1NkEZcmDOn6AIDLtQxqd3/LzL4o6bCa7XkPuvupogvZuWUtwQwAKXL1Ubv7tyV9u+RaAAApuDMRAIIjqAEgOIIaAIIjqAEgOHO/7N6Vhb+o2aSkn87z21dLerXAcspAjcWgxmJQYzGqrvHn3T31bsFSgnohzGzE3YeqrmMu1FgMaiwGNRYjco1c+gCA4AhqAAguYlAfqLqAHKixGNRYDGosRtgaw12jBgDMFPGMGgAwDUENAMGFCWozu9nMJszseTO7q+JaXjSzcTMbM7ORZGyVmR0xs+eSr++dtv/epO4JM9tRUk0PmtlZM3tm2ljbNZnZLyd/t+fN7K/NLG1hiCJrvM/MGsmxHDOzWyqu8Toz+56ZnTazU2Z2ezIe5ljOUWOYY2lmV5jZU2Z2Mqnx/mQ80nHMqjHMcczN3St/qPnxqS9Ier+kJZJOSrq+wnpelLR61tifS7or2b5L0peS7euTepdKel/y9+groaaPSLpR0jMLqUnSU5I+rObKPd+R9MmSa7xP0p0p+1ZV49WSbky2l0v6j6SWMMdyjhrDHMvk9a5MtvslHZe0NdhxzKoxzHHM+4hyRl36AroF+Iykh5LthyTtnDb+DXd/093/U9Lzav59CuXuT0p6bSE1mdnVkla4+zFv/uv7p2nfU1aNWaqq8WV3/3GyfV7SaTXXAA1zLOeoMUsVNbq7v5H8sT95uGIdx6was1TybzKPKEGdawHdDnJJ3zWzE2a2Kxm7yt1flpo/SJLWJONV1t5uTWuT7dnjZfuimT2dXBq59Ktw5TWa2XpJW9Q80wp5LGfVKAU6lmbWZ2Zjks5KOuLu4Y5jRo1SoOOYR5SgzrWAbgdtc/cbJX1S0hfM7CNz7Butdim7pipq/TtJGyTdIOllSX+ZjFdao5ldKelbku5w99fn2jWjntLrTKkx1LF094vufoOa66jeZGYfmGP3SDWGOo55RAnq0hfQbYe7n0m+npX0mJqXMl5JfgVS8vVssnuVtbdb00vJ9uzx0rj7K8kPy9uSDurdy0KV1Whm/WoG4MPufigZDnUs02qMeCyTus5JOirpZgU7jmk1Rj2Oc4kS1O8soGtmS9RcQPfxKgoxs2VmtvzStqRPSHomqedzyW6fk/Qvyfbjkj5rZkutuQDwRjUnHjqhrZqSX0XPm9nWZNb696Z9Tyku/dAmfl3NY1lZjclrfkXSaXffP+2pMMcyq8ZIx9LMBs1sZbI9IOnjkp5VrOOYWmOk45hbJ2cu53pIukXN2e0XJN1dYR3vV3Pm96SkU5dqkfRzkv5N0nPJ11XTvufupO4JlTQbLOnrav6adkHN/+E/P5+aJA2p+Q/zBUl/o+Tu1BJr/GdJ45KeVvMH4eqKa/wVNX9tfVrSWPK4JdKxnKPGMMdS0mZJo0ktz0j64/n+nFRQY5jjmPfBLeQAEFyUSx8AgAwENQAER1ADQHAENQAER1ADQHAENQAER1ADQHD/D02xIB1C+UcKAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "y_inverse_quadric = interpolate.interp1d(y, x, kind=1) # normally it is x, y\n",
    "\n",
    "x_inverse_high_density = np.linspace(y.min(), y.max(), 1000)\n",
    "\n",
    "plt.scatter(x_inverse_high_density, y_inverse_quadric(x_inverse_high_density))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There exist also interpolations for 'cubic', 'nearest', etc. depending on the problem you like to solve. Normally you need to experiment a little, to find the one, which fit's the best. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}