Polynomial interpolation

Polynomial interpolation is the process of finding a polynomial

$$P(x)=\sum_{k=0}^{n} a_k x^k, n\in\mathbb{N},a_0, a_1, \dots, a_n\in \mathbb{F}$$

where $x\in \mathbb{F}$ is an independent variable

that for a definite set of datapoints $\{(x_i,y_i)| i\in [0,n], i \in N\}$ satisfies $P(x_i)=y_i$.

A polynomial of degree $n$ has $n+1$ arbitrary constants $a_k$ and therefore $n+1$ degrees of freedom, corresponding to $n+1$ datapoints $\{(x_i,y_i)| i\in [0,n], i \in N\}$.

It can be therefore stipulated $\exists P(x), \forall i, P(x_i)=y_i$.

Methods

Vandermonde polynomial

$$ P(x)=c_{0} x^{n}+c_{1} x^{n-1}+\cdots+c_{n} . \nonumber $$

Since $P(x_i)=y_1$, a linear system of $n$ equations and $n$N variables ($c_n$) can be constructed:

$$ \begin{gathered} y_{0}=c_{0} x_{0}^{n}+c_{1} x_{0}^{n-1}+\cdots+c_{n-1} x_{0}+c_{n} \\ y_{2}=c_{0} x_{1}^{n}+c_{1} x_{1}^{n-1}+\cdots+c_{n-1} x_{1}+c_{n} \\ \vdots \\ y_{n}=c_{0} x_{n}^{n}+c_{1} x_{n}^{n-1}+\cdots+c_{n-1} x_{n}+c_{n} \end{gathered} \nonumber $$

Which can also be written in matrix form:

$$ \left(\begin{array}{ccccc} x_{0}^{n} & x_{0}^{n-1} & \cdots & x_{0} & 1 \\ x_{1}^{n} & x_{1}^{n-1} & \cdots & x_{1} & 1 \\ \vdots & \vdots & \ddots & \vdots & \\ x_{n}^{n} & x_{n}^{n-1} & \cdots & x_{n} & 1 \end{array}\right)\left(\begin{array}{c} c_{0} \\ c_{1} \\ \vdots \\ c_{n} \end{array}\right)=\left(\begin{array}{c} y_{0} \\ y_{1} \\ \vdots \\ y_{n} \end{array}\right) \nonumber $$

Solving the system gives $P(x)$

Lagrange polynomial

$$ P(x)=\sum_{j=0}^n\left(y_j\prod_{\substack{k=0 \\ k\neq j}}^n\frac{x-x_k}{x_j-x_k}\right) $$

Since $\prod_{\substack{k=0 \\ k\neq j}}^n\frac{x-x_k}{x_j-x_k}$ always produce a polynomial of leading term $x^n$ (degree $n$), P(x), which is a sum of those polynomials, is itself a polynomial of degree $n$.

It can be further verified $P(x_i)=y_i$ by substituting $x=x_i$ into the original expression.

Intuition: Over $(x_i,y_i)$, Lagrange polynomial has the same effects as $P(x_i)=\sum_{j=0}^m \delta_{ij} y_i$

Where $\delta_{ij}$ is the kronecker delta function:

$$ \delta_{ij} = \begin{cases} 1, & i=j,\\ 0, & i\ne j. \end{cases} $$

Python implementation

import numpy as np
import matplotlib.pyplot as plot

x=[1,3,5,7,8,9,10,12,13]
y=[50,-30,-20,20,5,1,30,80,-10]

n=len(x)-1
prange = np.linspace(min(x),max(x),500)

plot.plot(x,y,marker='o', color='r', ls='', markersize=10)

def f(o):
  sum = 0
  for i in range(n+1):
    prod = y[i]
    for j in range(n+1):
      if i!= j:
        prod=prod*(o-x[j])/(x[i]-x[j])
    sum = sum + prod
  return sum

plot.plot(prange,f(prange))
plot.show()