E211 Cubic Splines and their Integrals

Subroutines RCSPLN and DCSPLN compute a (vector-valued) cubic spline function $F(x)$ which interpolates between a given set of points. Entries RCSPNT and DCSPNT compute the first and second integral over $F(x)$ .

On computers other than CDC or Cray, only the double-precision versions DCSPLN and DCSPNT are available. On CDC and Cray computers, only the single-precision versions RCSPLN and RCSPNT are available.

Given an interval $[a,b]$ , a subdivision of this interval into $n\; \ge 2$

subintervals

$a\; =\; x$₀< x₁< ... < x_n-1< x_n= b,

and n+1 function values $Y$_k= y_k1,...,y_km on the n+1 abscissae (called `knots') $x$<sub>k</sub> ($k=0,1,...,n$ ), RCSPLN and DCSPLN compute a function $F(x)$ of class $C2$ , defined on $[a,b]$ , which assumes the given value $Y$<sub>k</sub> at the knot $x$<sub>k</sub> (i.e. $F(x$<sub>k</sub>) = Y<sub>k</sub> ), and which, when restricted to the ith sub-interval $[x$<sub>i-1</sub>, x<sub>i</sub>] is identical with a set of m polynomials $p$<sub>i1</sub>,...,p<sub>im</sub> , each of degree at most 3. Any function $F(x)$ which satisfies the above two conditions is called a `cubic spline' through the n + 1 points $(x$_k, Y_k) . To define the spline function $F(x)$ uniquely the subroutines impose an additional boundary condition, specified by their MODE parameter:

MODE = 0: $F\text{'}\text{'}(x$₀) = F''(x_n)=0 (the so-called natural spline).

MODE = 1: $F\text{'}\text{'}(x$₀) = F''(x₁) and $F\text{'}\text{'}(x$_n-1) = F''(x_n) .

Structure:

SUBROUTINE subprograms
User Entry Names: RCSPLN, RCSPNT, DCSPLN, DCSPNT
Files referenced: Unit 6

Usage:

For $t=R$ (type REAL), $t=D$ (type DOUBLE PRECISION),

Spline:	CALL tCSPLN(N,X,M,Y,NDIM,MODE,A,B,C,D)
Integrals:	CALL tCSPNT(N,X,M,Y,NDIM,MODE,A,B,C,D)

N

( INTEGER) Number n of subintervals $[x$_i-1, x_i] . Must contain a value of at least 2 on entry. Unchanged on return.

X

(type according to t) One-dimensional array of dimension (0:d) of at least n + 1 elements. On entry, X(k) must contain the abscissa $x$_k , $(k=0,1,...,n)$ . Unchanged on return.

M

( INTEGER) Number m of components of the vector-valued spline function $F(x)$ . Must contain a value of at least 1 on entry. Unchanged on return.

Y

(type according to t) Two-dimensional array of dimension (0:NDIM,d) where d is a value not less than m. On entry, Y(k,j) must contain the value $y$_kj of the jth component of the vector $Y$_k , $(k=0,1,...,n$ ; $j=1,...,m)$ . Unchanged on return.

NDIM

( INTEGER) Upper bound of the first dimension of arrays A, B, C, D and Y. Must contain a value of at least n on entry. Unchanged on return.

MODE

( INTEGER) Type of boundary condition (see description above). Must contain the value 0 or 1 on entry. Unchanged on return.

A,B,C,D

(type according to t) Two-dimensional arrays of dimension (NDIM,d), where $d\ge m$ .
On return from RCSPLN, A(i,j), B(i,j), C(i,j) and D(i,j) will contain the four coefficients $a$_ij, b_ij, c_ij, and $d$_ij of the polynomial

$p$_ij= a_ij+ b_ij(x - x_i-1) + c_ij(x - x_i-1)²+ d_ij(x - x_i-1)³

that determines the jth component $f$_j(x) of the spline in the ith subinterval $[x$_i-1, x_i] , $i=l,...,n$ , $j=1,...,m$ .
On return from RCSPNT,
A(i,j) = $\int$_a^x_if_j(t) dt and B(i,j) = $\int$_a^x_i∫_a^xf_j(t) dt dx ,
with $i=1,...,n;$ $j=1,...,m$ .
Arrays C and D have been used as working space.

Restrictions:

$N\; \ge 2$ , $M\; \ge 1$ , $NDIM\; \ge N$ , $MODE\; =\; 0$ or 1.

Error handling:

Error E211.1: .
Error E211.2: .
Error E211.3: .
Error E211.4: $MODE\; \ne 0$ and $MODE\; \ne 1$ .
A message is written on Unit 6, unless subroutine MTLSET (N002) has been called.
$•$

E222

K.S. Kölbig