polynomial_interpolation_equations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derive the formula for a polynomial that passes through three equally spaced points\n",
    "\n",
    "We assume the polynomial passes through $(-1, y_1)$, $(0, y_2)$, $(1, y_3)$ and $(2, y_4)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A third degree polynomial can be written as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "a \\cdot x^3 + b \\cdot x^2 + c \\cdot x + d = y\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we know this equation to hold for the four points above, we get four equations that will determine the four coefficients $a$, $b$, $c$ and $d$ that we need. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "a \\cdot (-1)^3 + b\\cdot(-1)^2 + c\\cdot -1 + d &= y_1\\\\\n",
    "d &= y_2\\\\\n",
    "a \\cdot 1^3 + b\\cdot 1 ^2 + c\\cdot 1 + d &= y_3\\\\\n",
    "a \\cdot 2^3 + b\\cdot 2^2 + c\\cdot 2 + d &= y_4\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The formula for $d$ is the easiest. We can immediately combine this with the first equation to get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "-a + b - c + d &= y_1\\\\\n",
    "-a + b - c &= y_1 - y_2\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now eliminate $d$ in the same way from the third equation, and add the result to the formula above for $-a+b-c$. This leads to an expression for $b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "a + b + c &= y_3 - y_2\\\\\n",
    "2\\cdot b &= y_1 - 2\\cdot y_2 + y_3\\\\\n",
    "b &= \\frac{1}{2}\\cdot y_1 - y_2 + \\frac{1}{2}\\cdot y_3\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, starting with the fourth equation and subtracting twice the formula for $a + b + c$ above, we can isolate $a$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "8 a + 4b + 2c &= y_4 - y_2\\\\\n",
    "6a + 2b &= y_4 - y_2 - 2 y_3 + 2 y_2\\\\\n",
    "6a &= y_4 + y_2 - 2y_3 - 2\\cdot\\left(\\frac{1}{2}y_1 - y_2 + \\frac{1}{2}y_3\\right)\\\\\n",
    "6a &= y_4 + 3 y_2 - 3 y_3 - y_1\\\\\n",
    "a &= -\\frac{1}{6}y_1 + \\frac{1}{2}y_2 - \\frac{1}{2}y_3 +  \\frac{1}{6} y_4\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need $c$, but you notice, this is hard work. We can instead use matrix inversion to let the computer derive all of these equations for us. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-1,  1, -1,  1],\n",
       "       [ 0,  0,  0,  1],\n",
       "       [ 1,  1,  1,  1],\n",
       "       [ 8,  4,  2,  1]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.vstack([[x ** 3, x ** 2, x, 1] for x in [-1, 0, 1, 2]])\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.16666667,  0.5       , -0.5       ,  0.16666667],\n",
       "       [ 0.5       , -1.        ,  0.5       ,  0.        ],\n",
       "       [-0.33333333, -0.5       ,  1.        , -0.16666667],\n",
       "       [ 0.        ,  1.        ,  0.        ,  0.        ]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.inv(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the coefficients on the first row are exactly the coefficients we derived for $a$ in terms of $y_1$, $y_2$, $y_3$ and $y_4$. The second row is for $b$ and so forth. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some other ideas, for example if we instead know the values at -2, -1, 0, 1 and 2, and want a 4th degree polynomial, these are the coefficients. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.04166667, -0.16666667,  0.25      , -0.16666667,  0.04166667],\n",
       "       [-0.08333333,  0.16666667,  0.        , -0.16666667,  0.08333333],\n",
       "       [-0.04166667,  0.66666667, -1.25      ,  0.66666667, -0.04166667],\n",
       "       [ 0.08333333, -0.66666667, -0.        ,  0.66666667, -0.08333333],\n",
       "       [ 0.        ,  0.        ,  1.        ,  0.        ,  0.        ]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.vstack([[x ** i for i in [4, 3, 2, 1, 0]] for x in [-2, -1, 0, 1, 2]])\n",
    "np.linalg.inv(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or if we know the coefficients at 0, $\\frac{1}{3}$, $\\frac{2}{3}$ and 1,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -4.5,  13.5, -13.5,   4.5],\n",
       "       [  9. , -22.5,  18. ,  -4.5],\n",
       "       [ -5.5,   9. ,  -4.5,   1. ],\n",
       "       [  1. ,   0. ,   0. ,   0. ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = np.vstack([[x ** i for i in [3, 2, 1, 0]] for x in [0, 1/3, 2/3, 1]])\n",
    "np.linalg.inv(X)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}