e01da:: Interpolation (NAG Toolbox)

nag_interp_2d_spline_grid (e01da) determines a bicubic spline interpolant to the set of data points

(x_{q}, y_{r}, f_{q, r})

, for

q = 1, 2, \dots, m_{x}

and

r = 1, 2, \dots, m_{y}

. The spline is given in the B-spline representation

s (x, y) = \sum_{i = 1}^{m_{x}} \sum_{j = 1}^{m_{y}} c_{i j} M_{i} (x) N_{j} (y),

such that

s (x_{q}, y_{r}) = f_{q, r},

where

M_{i} (x)

and

N_{j} (y)

denote normalized cubic B-splines, the former defined on the knots

λ_{i}

λ_{i + 4}

and the latter on the knots

μ_{j}

μ_{j + 4}

, and the

c_{i j}

are the spline coefficients. These knots, as well as the coefficients, are determined by the function, which is derived from the function B2IRE in Anthony et al. (1982). The method used is described in Outline of Method Used.

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

The main sources of rounding errors are in steps

2

3

6

and

7

of the algorithm described in Outline of Method Used. It can be shown (see Cox (1975)) that the matrix

A_{x}

formed in step

2

has elements differing relatively from their true values by at most a small multiple of

3 ε

, where

ε

is the machine precision.

A_{x}

is ‘totally positive’, and a linear system with such a coefficient matrix can be solved quite safely by elimination without pivoting. Similar comments apply to steps

6

and

7

. Thus the complete process is numerically stable.

Further Comments

Timing

Outline of Method Used

The process of computing the spline consists of the following steps:

choice of the interior

x

-knots

λ_{5}

λ_{6}, \dots, λ_{m_{x}}

λ_{i} = x_{i - 2}

, for

i = 5, 6, \dots, m_{x}

formation of the system

A_{x} E = F,

where

A_{x}

is a band matrix of order

m_{x}

and bandwidth

4

, containing in its

q

th row the values at

x_{q}

of the B-splines in

x

f

is the

m_{x}

m_{y}

rectangular matrix of values

f_{q, r}

, and

E

denotes an

m_{x}

m_{y}

rectangular matrix of intermediate coefficients,

use of Gaussian elimination to reduce this system to band triangular form,

solution of this triangular system for

E

choice of the interior

y

knots

μ_{5}

μ_{6}, \dots, μ_{m_{y}}

μ_{i} = y_{i - 2}

, for

i = 5, 6, \dots, m_{y}

formation of the system

A_{y} C^{T} = E^{T},

where

A_{y}

is the counterpart of

A_{x}

for the

y

variable, and

C

denotes the

m_{x}

m_{y}

rectangular matrix of values of

c_{i j}

use of Gaussian elimination to reduce this system to band triangular form,

solution of this triangular system for

C^{T}

and hence

C

Evaluation of Computed Spline

The values of the computed spline at the points

(x_{k}, y_{k})

, for

k = 1, 2, \dots, m

, may be obtained in the double array ff (see nag_fit_2dspline_evalv (e02de)), of length at least

m

, by the following call:

[ff, ifail] = e02de(x, y, lamda, mu, c);

where

M = m

and the coordinates

x_{k}

y_{k}

are stored in

X (k)

Y (k)

. LAMDA, MU and C have the same values as lamda, mu and c output from nag_interp_2d_spline_grid (e01da). (See nag_fit_2dspline_evalv (e02de).)

To evaluate the computed spline on an

m_{x}

m_{y}

rectangular grid of points in the

x

y

plane, which is defined by the

x

coordinates stored in

X (j)

, for

j = 1, 2, \dots, m_{x}

, and the

y

coordinates stored in

Y (k)

, for

k = 1, 2, \dots, m_{y}

, returning the results in the double array ff (see nag_fit_2dspline_evalm (e02df)) which is of length at least

mx \times my

, the following call may be used:

[fg, ifail] = e02df(x, y, lamda, mu, c);

where

MX = m_{x}

MY = m_{y}

. LAMDA, MU and C have the same values as lamda, mu and c output from nag_interp_2d_spline_grid (e01da).

Example

function e01da_example


fprintf('e01da example results\n\n');

x = [1.0   1.10    1.30    1.50    1.60    1.80    2.0];
f = [1.0   1.21    1.69    2.25    2.56    3.24    4.0;   
     1.1   1.31    1.79    2.35    2.66    3.34    4.1; 
     1.4   1.61    2.09    2.65    2.96    3.64    4.4; 
     1.7   1.91    2.39    2.95    3.26    3.94    4.7; 
     1.9   2.11    2.59    3.15    3.46    4.14    4.9;
     2.0   2.21    2.69    3.25    3.56    4.24    5.0];
y = [0.0;
     0.1;
     0.4;
     0.7;
     0.9;
     1.0];

[px, py, lamda, mu, c, ifail] = e01da( ...
                                       x, y, f);

% Display the knot sets, LAMDA and MU.

 fprintf('\n               i   knot lamda(i)      j     knot mu(j)\n');
j = 4:min(px,py)-3;
fprintf('%16d%12.4f%11d%12.4f\n',[j' lamda(j) j' mu(j)]');
if (px>py);
  j = py-2:px-3;
  fprintf('%16d%12.4f\n',[j' lamda(j)]');
elseif (px<py);
  j = px-2:py-3;
  fprintf('%16d%12.4f\n',[j' mu(j)]')
end

% Display the spline coefficients.
c = reshape(c, size(f')); 
fprintf('\n');
disp('The B-Spline coefficients:');
disp(c');

% Evaluate spline on regular 6-by-6 mesh
dx = (x(end)-x(1))/5;
dy = (y(end)-y(1))/5;
tx = [x(1):dx:x(end)];
ty = [y(1):dy:y(end)]';

[ff, ifail] = e02df( ...
                     tx, ty, lamda, mu, c);

% Display evaluations as matrix
ff = reshape(ff,[6,6]);

matrix = 'General';
diag = 'Non-unit';
format = 'F8.3';
title  = 'Spline evaluated on a regular mesh (x across, y down):';
chlab  = 'Character';
rlabs  = cellstr(num2str(tx'));
clabs  = cellstr(num2str(ty));
ncols  = int64(80);
indent = int64(0);
[ifail] =  x04cb( ...
                  matrix, diag, ff, format, title, chlab, ...
                  rlabs, chlab, clabs, ncols, indent);

e01da example results


               i   knot lamda(i)      j     knot mu(j)
               4      1.0000          4      0.0000
               5      1.0000          5      0.0000
               6      2.0000          6      1.0000
               7      2.0000          7      1.0000
               8      2.0000

The B-Spline coefficients:
    1.0000    1.1333    1.3667    1.7000    1.9000    2.0000    1.2000
    1.3333    1.5667    1.9000    2.1000    2.2000    1.5833    1.7167
    1.9500    2.2833    2.4833    2.5833    2.1433    2.2767    2.5100
    2.8433    3.0433    3.1433    2.8667    3.0000    3.2333    3.5667
    3.7667    3.8667    3.4667    3.6000    3.8333    4.1667    4.3667
    4.4667    4.0000    4.1333    4.3667    4.7000    4.9000    5.0000

 Spline evaluated on a regular mesh (x across, y down):
            0     0.2     0.4     0.6     0.8       1
   1    1.000   1.440   1.960   2.560   3.240   4.000
 1.2    1.200   1.640   2.160   2.760   3.440   4.200
 1.4    1.400   1.840   2.360   2.960   3.640   4.400
 1.6    1.600   2.040   2.560   3.160   3.840   4.600
 1.8    1.800   2.240   2.760   3.360   4.040   4.800
   2    2.000   2.440   2.960   3.560   4.240   5.000

On entry,	$mx < 4$ ,
or	$my < 4$ .

NAG Toolbox: nag_interp_2d_spline_grid (e01da)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Timing

Outline of Method Used

Evaluation of Computed Spline

Example