Rosio Pavoris a blog

You’ve had this problem before: you have a bunch of data points, and you want to interpolate between them.
For various reasons, higher order polynomial interpolation (where you try to find an nth-degree polynomial through n + 1 of your data points) can be a bad idea, so you decide that rather than using a simple equation, you’ll use a series of them to connect your data points. These equations are splines, and the simplest form of spline interpolation is just, well, connecting your data points directly:

That’s pretty ugly, though. Is there a way to achieve something like this:

instead?
Yes, obviously, and one of those ways is to use quadratic splines instead. Let’s use a simpler example, though. Suppose we only have four data points, (x₀, y₀) through (x₃, y₃):

The black dots are our actual data points, the red lines are our linear splines. What we’d actually like, though, is this:

Turns out that’s not that hard to do. As you can see, every spline is a quadratic equation, which obviously is of the form f(x) = ax² + bx + c. So each spline equation has three unknowns (a, b, and c), and there are three splines, for a total of nine unknowns (let’s call them a₁ through a₃ and so on).

Since two points are known for each spline equation, that gives us the following six equations:

To solve for nine unknowns, obviously we need nine equations. So what else do we know?

Well, the reason the linear spline interpolation looks like crap is because of the sharp breaks at the spline edges, so we would like our neighboring quadratic splines to have the same slope in the point that they share. In other words, if our spline equations are f, g, and h, we want the derivative of f to equal the derivative of g in the point (x₁, y₁), and we want the derivative of g to equal the derivative of h in the point (x₂, y₂).
The derivative is easy enough to find:

Filling in, this gives us two more equations:

Or equivalently:

Which brings our total to eight equations. We aren’t going to squeeze another legitimate equation out of this, so let’s just fill in one of the unknowns ourselves. If we make one of the as equal to 0, one of the quadratic splines becomes a linear spline, which is fine. Let’s take a₁ for simplicity’s sake.
This enables us to construct the following matrix:

The first three columns are the as, the next three the bs, the next three the cs, and the final column will hold the solutions after reduction.
Filled in and solved for our particular dataset:

Which gives us the following equations for our splines:

Obviously this is a lot of work, but it’s mechanical work that doesn’t require a lot of judgement. Which is why I’ve written this Python script to do it for you. Feed it data points and it’ll produce gnuplot code to plot your splines:

$ python qsi.py < data.txt 
plot 1.000000 <= x && x <= 3.000000 ? 0.000000 * x * x + 1.500000 * x + 1.500000 :\
     3.000000 <= x && x <= 5.000000 ? -1.250000 * x * x + 9.000000 * x + -9.750000 :\
     5.000000 <= x && x <= 9.000000 ? 1.125000 * x * x + -14.750000 * x + 49.625000 : 0/0 notitle

As you can tell, it’s not necessarily gorgeous, but it (probably) works, and it’s not like anyone has to see the code itself.
Format for the input file is as you’d expect: two numbers per line, first being x and second y, sorted. If gnuplot‘s output is jagged, increase the sampling (set sample 1000).
And if it doesn’t work, fix it.

Edit: In light of overwhelming demand, this is a script that interpolates using a higher-order polynomial, as mentioned above. Here’s how the approaches compare for our sample dataset:

This script will fail if you only have one datapoint and its x value is 0, but everything else should work.