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 p(t): # calculate the value P(t) for arbitrary t \in R
  sum = 0
  for i in range(n+1):
    prod = y[i]
    for j in range(n+1):
      if i!= j:
        prod=prod*(t-x[j])/(x[i]-x[j])
    sum = sum + prod
  return sum

plot.plot(prange,p(prange))
plot.show()

Newton polynomial

Although more complex, Newton polynomial is a somewhat smarter version of the Vandermonde polynomial:

$$ P_{n}(x)=c_{0}+c_{1}\left(x-x_{0}\right)+c_{2}\left(x-x_{0}\right)\left(x-x_{1}\right)+\cdots+c_{n}\left(x-x_{0}\right) \cdots\left(x-x_{n-1}\right) \nonumber $$

Which gives a polynomial of degree $n$. In fact, if we factor all the $(x-x_k)$ terms out, this is exactly the same polynomial that we would find via Vandermonde's method, but written differently.

One reason that Newton's form is better is that an updated $P_{n+1}(x)$ interpolating all the original points plus a further point $(x_{n+1},y_{n+1})$ can be constructed effortlessly:

$$ P_{n+1}(x)=P_n(x)+ c_{n+1}\left(x-x_{0}\right) \cdots\left(x-x_{n}\right) $$

Notice $\forall i \in [0,n], i \in N: c_{n+1}\left(x-x_{0}\right) \cdots\left(x-x_{n}\right)=0$ which means $\forall i \in [0,n], i \in N: P_{n+1}(x_i)=P_n(x_i)$. So $P_{n+1}(x)$ indeed interpolates for all the original datapoints.

Numerical solution

$$ \begin{aligned} y_{0} &=c_{0} \\ y_{1} &=c_{0}+c_{1}\left(x_{1}-x_{0}\right) \\ y_{2} &=c_{0}+c_{1}\left(x_{2}-x_{0}\right)+c_{2}\left(x_{2}-x_{0}\right)\left(x_{2}-x_{1}\right) \\ \vdots & \vdots \\ y_{n} &=c_{0}+c_{1}\left(x_{n}-x_{0}\right)+c_{2}\left(x_{n}-x_{0}\right)\left(x_{n}-x_{1}\right)+\cdots+c_{n}\left(x_{n}-x_{0}\right) \cdots\left(x_{n}-x_{n-1}\right) \end{aligned} \nonumber $$

Written in matrix form:

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

Notice that this is a lower-triangle matrice, meaning it is easier to solve than the Vandermonde polynomial. A second advantage of Newton's form.

Obviously, in practice, with modern computers, the difference is a few microseconds.

Python implementation

import numpy as np

def newton_interpolation(x: np.ndarray, y: np.ndarray) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    y = np.asarray(y, dtype=float)

    if x.ndim != 1 or y.ndim != 1:
        raise ValueError("x and y must be 1-D arrays")

    if len(x) != len(y):
        raise ValueError("the length of x and y arrays must be equal")

    if len(x) == 0:
        raise ValueError("x and y must not be empty")

    if len(np.unique(x)) != len(x):
        raise ValueError("x values must be distinct")

    m = len(x)
    A = np.zeros((m, m), dtype=float)

    for i in range(m):
        A[i, 0] = 1.0
        for j in range(1, m):
            A[i, j] = (x[i] - x[j - 1]) * A[i, j - 1]

    return np.linalg.solve(A, y)

And to calculate $P(t)$ for any $t \in \mathbb{R}$:

def newton_eval(x_nodes: np.ndarray, coeffs: np.ndarray, t: np.ndarray) -> np.ndarray:
    x_nodes = np.asarray(x_nodes, dtype=float)
    coeffs = np.asarray(coeffs, dtype=float)
    t = np.asarray(t, dtype=float)

    p_t = np.full_like(t, coeffs[-1], dtype=float)
    for k in range(len(coeffs) - 2, -1, -1):
        p_t = coeffs[k] + (t - x_nodes[k]) * p
    return p_t

Complete example with plotting of $(x_i,y_i)$ and $P(x)$:

# Add the import line to the top of the file
import matplotlib.pyplot as plt

# --- Data ---
x_data = np.array([0.0, 1.0, 2.0, 3.0, 4.0])
y_data = np.array([1.0, 2.5, 0.5, 3.0, 1.5])

# --- Compute Newton coefficients ---
coeffs = newton_interpolation(x_data, y_data)

# --- Evaluate on a fine grid ---
x_plot = np.linspace(x_data.min() - 0.5, x_data.max() + 0.5, 500)
y_plot = newton_eval(x_data, coeffs, x_plot)

# --- Plot ---
plt.figure(figsize=(8, 5))
plt.plot(x_plot, y_plot, 'b-', linewidth=2, label='Interpolating polynomial P(x)')
plt.plot(x_data, y_data, 'ro', markersize=8, label='Interpolation points')

# sanity check: mark the interpolated values and make sure they overlap visually with the original datapoints

plt.plot(x_data, y_data, 'ro', markersize=8, label='Interpolation points')

plt.xlabel('x')
plt.ylabel('y')
plt.title('Newton Interpolation')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.show()

Algebraic solution

Start with $i=0$,

$$\begin{gathered} y_0=P_0(x_0) \\ P_0(x_0)=c_0 \\ \Rightarrow c_0=y_0 \end{gathered}$$

And use the formula:

$$ \begin{aligned} P_{i+1}(x_{i+1}) &= P_i(x_{i+1}) + c_{i+1}\left(x_{i+1}-x_0\right)\cdots\left(x_{i+1}-x_i\right) \\[2mm] \Rightarrow\quad c_{i+1} &= \frac{y_{i+1}-P_i(x_{i+1})} {\left(x_{i+1}-x_0\right)\cdots\left(x_{i+1}-x_i\right)} \end{aligned} $$

to recursively workout $P_n(x)$