NAG Library Routine Document
e01baf (dim1_spline)
1
Purpose
e01baf determines a cubic spline interpolant to a given set of data.
2
Specification
Fortran Interface
Integer, Intent (In)  ::  m, lck, lwrk  Integer, Intent (Inout)  ::  ifail  Real (Kind=nag_wp), Intent (In)  ::  x(m), y(m)  Real (Kind=nag_wp), Intent (Out)  ::  lamda(lck), c(lck), wrk(lwrk) 

C Header Interface
#include nagmk26.h
void 
e01baf_ (const Integer *m, const double x[], const double y[], double lamda[], double c[], const Integer *lck, double wrk[], const Integer *lwrk, Integer *ifail) 

3
Description
e01baf determines a cubic spline
$s\left(x\right)$, defined in the range
${x}_{1}\le x\le {x}_{m}$, which interpolates (passes exactly through) the set of data points
$\left({x}_{\mathit{i}},{y}_{\mathit{i}}\right)$, for
$\mathit{i}=1,2,\dots ,m$, where
$m\ge 4$ and
${x}_{1}<{x}_{2}<\cdots <{x}_{m}$. Unlike some other spline interpolation algorithms, derivative end conditions are not imposed. The spline interpolant chosen has
$m4$ interior knots
${\lambda}_{5},{\lambda}_{6},\dots ,{\lambda}_{m}$, which are set to the values of
${x}_{3},{x}_{4},\dots ,{x}_{m2}$ respectively. This spline is represented in its Bspline form (see
Cox (1975)):
where
${N}_{i}\left(x\right)$ denotes the normalized Bspline of degree
$3$,
defined upon the knots
${\lambda}_{i},{\lambda}_{i+1},\dots ,{\lambda}_{i+4}$, and
${c}_{i}$ denotes its coefficient, whose value is to be determined by the routine.
The use of Bsplines requires eight additional knots ${\lambda}_{1}$, ${\lambda}_{2}$,
${\lambda}_{3}$, ${\lambda}_{4}$,
${\lambda}_{m+1}$, ${\lambda}_{m+2}$,
${\lambda}_{m+3}$ and ${\lambda}_{m+4}$ to be specified; e01baf sets the first four of these to ${x}_{1}$ and the last four to ${x}_{m}$.
The algorithm for determining the coefficients is as described in
Cox (1975) except that
$QR$ factorization is used instead of
$LU$ decomposition. The implementation of the algorithm involves setting up appropriate information for the related
routine
e02baf followed by a call of that routine. (See
e02baf for further details.)
Values of the spline interpolant, or of its derivatives or definite integral, can subsequently be computed as detailed in
Section 9.
4
References
Cox M G (1975) An algorithm for spline interpolation J. Inst. Math. Appl. 15 95–108
Cox M G (1977) A survey of numerical methods for data and function approximation The State of the Art in Numerical Analysis (ed D A H Jacobs) 627–668 Academic Press
5
Arguments
 1: $\mathbf{m}$ – IntegerInput

On entry: $m$, the number of data points.
Constraint:
${\mathbf{m}}\ge 4$.
 2: $\mathbf{x}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: ${\mathbf{x}}\left(\mathit{i}\right)$ must be set to ${x}_{\mathit{i}}$, the $\mathit{i}$th data value of the independent variable $x$, for $\mathit{i}=1,2,\dots ,m$.
Constraint:
${\mathbf{x}}\left(\mathit{i}\right)<{\mathbf{x}}\left(\mathit{i}+1\right)$, for $\mathit{i}=1,2,\dots ,{\mathbf{m}}1$.
 3: $\mathbf{y}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) arrayInput

On entry: ${\mathbf{y}}\left(\mathit{i}\right)$ must be set to ${y}_{\mathit{i}}$, the $\mathit{i}$th data value of the dependent variable $y$, for $\mathit{i}=1,2,\dots ,m$.
 4: $\mathbf{lamda}\left({\mathbf{lck}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the value of
${\lambda}_{\mathit{i}}$, the $\mathit{i}$th knot, for $\mathit{i}=1,2,\dots ,m+4$.
 5: $\mathbf{c}\left({\mathbf{lck}}\right)$ – Real (Kind=nag_wp) arrayOutput

On exit: the coefficient
${c}_{\mathit{i}}$ of the Bspline ${N}_{\mathit{i}}\left(x\right)$, for $\mathit{i}=1,2,\dots ,m$. The remaining elements of the array are not used.
 6: $\mathbf{lck}$ – IntegerInput

On entry: the dimension of the arrays
lamda and
c as declared in the (sub)program from which
e01baf is called.
Constraint:
${\mathbf{lck}}\ge {\mathbf{m}}+4$.
 7: $\mathbf{wrk}\left({\mathbf{lwrk}}\right)$ – Real (Kind=nag_wp) arrayWorkspace
 8: $\mathbf{lwrk}$ – IntegerInput

On entry: the dimension of the array
wrk as declared in the (sub)program from which
e01baf is called.
Constraint:
${\mathbf{lwrk}}\ge 6\times {\mathbf{m}}+16$.
 9: $\mathbf{ifail}$ – IntegerInput/Output

On entry:
ifail must be set to
$0$,
$1\text{ or}1$. If you are unfamiliar with this argument you should refer to
Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{ or}1$ is recommended. If the output of error messages is undesirable, then the value
$1$ is recommended. Otherwise, if you are not familiar with this argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{ or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{ifail}}=0$ or
$1$, explanatory error messages are output on the current error message unit (as defined by
x04aaf).
Errors or warnings detected by the routine:
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{m}}<4$, 
or  ${\mathbf{lck}}<{\mathbf{m}}+4$, 
or  ${\mathbf{lwrk}}<6\times {\mathbf{m}}+16$. 
 ${\mathbf{ifail}}=2$

The
xvalues fail to satisfy the condition
${\mathbf{x}}\left(1\right)<{\mathbf{x}}\left(2\right)<{\mathbf{x}}\left(3\right)<\cdots <{\mathbf{x}}\left({\mathbf{m}}\right)$.
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
See
Section 3.9 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
See
Section 3.8 in How to Use the NAG Library and its Documentation for further information.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
See
Section 3.7 in How to Use the NAG Library and its Documentation for further information.
7
Accuracy
The rounding errors incurred are such that the computed spline is an exact interpolant for a slightly perturbed set of ordinates ${y}_{i}+\delta {y}_{i}$. The ratio of the rootmeansquare value of the $\delta {y}_{i}$ to that of the ${y}_{i}$ is no greater than a small multiple of the relative machine precision.
8
Parallelism and Performance
e01baf is not threaded in any implementation.
The time taken by e01baf is approximately proportional to $m$.
All the
${x}_{i}$ are used as knot positions except
${x}_{2}$ and
${x}_{m1}$. This choice of knots
(see
Cox (1977)) means that
$s\left(x\right)$ is composed of
$m3$ cubic arcs as follows. If
$m=4$, there is just a single arc space spanning the whole interval
${x}_{1}$ to
${x}_{4}$. If
$m\ge 5$, the first and last arcs span the intervals
${x}_{1}$ to
${x}_{3}$ and
${x}_{m2}$ to
${x}_{m}$ respectively. Additionally if
$m\ge 6$, the
$\mathit{i}$th arc, for
$\mathit{i}=2,3,\dots ,m4$, spans the interval
${x}_{\mathit{i}+1}$ to
${x}_{\mathit{i}+2}$.
After the call
Call e01baf (m, x, y, lamda, c, lck, wrk, lwrk, ifail)
the following operations may be carried out on the interpolant
$s\left(x\right)$.
The value of
$s\left(x\right)$ at
$x={\mathbf{x}}$ can be provided in the real variable
s by the call
Call e02bbf (m+4, lamda, c, x, s, ifail)
(see
e02bbf).
The values of
$s\left(x\right)$ and its first three derivatives at
$x={\mathbf{x}}$ can be provided in the real array
s
of dimension
$4$, by the call
Call e02bcf (m+4, lamda, c, x, left, s, ifail)
(see
e02bcf).
Here
left
must specify whether the left or righthand value of the third derivative is required (see
e02bcf for details).
The value of the integral of
$s\left(x\right)$ over the range
${x}_{1}$ to
${x}_{m}$ can be provided in the
real variable
dint
by
Call e02bdf (m+4, lamda, c, dint, ifail)
(see
e02bdf).
10
Example
This example sets up data from $7$ values of the exponential function in the interval $0$ to $1$. e01baf is then called to compute a spline interpolant to these data.
The spline is evaluated by
e02bbf, at the data points and at points halfway between each adjacent pair of data points, and the spline values and the values of
${e}^{x}$ are printed out.
10.1
Program Text
Program Text (e01bafe.f90)
10.2
Program Data
None.
10.3
Program Results
Program Results (e01bafe.r)