NAG Library Routine Document

e01baf (dim1_spline)


    1  Purpose
    7  Accuracy


e01baf determines a cubic spline interpolant to a given set of data.


Fortran Interface
Subroutine e01baf ( m, x, y, lamda, c, lck, wrk, lwrk, ifail)
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)


e01baf determines a cubic spline sx, defined in the range x1xxm, which interpolates (passes exactly through) the set of data points xi,yi, for i=1,2,,m, where m4 and x1<x2<<xm. Unlike some other spline interpolation algorithms, derivative end conditions are not imposed. The spline interpolant chosen has m-4 interior knots λ5,λ6,,λm, which are set to the values of x3,x4,,xm-2 respectively. This spline is represented in its B-spline form (see Cox (1975)):
where Nix denotes the normalized B-spline of degree 3, defined upon the knots λi,λi+1,,λi+4, and ci denotes its coefficient, whose value is to be determined by the routine.
The use of B-splines requires eight additional knots λ1, λ2, λ3, λ4, λm+1, λm+2, λm+3 and λm+4 to be specified; e01baf sets the first four of these to x1 and the last four to xm.
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.


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


1:     m – IntegerInput
On entry: m, the number of data points.
Constraint: m4.
2:     xm – Real (Kind=nag_wp) arrayInput
On entry: xi must be set to xi, the ith data value of the independent variable x, for i=1,2,,m.
Constraint: xi<xi+1, for i=1,2,,m-1.
3:     ym – Real (Kind=nag_wp) arrayInput
On entry: yi must be set to yi, the ith data value of the dependent variable y, for i=1,2,,m.
4:     lamdalck – Real (Kind=nag_wp) arrayOutput
On exit: the value of λi, the ith knot, for i=1,2,,m+4.
5:     clck – Real (Kind=nag_wp) arrayOutput
On exit: the coefficient ci of the B-spline Nix, for i=1,2,,m. The remaining elements of the array are not used.
6:     lck – IntegerInput
On entry: the dimension of the arrays lamda and c as declared in the (sub)program from which e01baf is called.
Constraint: lckm+4.
7:     wrklwrk – Real (Kind=nag_wp) arrayWorkspace
8:     lwrk – IntegerInput
On entry: the dimension of the array wrk as declared in the (sub)program from which e01baf is called.
Constraint: lwrk6×m+16.
9:     ifail – IntegerInput/Output
On entry: ifail must be set to 0, -1​ 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​ 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 -1​ or ​1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

Error Indicators and Warnings

If on entry 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:
On entry,m<4,
The x-values fail to satisfy the condition
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.
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.
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.


The rounding errors incurred are such that the computed spline is an exact interpolant for a slightly perturbed set of ordinates yi+δyi. The ratio of the root-mean-square value of the δyi to that of the yi is no greater than a small multiple of the relative machine precision.

Parallelism and Performance

e01baf is not threaded in any implementation.

Further Comments

The time taken by e01baf is approximately proportional to m.
All the xi are used as knot positions except x2 and xm-1. This choice of knots (see Cox (1977)) means that sx is composed of m-3 cubic arcs as follows. If m=4, there is just a single arc space spanning the whole interval x1 to x4. If m5, the first and last arcs span the intervals x1 to x3 and xm-2 to xm respectively. Additionally if m6, the ith arc, for i=2,3,,m-4, spans the interval xi+1 to xi+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 sx.
The value of sx at x=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 sx and its first three derivatives at x=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 right-hand value of the third derivative is required (see e02bcf for details).
The value of the integral of sx over the range x1 to xm can be provided in the real variable dint by
Call e02bdf (m+4, lamda, c, dint, ifail)
(see e02bdf).


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 ex are printed out.

Program Text

Program Text (e01bafe.f90)

Program Data


Program Results

Program Results (e01bafe.r)

© The Numerical Algorithms Group Ltd, Oxford, UK. 2017