Connecting dots naturally

How to draw a nice curve through given points ${(x_i,y_i)}$ ?

One way is to connect them with straight lines, creating a piecewise linear function:

This is the shortest graph of a function that interpolates the data. In other words, the piecewise linear function minimizes the integral

$\displaystyle \int_0^{10} \sqrt{1+(f'(x))^2}\,dx$

among all functions with ${f(x_i)=y_i}$ . As is often the case, the length functional can be replaced with the elastic energy

$\displaystyle \mathcal{E}(f) = \int_0^{10} (f'(x))^2\,dx$

because the piecewise linear ${f}$ (and only it) minimizes it too.

Of course, it is not natural for the connecting curve to take such sharp turns at the data points. One could try to fit a polynomial function to these points, which is guaranteed to be smooth. With 11 points we need a 10th degree polynomial. The result is disappointing:

It is not natural for a curve connecting the points with ${1\le y\le 9}$ to shoot up over ${y=40}$ . We want a connecting curve that does not wiggle more than necessary.

To reduce the wiggling and remove sharp turns at the same time, one can minimize the bending energy of the function, thinking of its graph as a thin metal rod. This energy is

$\displaystyle \mathcal{B}(f) = \int_0^{10} (f''(x))^2\,dx$

and the function that minimizes it subject to conditions ${f(x_i)=y_i}$ looks very nice indeed:

The Euler-Lagrange equation for the functional ${\mathcal{B}}$ dictates that the fourth derivative of ${f}$ is zero in the intervals between the knots ${x_i}$ . Thus, ${f}$ is a piecewise cubic polynomial. Also, both ${f}$ and ${f'}$ must be continuous for any function with integrable second derivative. More delicate analysis is required for ${f''}$ , but it also can be shown to be continuous for minimizing function ${f}$ ; moreover, ${f''}$ must vanish at the endpoints ${0}$ and ${10}$ . Taken together, these properties (all derived from the variational problem) complete the description of a natural cubic spline.

It remains to actually construct one. I prefer to think of this process as adding a correction term ${C}$ to the piecewise linear interpolant ${L}$ . Here the spline is shown together with ${L}$ (green) and ${C}$ (magenta).

PL interpolant, correction term, and their sum: the cubic spline — PL interpolant, correction term, and their sum: cubic spline

On each interval ${[x_i,x_{i+1}]}$ the correction term ${C}$ is a cubic polynomial vanishing at both endpoints. The space of such polynomials is two-dimensional: thus,

$\displaystyle C(x) = \alpha_i (x_{i+1}-x)^2(x-x_i) + \beta_i (x_{i+1}-x)(x-x_i)^2$

on this interval. There are 20 coefficients ${\alpha_i}$ , ${\beta_i}$ to find. At each of 9 knots 1,2,…9 we have two conditions: ${C''}$ must have removable singularity and ${C'}$ must jump by the amount opposite to the jump of ${L'}$ . Since ${C''}$ also vanishes at ${1,10}$ , there are 20 linear equations for 20 unknowns.

It is easier to set up a linear system in terms of ${z_i=C''(x_i)}$ . Indeed, the values of ${C''}$ at two consecutive knots determine the correction term within: ${\alpha_i= \dfrac{z_{i+1}+2z_i}{6} }$ and ${\beta_i = \dfrac{2z_{i+1}+z_i}{6}}$ . This leaves ${n-1}$ equations (from the jumps of ${C'}$ ) for ${n-1}$ unknowns. The best part is that the matrix of this system is really nice: tridiagonal with dominant diagonal.

$\displaystyle \frac{1}{6}\begin{pmatrix} 4 & 1 & 0 & 0 & \ldots & 0 & 0 \\ 1 & 4 & 1 & 0 & \ldots & 0 & 0 \\ 0 & 1 & 4 & 1 & \ldots & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \ddots & \ldots & \ldots \\ 0 & 0 & 0 & 0 & \ldots & 1 & 4 \end{pmatrix}$

One can solve the system for ${z_i}$ within a for loop, but I used the Scilab solver instead. Here is the Scilab code for the most interesting part: the spline. The jumps of the derivative of the piecewise linear interpolant are obtained from the second order difference of the sequence of y-values.

a = 0; b = 10
y = [3 1 4 1 5 9 2 6 5 3 5]
n = length(y)-1
h = (b-a)/n
jumps = diff(y,2)/h
A = (h/6)*(diag(4*ones(1,n-1))+diag(ones(1,n-2),1)+diag(ones(1,n-2),-1))
z = [0,-jumps/A,0] 
allx = []; spl = []
for i=1:n  
   xL = a+h*(i-1)
   xR = a+h*i
   x = linspace(xL,xR,100)
   linear = y(i)*(xR-x)/h + y(i+1)*(x-xL)/h
   cor = ((z(i+1)+2*z(i))*(xR-x)+(2*z(i+1)+z(i))*(x-xL)).*(xR-x).*(x-xL)/(6*h)
   allx = [allx, x]   
   spl = [spl, linear + cor] 
end
plot(allx, spl) 
plot(a:h:b, y, 'r*')

Connecting dots naturally

Related

Leave a Reply Cancel reply