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


Cover Contents Index Main page

last modified 21th August 2001
Copyright 2000-2001 Vladimir Zabrodsky
Czech Republic