Given n+1 values of a function y(x), obtained, for example, by experiments in a laboratory, find an interpolation function Y(x) which passes through the points (x,y) and permits evaluation of approximations to the function y(x) for arbitrary values of x. Naturally, one wants to use interpolation functions which are readily evaluated and differentiated, the last being important for later work on differential equations.
One's first thought is of polynomials, so that one will arrive at an expansion
where the coefficients Ai have to be evaluated from the system of n+1 simultaneous equations, linear in the Ai:
and it must at first the conditions be discovered under which this system has a solution, i.e., one must find out when the determinant
This determinant equals the product
as is shown by the algorithm which utilizes the fact that the value of a determinant is not changed by addition of any row or column multiplied by a constant:
Thus, as long as all the xi are different, i.e., the required approximation is not double-valued, the system has a unique solution. For example, it could be solved by Cramer's Rule which employs altogether n + 2 determinants.
Unfortunately, determinants must be handled with great care, because they involve differences of products which may lead to a loss of accuracy, as is shown by the example
Thus, it
is natural to look for an alternative method of finding
polynomial interpolating functions. Such a method was found 200
years ago by the French mathematician Lagrange, who wrote:
He derived these rational functions in the form
Fit a quadratic to ex in the range 0 < x < 1 using the points:
One has
Evaluate now Y(x) for x = 1/4 and x = 3/4:
Y(1/4) = 3/8*1 +
3/4*1.649 - (1/8)*2.718 = 0.375 + 1.237 - 0.340 = 1.272,
Y(3/4) = (-1/4)*(1/2)*1+ 4*(3/4)*(1/4)*1.649 +
(3/4)*(1/2)*2.718 = 2.131,
while the exponential function yields
y(1/4) = 1.284,
y(3/4) = 2.117.
Hence
Y(1/4) - y(1/4) =
- 0.012,
Y(3/4) - y(3/4) = 0.014,
i.e., Y(x) oscillates around y(x). If a higher order polynomial is used, the oscillations will become larger, which is a good reason not to go to very high order polynomials. As a rule, given a larger number of data points, it is better to subdivide them into groups of at most five points and to interpolate by means of fourth order polynomials.
Interpolate the function
using the data
one finds
Evaluate Y(x) and y(x):
Hence it is seen that in between data points the interpolation function deviates substantially from the straight lines joining individual pairs of consecutive data points. This insight is the beginning of
During the design of ship, engineers compute a number of points of a level line of its body and then draw curves through them by means of flexible rulers weighted down at the data points. In this way, the gradients of the curves linking adjacent segments at the data points are the same. Theoretically, one proceeds as follows:
Consider a segment <xi,xi+1> in which, by Lagrange, the straight line is given by:
where
Now assume that one is also given y" at these points, i.e., that one must fit a cubic to the data yi,yi+1, y"i, y"i+1, so that
One must then find C and D in terms of x so that Y" will have the values y"j and y"j+1, respectively. It can be assumed that C and D will be functions of A which are linear in x, i.e.,
However, C(A) must also vanish when A = 0, whence d = 0. Differentiation with respect to A yields
whence b = 0, a = 1/6 and
However, for A = 1, one must have C(A) = 0, whence c = -1/6 and
In a similar manner, one finds
One must now change from A, B to x:
whence
and the final interpolation formula is
Check now this formula. Since A = 1, B = 0 at xj, A = 0, B = 1 at xj+1, the formula is correct for the original interpolation function. Next, differentiate Y(x) with respect to x:
as was intended. If one chooses the values of y" at the data points arbitrarily, the first derivative of Y(x) will be discontinuous at these points. There are n+1 data points, on n-1 of which one can impose continuity conditions for :
There are n-1 such conditions, but n+1 values of Y"(x). Thus, one gives two of them, namely y"0 and y"n arbitrary values. If one sets these two values equal to 0, we arrives a natural cubic spline. Alternatively, onecan specify the value of y' there and compute the values of y" from the formula for Y'.
Compute the natural spline for the data
Starting from j = 3, one finds the system of three linear simultaneous equations
Cramer's rule yields
Hence the interpolation function has the expressions:
Complete these calculations and graph the splines, using Pascal or any other language you know.
In experimental work, one often has repeated runs and many data to which one tries to fit a simple curve, a straight line, an exponential, etc. One uses then the least square method.
Its principle is to form a functional involving the data and parameters specifying the approximation and to minimize it with respect to the parameters of the approximation:
Here, the squares of the errors are used so that they cannot cancel each other. The choice of the approximation Y(x) determines the type of least square approach used. In general, one writes
where the functions fk(x) must satisfy specific conditions, socalled conditions of linear independence. Powers of x are such functions, because you cannot express xl in terms of xk, k=0,1,2, , l-1, whence
The conditions of minimizing the function I are then
They yield the system of N+1 equations which is linear in the parameters ak:
or 
.
Given the data
fit a straight line by the least square method!
In this case, n = 3, N = 2:
The system of two linear equations for a0, a1 is
whence a1 = 2, a0 = 2/3 and the least square line through these data is
which yields at the data points the errors y - Y = 1/3, -2/3, 1/3.
.
