NAG FL Interface
e01baf (dim1_spline)

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 <nag.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)

The routine may be called by the names e01baf or nagf_interp_dim1_spline.

3 Description

e01baf determines a cubic spline

s (x)

, defined in the range

x_{1} \leq x \leq x_{m}

, which interpolates (passes exactly through) the set of data points

(x_{i}, y_{i})

, for

i = 1, 2, \dots, m

, where

m \geq 4

and

x_{1} < x_{2} < \dots < x_{m}

. Unlike some other spline interpolation algorithms, derivative end conditions are not imposed. The spline interpolant chosen has

m - 4

interior knots

λ_{5}, λ_{6}, \dots, λ_{m}

, which are set to the values of

x_{3}, x_{4}, \dots, x_{m - 2}

respectively. This spline is represented in its B-spline form (see Cox (1975)):

s (x) = \sum_{i = 1}^{m} c_{i} N_{i} (x),

where

N_{i} (x)

denotes the normalized B-spline of degree

3

, defined upon the knots

λ_{i}, λ_{i + 1}, \dots, λ_{i + 4}

, and

c_{i}

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

x_{1}

and the last four to

x_{m}

The algorithm for determining the coefficients is as described in Cox (1975) except that

Q R

factorization is used instead of

L U

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: $m$ – Integer Input: On entry: $m$ , the number of data points.

Constraint: $m \geq 4$ .
2: $x (m)$ – Real (Kind=nag_wp) array Input: On entry: $x (i)$ must be set to $x_{i}$ , the $i$ th data value of the independent variable $x$ , for $i = 1, 2, \dots, m$ .

Constraint: $x (i) < x (i + 1)$ , for $i = 1, 2, \dots, m - 1$ .
3: $y (m)$ – Real (Kind=nag_wp) array Input: On entry: $y (i)$ must be set to $y_{i}$ , the $i$ th data value of the dependent variable $y$ , for $i = 1, 2, \dots, m$ .
4: $lamda (lck)$ – Real (Kind=nag_wp) array Output: On exit: the value of $λ_{i}$ , the $i$ th knot, for $i = 1, 2, \dots, m + 4$ .
5: $c (lck)$ – Real (Kind=nag_wp) array Output: On exit: the coefficient $c_{i}$ of the B-spline $N_{i} (x)$ , for $i = 1, 2, \dots, m$ . The remaining elements of the array are not used.
6: $lck$ – Integer Input: On entry: the dimension of the arrays lamda and c as declared in the (sub)program from which e01baf is called.

Constraint: $lck \geq m + 4$ .
7: $wrk (lwrk)$ – Real (Kind=nag_wp) array Workspace
8: $lwrk$ – Integer Input: On entry: the dimension of the array wrk as declared in the (sub)program from which e01baf is called.

Constraint: $lwrk \geq 6 \times m + 16$ .
9: $ifail$ – Integer Input/Output: On entry: ifail must be set to $0$ , $- 1$ or $1$ to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of $0$ causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of $- 1$ means that an error message is printed while a value of $1$ means that it is not.

If halting is not appropriate, the value $- 1$ or $1$ is recommended. If message printing is undesirable, then the value $1$ is recommended. Otherwise, the value $0$ is recommended. 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).

6 Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: On entry, $lck = 〈value〉$ and $m = 〈value〉$ .
Constraint: $lck \geq m + 4$ .

On entry, lwrk is too small. $lwrk = 〈value〉$ . Minimum possible dimension: $〈value〉$ .

On entry, $m = 〈value〉$ .
Constraint: $m \geq 4$ .

$ifail = 2$: On entry, $I = 〈value〉$ , $x (I) = 〈value〉$ and $x (I - 1) = 〈value〉$ .
Constraint: $x (I) > x (I - 1)$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface 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} + δ y_{i}

. The ratio of the root-mean-square value of the

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

9 Further Comments

The time taken by e01baf is approximately proportional to

m

All the

x_{i}

are used as knot positions except

x_{2}

and

x_{m - 1}

. This choice of knots (see Cox (1977)) means that

s (x)

is composed of

m - 3

cubic arcs as follows. If

m = 4

, there is just a single arc space spanning the whole interval

x_{1}

x_{4}

. If

m \geq 5

, the first and last arcs span the intervals

x_{1}

x_{3}

and

x_{m - 2}

x_{m}

respectively. Additionally if

m \geq 6

, the

i

th arc, for

i = 2, 3, \dots, m - 4

, spans the interval

x_{i + 1}

x_{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 (x)

The value of

s (x)

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

s (x)

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

s (x)

over the range

x_{1}

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

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