Cubic Spline interpolation in C++

Aims

Download Source Code

It is implemented as a single header file:

Usage

#include <cstdio>
#include <cstdlib>
#include <vector>
#include "spline.h"

int main(int argc, char** argv) {

   std::vector<double> X(5), Y(5);
   X[0]=0.1; X[1]=0.4; X[2]=1.2; X[3]=1.8; X[4]=2.0;
   Y[0]=0.1; Y[1]=0.7; Y[2]=0.6; Y[3]=1.1; Y[4]=0.9;

   tk::spline s;
   s.set_points(X,Y);    // currently it is required that X is already sorted

   double x=1.5;

   printf("spline at %f is %f\n", x, s(x));

   return EXIT_SUCCESS;
}
Compile:
$ g++ -Wall demo.cpp -o demo
$ ./demo
spline at 1.500000 is 0.915345

Comparison

The result of this library is compared against alglib. Interpolation results are identical but extrapolation differs, as this library is designed to extrapolate linearly.
spline compare with alglib
Source:

Benchmark

A simple benchmark shows that performance is comparable with alglib. The following operations are measured: The output "cycl" is the number of cpu cycles required for a single operation.
$ ./bench 
usage: ./bench <cpu mhz> <spline size> <loops in thousand>

$ ./bench 2666 10 10000	    # splines based on 10 grid points, 10million loops
                                tk                      alglib
random access:   loops=1e+07,   0.198s (  53 cycl)      0.400s ( 107 cycl)
spline creation: loops=1e+06,   2.279s (6075 cycl)      3.566s (9507 cycl)
grid transform:  loops=1e+06,   2.685s (7157 cycl)      3.356s (8948 cycl)

accuracy: max difference = 5.55e-16, l2-norm difference = 2.46e-20

$ ./bench 2666 10000 10000  # splines based on 10000 grid points, 10million loops
                                tk                      alglib
random access:   loops=1e+07,   1.370s ( 365 cycl)      1.538s ( 410 cycl)
spline creation: loops=1e+03,   1.726s (4.6e+06 cycl)   1.691s (4.5e+06 cycl)
grid transform:  loops=1e+03,   2.258s (6.0e+06 cycl)   1.414s (3.8e+06 cycl)

accuracy: max difference = 4.41e-13, l2-norm difference = 7.85e-19
Source: Compile as follows where the alglib source and compiled lib needs to be in the appropriate directory:
$ g++ -Wall -O2 -I../alglib/src -c bench.cpp -o bench.o
$ g++ bench.o -o bench -L../alglib/ -lalglib

Maths

$ \newcommand{\R}{\mathbb{R}} \newcommand{\set}[1]{\left\{ #1 \right\}} $ Given a set of inputs $\set{(x_1,y_1),\dots,(x_n,y_n)}$ with $x_1 \lt x_2 \lt \dots \lt x_n$. Define piecewise cubic polynomials $f_1,\dots,f_{n-1}$ with $f_i :[x_i,x_{i+1}]\to\R$ by \begin{equation} f_i(x) := a_i (x-x_i)^3 + b_i (x-x_i)^2 + c_i (x-x_i) + y_i, \qquad x\in[x_i, x_{i+1}], \end{equation} with derivatives \begin{equation*} \begin{split} f'_i(x) & = 3 a_i (x-x_i)^2 + 2 b_i (x-x_i) + c_i ,\\ f''_i(x) & = 6 a_i (x-x_i) + 2 b_i .\\ \end{split} \end{equation*} In order for the overall function to be twice continuously differentiable we require \begin{equation*} \begin{aligned} f_i(x_{i+1}) & =y_{i+1}: & a_i h_i^3 + b_i h_i^2 + c_i h_i & = y_{i+1}-y_i,\\ f'_{i-1}(x_i) & =f'_i(x_i): & 3 a_{i-1} h_{i-1}^2 + 2 b_{i-1} h_{i-1} + c_{i-1} & = c_i,\\ f''_{i-1}(x_i) & =f''_i(x_i): & 6 a_{i-1} h_{i-1} + 2 b_{i-1} & = 2 b_i, \end{aligned} \end{equation*} with $h_i:=x_{i+1}-x_i$. This gives $3(n-1)$ equations which can be simplified by first expressing $a$ and $c$ in terms of $b$ and then solving the remaining equation system for $b$: from the continuity of the second derivative it follows \begin{align} f''_{i}(x_{i+1}) & =f''_{i+1}(x_{i+1}) & \Rightarrow & \quad a_i =\frac{b_{i+1}-b_i}{3 h_i}, \label{eq:spline_a} \\ f_i(x_{i+1}) & =y_{i+1} & \Rightarrow & \quad c_i = \frac{y_{i+1}-y_i}{h_i} - \frac{1}{3} (2 b_i + b_{i+1}) h_i, \label{eq:spline_c} \\ f'_{i-1}(x_i) & =f'_i(x_i) & \Rightarrow & \quad \frac{1}{3} h_{i-1} b_{i-1} + \frac{2}{3} (h_{i-1}+h_i) b_i + \frac{1}{3} h_i b_{i+1} = \frac{y_{i+1}-y_i}{h_i} - \frac{y_i-y_{i-1}}{h_{i-1}} . \label{eq:spline_b} \end{align} First solve the tridiagonal equation system \eqref{eq:spline_b} then using \eqref{eq:spline_a} and \eqref{eq:spline_c} to get the remaining coefficients. For extrapolation we use second order polynomials (as we have one condition less). From the continuity conditions it follows: \begin{equation*} \begin{split} f_0(x) & := b_1 (x-x_1)^2 + c_1 (x-x_1) + y_1, \quad x \leq x_1, \\ f_n(x) & := b_n (x-x_n)^2 + c_n (x-x_n) + y_n, \quad x \geq x_n. \end{split} \end{equation*} By defining $f_n$, the right boundary condition simplifies. For example for zero curvature at $x_1$ and $x_n$ we require $b_1=0$ and $b_n=0$.

Math symbols are rendered by MathJax which requires JavaScript.