02 Polynomial Interpolation

The Core Problem: What Happens Between the Measurements?

Let’s begin with a simple, practical scenario. Imagine you have a set of precise measurements. These could be the altitude of a drone at specific times, the pressure inside an engine cylinder at different crank angles, or the value of a stock at the close of each day. You have a set of points:

$$(x_0,y_0),(x_1,y_1),\dots ,(x_n,y_n)$$

You know the value at these points. But the crucial question is: what is the value between these points? If you plot these points on a graph, your brain almost automatically connects them, imagining a smooth curve passing through them. The task of interpolation is to make this imagined curve mathematically concrete.

The Interpolation Task: A Formal Statement

Given a set of $n+1$ data points $(x_i,y_i)$ for $i=0,\dots ,n$, where all the $x_i$ are distinct, the goal is to find a function $f(x)$ that satisfies the interpolation conditions:

$$f(x_i)=y_i\text{for all }i=0,\dots ,n$$

We are assuming our measurements are exact for now. The case of noisy data, where the curve shouldn’t necessarily pass through every point, is a problem of regression or approximation, which we’ll discuss later.

The First Idea: Just Connect the Dots

The most straightforward way to create a curve is to draw straight lines between consecutive points. This is called piecewise linear interpolation. It’s simple and intuitive.

However, as was pointed out in the lecture, this approach has a significant drawback: the resulting function has sharp corners at each data point. It is continuous, but it is not differentiable at the points $x_i$. For many physical systems, like the trajectory of a drone, we expect a smooth path without instantaneous changes in direction. We need a smoother solution.

The Framework: Choosing Our Reality from Function Spaces

The problem is that there are infinitely many possible smooth curves that can pass through a given set of points. How do we choose just one? We need a systematic approach.

This is where the idea of function spaces comes in. We assume that our discrete measurements are just samples from some “true” but unknown underlying function, let’s call it $\tilde{f}(x)$. This function lives in some vast, infinite-dimensional space of functions, $V$. For example, $V$ could be the space of all continuous functions, $C([a,b])$.

Since a computer cannot work with an infinite-dimensional object, we must simplify. We choose a much smaller, finite-dimensional subspace, let’s call it $V_n$, and search for our solution only within that subspace. This is a modeling choice. We are saying, “I don’t know what the true function is, but I’m going to assume it can be well-approximated by a function from my chosen simple family, $V_n$.”

Basis Functions: The Building Blocks of Our Model

How do we define this simpler space $V_n$? We define it with a basis. Just as the vectors $\mathbf{i},\mathbf{j},\mathbf{k}$ form a basis for 3D space, we can choose a set of basis functions $\{b_0(x),b_1(x),\dots ,b_n(x)\}$ to span our function space $V_n$.

Any function $f_n(x)$ in our chosen space can then be written as a unique linear combination of these basis functions:

$$f_n(x)={\alpha }_0{b}_0(x)+{\alpha }_1{b}_1(x)+\cdots +{\alpha }_n{b}_n(x)=\sum _{j=0}^n{\alpha }_j{b}_j(x)$$

The interpolation problem is now transformed from a vague search for a “function” into a concrete, algebraic problem: find the specific coefficients ${\alpha }_j$ that make our model $f_n(x)$ satisfy the interpolation conditions.

The choice of basis is the most critical decision we will make. It defines the character of our solution and the difficulty of finding it.

• Polynomials: A natural choice for smooth functions.
• Trigonometric Functions (Sines and Cosines): The ideal choice for periodic phenomena. This leads to Fourier analysis and the Fast Fourier Transform (FFT), one of the most important algorithms ever developed.
• Splines: Piecewise polynomials that are smoothly connected. These are used everywhere in computer graphics and engineering design.

For now, we will focus on the simplest and most fundamental choice: polynomials.

Polynomial Interpolation: The Naive Approach (and Why It Fails)

Let’s try to find a single polynomial that passes through all $n+1$ of our data points. To satisfy $n+1$ conditions, we need $n+1$ degrees of freedom. This leads us to choose a polynomial of degree at most $n$.

The Monomial Basis

The most obvious basis for the space of polynomials of degree $n$ is the monomial basis: $\{1,x,x^2,\dots ,x^n\}$. Our interpolating polynomial is then:

$$p_n(x)={\alpha }_0+{\alpha }_{1}x+{\alpha }_2{x}^2+\cdots +{\alpha }_n{x}^n$$

To find the coefficients ${\alpha }_j$, we enforce the interpolation conditions $p_n(x_i)=y_i$ for each point. This gives us a system of $n+1$ linear equations.

This system, $V\alpha =y$, involves the Vandermonde matrix. In theory, this matrix is invertible as long as the points $x_i$ are distinct, which means a unique solution for the coefficients ${\alpha }_j$ exists.

What is a Vandermonde Matrix?

A Vandermonde matrix is a very structured matrix that shows up whenever you express a polynomial in the monomial basis and enforce interpolation conditions. For interpolation points $x_0,x_1,\dots ,x_n$, it looks like this:

$$V=\left[\begin{array}{ccccc}1& x_0& x_0^2& \dots & x_0^n\\ 1& x_1& x_1^2& \dots & x_1^n\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_n& x_n^2& \dots & x_n^n\end{array}\right]$$

Each row corresponds to plugging one $x_i$ into the polynomial, and each column corresponds to one power of $x$.

Why does it show up here? Because when we write

$$p_n(x)={\alpha }_0+{\alpha }_{1}x+{\alpha }_2{x}^2+\cdots +{\alpha }_n{x}^n$$

and demand $p_n(x_i)=y_i$, we are literally saying:

$${\alpha }_0+{\alpha }_1{x}_i+{\alpha }_2{x}_i^2+\cdots +{\alpha }_n{x}_i^n=y_i$$

for each $i$. If you stack these equations for all $i$, you get the matrix-vector system $V\alpha =y$.

Why are Vandermonde matrices “bad”? Two main reasons:

• Ill-conditioning: The determinant of $V$ involves products of differences $(x_j-x_i)$. If the $x_i$ are close to each other or $n$ is large, the determinant becomes very small, meaning the matrix is nearly singular. That makes numerical inversion unstable: small errors in input create huge errors in output.
• Growth of powers: In each row, you have powers like $x_i^n$. For moderately large $n$, these numbers can explode or vanish depending on the size of $x_i$. This amplifies rounding errors when working in floating-point arithmetic.

In short: Vandermonde matrices are the unavoidable result of using the monomial basis in interpolation, but they are numerically treacherous.

Problem solved? Not at all. This approach is a numerical disaster.

Why the Monomial Basis is a Trap

1. Computational Cost: Solving the dense Vandermonde system with Gaussian elimination costs $O(n^3)$ operations. For many points, this is far too slow.
2. Numerical Instability: The Vandermonde matrix is famously ill-conditioned. This means it is “almost singular.” Tiny rounding errors in your input data get magnified into enormous errors in the computed coefficients. As the lecturer stated, for even a moderate number of points, the results are “brutal schlecht” (brutally bad) and “völlig sinnlos” (completely senseless). Standard numerical software will even warn you not to trust the solution.

Even if we could find the coefficients exactly, evaluating the polynomial in its standard form is also prone to numerical issues. A much more stable and efficient method for evaluation is Horner’s Schema, which uses a nested form:

$$p_n(x)={\alpha }_0+x({\alpha }_1+x({\alpha }_2+\cdots +x({\alpha }_{n-1}+x{\alpha }_n)\dots ))$$

This reduces the number of multiplications and is less susceptible to overflow and rounding errors.

A Better Way: The Newton Basis

The failure of the monomial basis does not mean polynomial interpolation is a bad idea. It means we chose the wrong building blocks. We need a smarter basis.

The Newton form of the interpolating polynomial is built on a beautifully simple, incremental idea. Imagine you are receiving your data points one by one.

1. First point $(x_0,y_0)$: The simplest polynomial that passes through this point is a constant function, a polynomial of degree 0: $p_0(x)=y_0$
2. Second point $(x_1,y_1)$: We want to find $p_1(x)$ that passes through both points. We can write it as a correction to our previous polynomial: $p_1(x)=p_0(x)+\text{correction term}$ The correction term must be zero at $x_0$ so we don’t mess up our first condition. The simplest way to ensure this is to include a factor of $(x-x_0)$. So, we try: $p_1(x)=p_0(x)+{\beta }_1(x-x_0)$ We find the constant ${\beta }_1$ by enforcing $p_1(x_1)=y_1$.
3. Third point $(x_2,y_2)$: We update again. The correction term must now be zero at both $x_0$ and $x_1$. The simplest way is to include factors of $(x-x_0)$ and $(x-x_1)$: $p_2(x)=p_1(x)+{\beta }_2(x-x_0)(x-x_1)$ We find ${\beta }_2$ by enforcing $p_2(x_2)=y_2$.

This incremental process defines the Newton basis functions:

$$\begin{array}{rl}{N}_0(x)& =1\\ N_1(x)& =(x-x_0)\\ N_2(x)& =(x-x_0)(x-x_1)\\ & ⋮\\ N_k(x)& =\prod _{i=0}^{k-1}(x-x_i)\end{array}$$

The interpolating polynomial is then written in the Newton form:

$$p_n(x)={\beta }_0{N}_0(x)+{\beta }_1{N}_1(x)+\cdots +{\beta }_n{N}_n(x)$$

The coefficients ${\beta }_j$ are called divided differences.

The divided differences are defined recursively and represent the leading coefficient of the polynomial interpolating a subset of the points.

• 0-th order: $y[x_i]=y_i$
• k-th order: $y[x_i,\dots ,x_{i+k}]=\frac{y[x_{i+1},\dots ,x_{i+k}]-y[x_i,\dots ,x_{i+k-1}]}{x_{i+k}-x_i}$

The coefficients in the Newton form are the divided differences along the top diagonal of the computation table: ${\beta }_k=y[x_0,x_1,\dots ,x_k]$.

Good video to grasp the concept…

Why the Newton Basis is Superior

This approach elegantly solves the problems of the monomial basis.

1. Efficient and Stable Coefficient Calculation: If we write out the interpolation conditions using the Newton basis, the resulting system matrix is lower triangular. A triangular system is trivial to solve. We can find the coefficients using forward substitution in only $O(n^2)$ operations, which is much faster than $O(n^3)$ and, more importantly, is numerically stable.
2. Adaptive: If a new data point $(x_{n+1},y_{n+1})$ arrives, we don’t have to start from scratch. We simply keep our existing polynomial and add one new term, ${\beta }_{n+1}{N}_{n+1}(x)$. We only need to compute one new row in the divided difference table. This is impossible with the monomial basis, where every single coefficient would change.

Summary: Monomial vs. Newton Basis

Feature	Monomial Basis ($x^j$)	Newton Basis ($\prod (x-x_i)$)
System Matrix	Vandermonde (Dense, Ill-conditioned)	Lower Triangular (Well-conditioned)
Cost to find coeffs.	$O(n^3)$	$O(n^2)$
Adding a new point	Requires complete re-computation	Easy and efficient update
Practical Use	Avoid! For theoretical use only.	Recommended. Stable and efficient.

The Newton form of the interpolating polynomial provides a practical, stable, and efficient method for polynomial interpolation, addressing the severe shortcomings of the naive monomial approach.

Continue here: 03 Polynomial Interpolation using Lagrange Form, Barycentric Lagrange Form, Runge’s Phenomenon, Chebyshev Nodes