This is Words and Buttons Online — a collection of interactive #tutorials, #demos, and #quizzes about #mathematics, #algorithms and #programming.

Programmer's guide to polynomials and splines

Mathematical analysis explained with Python, blood, and TNT shows how you can analyze and synthesize arbitrary functions as polynomials. But you can synthesize polynomials from different inputs, too. You can put a spline through some points, make it sharp at one point and smooth at another, or make it flat near some point and explicitly curvy near another. Polynomials are your digital clay.

Polynomials or even polynomial splines might not always be the best tool for the job, but they possess some traits that programmers value. They are inherently simple, versatile, and, most noticeably, highly effective in terms of performance.

Let's say we have a polynomial written like this:

P(x) = ax³ + bx² + cx + d

This polynomial is defined for every x, it has a local minimum, a local maximum, and it covers all the numbers with its image too. It only takes 6 multiplications and 3 additions to calculate, and you can already describe non-trivial things with it. But even this can be reduced by using Horner's method. The same polynomial may be written as:

P(x) = ((ax+ b)x + c)x + d

And now it's only 3 multiplications and 3 additions. See, we barely started, and you already learned how to speed things up.

Polynomial interpolation

Fitting an n-degree polynomial in exactly n+1 points is called polynomial interpolation. There are several ways to do that. You can use Newton polynomials or Lagrange polynomials. But the very basic way you can get an interpolating polynomial is by solving a linear system.

If a polynomial P(x) goes through the point (x_i, y_i) then, obviously enough, we can claim that P(x_i) = y_i. Say we want to fit a polynomial into a set of three points. This means that:

P(x₁) = y₁
P(x₂) = y₂
P(x₃) = y₃

Generally, we can not fit a straight line into three arbitrary points. We would have to bend it forming a parabola. Or, put in other words, 2-nd degree polynomial also known as quadric.

P(x₁) = ax₁² + bx₁ + c = y₁
P(x₂) = ax₂² + bx₂ + c = y₂
P(x₃) = ax₃² + bx₃ + c = y₃

Since we know xs and ys, we only have to solve the system for (a, b, c) coefficients and since it's a three-equation system with three variables, it should generally give us one and only one solution.

Here, try it yourself. This is a canvas with three points on and a parabola that runs through them. You can drag them around to see what happens.

This interplation method is also a good mental model for understanding linear systems. Generally, you can't fit a straight line into three points just like you can't find a solution for (n-1) variables for n-equations system. But sometimes you can. It's when some of the points coincide or they all deliberately lie on a straight line.

And the vise versa, we can fit an infinite number of parabolas into a set of only two points. We just can't choose between them, that's why we can't unambiguously solve an n-equations system for (n+1) variables although there is an infinite number of solutions that fit.

But if the problem is in choosing the right function, maybe we can add some additional criteria to make the choise unambiguous?

Synthesis

And this brings us to polynomial synthesis. It is a crossover between the polynomial series and polynomial interpolation. With polynomial series, we can model a function based on its derivatives at some point, and with synthesis, we can use both points and derivatives. And even more than that actually, but that's a whole different story.

Function's derivative is closely related to the geometric properties of its plot. The first derivative sets a tangent, and the second one sets a curvature.

Let than say we want a function that goes through a pair of points and has a specific tangent in both points. We can easily synthesize this function as a polynomial.

Just like before, we should prepare a system of equations. Now we want four conditions, so we should pick a 3-rd degree polynomial — a cubic polynomial.

P(x₁)' = 3ax₁² + 2bx₁ + c = dy₁/dx
P(x₁) = ax₁³ + bx₁² + cx₁ + d = y₁
P(x₂) = ax₂³ + bx₂² + cx₂ + d = y₂
P(x₂)' = 3ax₂² + 2bx₂ + c = dy₂/dx

You see, some of the equations are formed from points and some from derivatives. You can even add integrals into the mix if you want to constraint some integral properties. This is a rather powerful technique.

Our four equations form a smooth continuous function between two points with tangential constraints at these points.

Now you have a nice obedient function between two points. But what if you have to fit a polynomial into many more points?

Runge's phenomenon

One unpleasant property of polynomial interpolation is that the function tends to oscillate at both ends of the interval more and more with adding more points. This is called Runge's phenomenon.

The other thing is that the interpolating polynomials are global, meaning that touching a single point causes the whole function to change. Now combine this with the oscillations, and you have yourself a recipe for chaos.

Chebyshev nodes

One way to mitigate this chaos would be to select a special grid for an interpolation: Chebyshev nodes. These are specially calculated x values that represent a projection of equidistant intervals from a semicircle with radius 1 onto the x-axis.

There is some mathematical magic behind this; pragmatically this disposition minimizes the Runge's phenomenon. In the range (-1; 1) the interpolation now works more or less predictable.

Of course, you can translate this range to any other range on x-axis. It doesn't have to be always (-1; 1).

But the interpolation is still global. The change in the first point still affects the function near the last one, although not so drastically.

Splines

There are quite a lot of different types of splines, but they all share the common motivation. When, for some reason, global interpolation doesn't work for us, we can divide our domain into pieces, and we can then define separate interpolating functions for each of them. We only have to make sure they conjoin in their ends so the function appears continuous. If we guarantee continuity not only for our resulting piecewise function but for its first derivative as well, then the tangents of every piece would match, and the function will become smooth, too.

There is a definite classification for the types of splines. For instance, here is a 2-piece polynomial spline. It is cubic, meaning that the polynomial in every piece is of the third degree. It has 1-st derivative continuity since we make the pieces conjoint along with their tangents. It is non-uniform since the length of the pieces' intervals is variable. It is not natural since you can rule the derivatives at the end. And it is, of course, an interpolating spline since it goes exactly through the grid points we set up.

Conclusion

It is highly unlikely that you would ever have to implement your own interpolation. There are a lot of ready-made solutions out there, most of the time you would only have to pick the right tool for the job. It's not that hard, but the amount of unknown words and names like Runge or Chebyshev * (fun fact, Chebyshev's real name is Chebyshóv. The confusion comes from that weird Russian letter «ё» that, unlike in other languages, doesn't preserve the sound “e” but rather changes it to “o” and softens the previous consonant. Russians don't like this letter very much themselves and tend to omit it when possible. But they substitute it with «е» not «о» thus breeding confusion)← may be overwhelming.

The purpose of this guide is to provide a very basic understanding of polynomials and splines. It is not at all comprehensive, in fact, there are whole books being written on every small piece of this guide. I hope, however, that it gives enough introduction to practicing programmers so they could choose the best solution for a problem without reading all those books.