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
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
|
