Cubic spline interpolation

PROBLEM
We know the value of a function at a set of points A.1<A.2<...<A.N. Estimate the value of a function for arbitrary V in [A.1,A.N].

DEFINITION
A function S is called a cubic spline on [A.1,A.N] if
• S is defined on [A.1,A.N];
• S, first and second derivative of S are all continuous functions on [A.1,A.N];
• There are points (knots of the spline S) such that A.1<A.2<...<A.N and S is a polynomial of degree <=3 on each subinterval [A.J,A.Jp1], where Jp1=J+1.

ALGORITHM
Given arrays A.1,...,A.N and B.1,...,B.N containing a tabulated function with A.1<A.2<...<A.N, and given values B1A1 and B1AN for the first derivative of the interpolating function at points A.1 and A.N. The routine SPLINE computes array B2. that contains the second derivatives of the interpolating function at the tabulated points. The routine SPLINE is called only once. Values of the interpolated function for any value V are obtained by calls to a spline interpolation function SPLINT.

IMPLEMENTATION of SPLINE
Unit: internal subroutine

Global variables: input ascending sequence of values A.1,...,A.N, input array B. and output array B2. that contains the second derivatives of the interpolating function at the tabulated points

Parameters: positive integer N, values B1A1 and B1AN for the first derivative of the interpolating function at points A.1 and A.N

 SPLINE: procedure expose A. B. B2. N = ARG(1); B1A1 = ARG(2); B1AN = ARG(3) if B1A1 > 0.99E+99 then do; B2.1=0; U.1=0; end   else do     B2.1 = -0.5     U.1=(3 / (A.2 - A.1)) *,       ((B.2 - B.1) / (A.2 - A.1) - B1A1)   end do J = 2 to N - 1   Jm1 = J - 1; Jp1 = J + 1   Sig = (A.J - A.Jm1) / (A.Jp1 - A.Jm1)   P = Sig * B2.Jm1 + 2   B2.J = (Sig - 1) / P   U.J = (B.Jp1 - B.J) / (A.Jp1 - A.J) -,     (B.J - B.Jm1) / (A.J - A.Jm1)   U.J = (6 * U.J / (A.Jp1 - A.Jm1) - Sig * U.Jm1) / P end Nm1 = N - 1 if B1AN > 0.99E+30 then do; Qn = 0; Un = 0; end   else do     Qn = 0.5;     Un = (3 / (A.N - A.Nm1)) *,       (B1AN - (B.N - B.Nm1) / (A.N - A.Nm1))   end B2.N = (Un - Qn * U.Nm1) / (Qn * B2.Nm1 + 1) do J = N - 1 to 1 by -1   Jp1 = J + 1   B2.J = B2.J * B2.Jp1 + U.J end return

IMPLEMENTATION of SPLINT
Unit: internal function

Global variables: input arrays: ascending sequence of values A., array B. and array B2. computed by SPLINE

Parameters: positive integer N, arbitrary value V

Returns: cubic-spline interpolated value of B.V

 SPLINT: procedure expose A. B. B2. N = ARG(1); V = ARG(2) L = 1; H = N do while (H - L) > 1   K = (H + L) % 2   if A.K > V then H = K     else L = K end Q = A.H - A.L if Q = 0 then   call ERROR "SPLINT: Input has to be",     "an ascending sequence" X = (A.H - V) / Q; Y = (V - A.L) / Q BV = (X**3 - X) * B2.L + (Y**3 - Y) * B2.H BV = X * B.L + Y * B.H + BV * (Q**2) / 6 return BV   ERROR: say ARG(1); exit

EXAMPLE
The following program displays on the screen

2.52386364
2.71270431

 do J = 1 to 6; A.J = J; end B.1 = 1.1; B.2 = 2.5; B.3 = 2.6 B.4 = 3.0; B.5 = 5.0; B.6 = 4.0 call SPLINE 6, 0, 0 say SPLINT(6, 3.5); say SPLINT(6, 3.8) exit SPLINE: procedure expose A. B. B2. ... SPLINT: procedure expose A. B. B2. ...

CONNECTIONS

CO-AUTHOR
Doug Rickman - the SPLINE, SPLINT functions and example

Literature
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing
- 2nd ed. University Press, Cambridge, 1992