How to draw a nice curve through given points ?

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
among all functions with . As is often the case, the length functional can be replaced with the elastic energy
because the piecewise linear (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 to shoot up over
. 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
and the function that minimizes it subject to conditions looks very nice indeed:

The Euler-Lagrange equation for the functional dictates that the fourth derivative of
is zero in the intervals between the knots
. Thus,
is a piecewise cubic polynomial. Also, both
and
must be continuous for any function with integrable second derivative. More delicate analysis is required for
, but it also can be shown to be continuous for minimizing function
; moreover,
must vanish at the endpoints
and
. 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 to the piecewise linear interpolant
. Here the spline is shown together with
(green) and
(magenta).

On each interval the correction term
is a cubic polynomial vanishing at both endpoints. The space of such polynomials is two-dimensional: thus,
on this interval. There are 20 coefficients ,
to find. At each of 9 knots 1,2,…9 we have two conditions:
must have removable singularity and
must jump by the amount opposite to the jump of
. Since
also vanishes at
, there are 20 linear equations for 20 unknowns.
It is easier to set up a linear system in terms of . Indeed, the values of
at two consecutive knots determine the correction term within:
and
. This leaves
equations (from the jumps of
) for
unknowns. The best part is that the matrix of this system is really nice: tridiagonal with dominant diagonal.
One can solve the system for 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*')