Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox: nag_fit_2dspline_ts_evalm (e02jf)

## Purpose

nag_fit_2dspline_ts_evalm (e02jf) calculates a mesh of values of a spline computed by nag_fit_2dspline_ts_sctr (e02jd).

## Syntax

[fevalm, ifail] = e02jf(xevalm, yevalm, coefs, iopts, opts, 'nxeval', nxeval, 'nyeval', nyeval)
[fevalm, ifail] = nag_fit_2dspline_ts_evalm(xevalm, yevalm, coefs, iopts, opts, 'nxeval', nxeval, 'nyeval', nyeval)

## Description

nag_fit_2dspline_ts_evalm (e02jf) calculates values on a rectangular mesh of a bivariate spline computed by nag_fit_2dspline_ts_sctr (e02jd). The points in the mesh are defined by x$x$ coordinates (xi${x}_{\mathit{i}}$), for i = 1,2,,nx$\mathit{i}=1,2,\dots ,{n}_{x}$, and y$y$ coordinates (yj${y}_{\mathit{j}}$), for j = 1,2,,ny$\mathit{j}=1,2,\dots ,{n}_{y}$. This function is derived from the TSFIT package of O. Davydov and F. Zeilfelder.

## References

Davydov O and Zeilfelder F (2004) Scattered data fitting by direct extension of local polynomials to bivariate splines Advances in Comp. Math. 21 223–271
Farin G and Hansford D (2000) The Essentials of CAGD Natic, MA: A K Peters, Ltd.

## Parameters

### Compulsory Input Parameters

1:     xevalm(nxeval) – double array
nxeval, the dimension of the array, must satisfy the constraint nxeval1${\mathbf{nxeval}}\ge 1$.
The (xi)$\left({x}_{i}\right)$ values forming the mesh on which the spline is to be evaluated.
Constraint: for all i$i$, xevalm(i)${\mathbf{xevalm}}\left(i\right)$ must lie inside, or on the boundary of, the spline's bounding box as determined by nag_fit_2dspline_ts_sctr (e02jd).
2:     yevalm(nyeval) – double array
nyeval, the dimension of the array, must satisfy the constraint nyeval1${\mathbf{nyeval}}\ge 1$.
The (yj)$\left({y}_{j}\right)$ values forming the mesh on which the spline is to be evaluated.
Constraint: for all j$j$, yevalm(j)${\mathbf{yevalm}}\left(j\right)$ must lie inside, or on the boundary of, the spline's bounding box as determined by nag_fit_2dspline_ts_sctr (e02jd).
3:     coefs(lcoefs$\mathit{lcoefs}$) – double array
The computed spline coefficients coefs as output from nag_fit_2dspline_ts_sctr (e02jd).
4:     iopts(liopts$\mathit{liopts}$) – int64int32nag_int array
Must be the same array iopts supplied in a previous call to nag_fit_2dspline_ts_sctr (e02jd). The contents of the array must not have been modified either directly or indirectly, by a call to nag_fit_opt_set (e02zk), between calls to nag_fit_2dspline_ts_sctr (e02jd) and nag_fit_2dspline_ts_evalm (e02jf).
5:     opts(lopts$\mathit{lopts}$) – double array
Must be the same array opts supplied in a previous call to nag_fit_2dspline_ts_sctr (e02jd). The contents of the array must not have been modified either directly or indirectly, by a call to nag_fit_opt_set (e02zk), between calls to nag_fit_2dspline_ts_sctr (e02jd) and nag_fit_2dspline_ts_evalm (e02jf).

### Optional Input Parameters

1:     nxeval – int64int32nag_int scalar
Default: The dimension of the array xevalm.
nx${n}_{x}$, the number of values in the x$x$ direction forming the mesh on which the spline is to be evaluated.
Constraint: nxeval1${\mathbf{nxeval}}\ge 1$.
2:     nyeval – int64int32nag_int scalar
Default: The dimension of the array yevalm.
ny${n}_{y}$, the number of values in the y$y$ direction forming the mesh on which the spline is to be evaluated.
Constraint: nyeval1${\mathbf{nyeval}}\ge 1$.

None.

### Output Parameters

1:     fevalm(nxeval,nyeval) – double array
If ${\mathbf{ifail}}={\mathbf{0}}$ on exit fevalm(i,j)${\mathbf{fevalm}}\left(i,j\right)$ contains the computed spline value at (xi,yj)$\left({x}_{i},{y}_{j}\right)$.
2:     ifail – int64int32nag_int scalar
${\mathrm{ifail}}={\mathbf{0}}$ unless the function detects an error (see [Error Indicators and Warnings]).

## Error Indicators and Warnings

Errors or warnings detected by the function:
ifail = 6${\mathbf{ifail}}=6$
Constraint: nxeval1${\mathbf{nxeval}}\ge 1$.
ifail = 7${\mathbf{ifail}}=7$
Constraint: nyeval1${\mathbf{nyeval}}\ge 1$.
ifail = 9${\mathbf{ifail}}=9$
Option arrays are not initialized or are corrupted.
ifail = 10${\mathbf{ifail}}=10$
The fitting routine has not been called, or the array of coefficients has been corrupted.
ifail = 13${\mathbf{ifail}}=13$
Constraint: _xevalm(i)_$_\le {\mathbf{xevalm}}\left(i\right)\le _$ for all i$i$.
ifail = 14${\mathbf{ifail}}=14$
Constraint: _yevalm(j)_$_\le {\mathbf{yevalm}}\left(j\right)\le _$ for all j$j$.
ifail = 999${\mathbf{ifail}}=-999$
Dynamic memory allocation failed.

## Accuracy

nag_fit_2dspline_ts_evalm (e02jf) uses the de Casteljau algorithm and thus is numerically stable. See Farin and Hansford (2000) for details.

A real array of length O(1)$\mathit{O}\left(1\right)$ is dynamically allocated by each invocation of nag_fit_2dspline_ts_evalm (e02jf).

## Example

```function nag_fit_2dspline_ts_evalm_example
xdata = [0; 0.5; 1; 1.5; 2; 2.5; 3; 4; 4.5; 5; 5.5; 6; 7; 7.5; 8];
ydata = [-1.1; -0.372; 0.431; 1.69; 2.11; 3.1; 4.23; 4.35; 4.81; 4.61; 4.79; ...
5.23; 6.35; 7.19; 7.97];
wdata = [1; 1; 1.5; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1];
cstart = 'c';
sfac  = 0.001;
x     = [6.5178; 7.2463; 1.0159; 7.3070; 5.0589; 0.7803; 2.2280; 4.3751; ...
7.6601; 7.7191; 1.2609; 7.7647; 7.6573; 3.8830; 6.4022; 1.1351; ...
3.3741; 7.3259; 6.3377; 7.6759];
nest  = int64(numel(xdata) + 4);
ixloc = zeros(numel(x), 1, 'int64');
wrk   = zeros(4*numel(xdata) + 16*nest + 41, 1);
iwrk1 = zeros(nest, 1, 'int64');
iwrk2 = zeros(3+3*numel(x), 1, 'int64');
lamda = zeros(nest, 1);
xord  = int64(0);
start = int64(0);
deriv = int64(3);

% Generate the data to fit.
[x, y, f, lsminp, lsmaxp, nxcels, nycels] = generate_data();

% Initialize the options arrays and set/get some options.
[iopts, opts] = handle_options();

% Compute the spline coefficients.
[coefs, iopts, opts, ifail] = ...
nag_fit_2dspline_ts_sctr(x, y, f, lsminp, lsmaxp, nxcels, nycels, iopts, opts);

% pmin and pmax form the bounding box of the spline. We must not attempt to
% evaluate the spline outside this box.
pmin = [min(x); min(y)];
pmax = [max(x); max(y)];

% Evaluate the approximation at a vector of values.
evaluate_at_vector(coefs, iopts, opts, pmin, pmax);

% Evaluate the approximation on a mesh.
evaluate_on_mesh(coefs, iopts, opts, pmin, pmax);

function [x, y, f, lsminp, lsmaxp, nxcels, nycels] = generate_data()
% Generates an x and a y vector of n pseudorandom uniformly distributed
% values on (0,1]. These are passed to the bivariate function of R. Franke
% to create the data set to fit.  The remaining input data for
% nag_fit_2dspline_ts_sctr are set to suitable values for this problem,
% as discussed by Davydov and Zeilfelder.

n = int64(100);

% Initialize the generator to a repeatable sequence
[state, ifail] = nag_rand_init_repeat(int64(1), int64(0), int64(32958));

% Generate x values
[state, x, ifail] = nag_rand_dist_uniform01(n, state);

% Generate y values
[state, y, ifail] = nag_rand_dist_uniform01(n, state);

% Ensure that the bounding box stretches all the way to (0,0) and (1,1)
x(1) = 0;
y(1) = 0;
x(n) = 1;
y(n) = 1;

f = 0.75*exp(-((9*x(:)-2).^2+(9*y(:)-2).^2)/4)+0.75*exp(-(9*x(:)+ 1).^2/49 ...
-(9*y(:)+1)/10) + 0.5*exp(-((9*x(:)-7).^2+(9*y(:)- 3).^2)/4) - ...
0.2*exp(-(9*x(:)- 4).^2-(9.*y(:)-7).^2);

% Grid size for the approximation
nxcels = int64(6);
nycels = int64(6);

% Identify the computation.
fprintf('\nComputing the coefficients of a C^1 spline approximation to Franke''s function\n');
fprintf(' Using a %d by %d grid\n', nxcels, nycels);

% Local-approximation control parameters.
lsminp = int64(3);
lsmaxp = int64(100);
function [iopts, opts] = handle_options()
% Initialize the options arrays and demonstrate how to set and get
% optional parameters.
opts  = zeros(100, 1);
iopts = zeros(100, 1, 'int64');

[iopts, opts, ifail] = ...
nag_fit_opt_set('Initialize = nag_fit_2dspline_ts_sctr', iopts, opts);

%  Set some non-default parameters for the local approximation method.
optstr = strcat('Minimum Singular Value LPA = ', num2str(1/32));
[iopts, opts, ifail] = nag_fit_opt_set(optstr, iopts, opts);
[iopts, opts, ifail] = ...
nag_fit_opt_set('Polynomial Starting Degree = 3', iopts, opts);

% Set a non-default parameter for the global approximation method.
[iopts, opts, ifail] = nag_fit_opt_set('Averaged Spline = Yes', iopts, opts);

% As an example of how to get the value of an optional parameter,
% display whether averaging of local approximations is in operation.
[~, ~, cvalue, ~, ifail] = nag_fit_opt_get('Averaged Spline', iopts, opts);
if strcmp(cvalue, 'YES')
fprintf(' Using an averaged local approximation\n');
end
function evaluate_at_vector(coefs, iopts, opts, pmin, pmax)
% Evaluates the approximation at a (in this case trivial) vector of values.

xevalv = [0];
yevalv = [0];

% Force the points to be within the bounding box of the spline'
for i = 1:numel(xevalv)
xevalv(i) = max(xevalv(i),pmin(1));
xevalv(i) = min(xevalv(i),pmax(1));
yevalv(i) = max(yevalv(i),pmin(2));
yevalv(i) = min(yevalv(i),pmax(2));
end

[fevalv, ifail] = nag_fit_2dspline_ts_evalv(xevalv, yevalv, coefs, iopts, opts);

fprintf('\n Values of computed spline at (x_i,y_i):\n\n');
fprintf('         x_i          y_i   f(x_i,y_i)\n');
for i = 1:numel(xevalv)
fprintf('%12.2f %12.2f %12.2f\n', xevalv(i),yevalv(i),fevalv(i));
end
function evaluate_on_mesh(coefs,iopts,opts,pmin,pmax)
% Evaluates the approximation on a mesh of n_x * n_y values.
nxeval = 101;
nyeval = 101;

% Define the mesh by its lower-left and upper-right corners.
ll_corner = [0; 0];
ur_corner = [1; 1];

% Set the mesh spacing and the evaluation points.
% Force the points to be within the bounding box of the spline.
h = [(ur_corner(1)-ll_corner(1))/(nxeval-1); ...
(ur_corner(2)-ll_corner(2))/(nyeval-1)];

xevalm = ll_corner(1) + [0:nxeval-1]*h(1);
yevalm = ll_corner(2) + [0:nyeval-1]*h(2);

% Ensure that the evaluation points are in the bounding box
xevalm = max(pmin(1), xevalm);
xevalm = min(pmax(1), xevalm);
yevalm = max(pmin(2), yevalm);
yevalm = min(pmax(2), yevalm);

% Evaluate
[fevalm, ifail] = nag_fit_2dspline_ts_evalm(xevalm, yevalm, coefs, iopts, opts);

print_mesh = false;

if print_mesh
fprintf('\nValues of computed spline at (x_i,y_j):\n\n');
fprintf('         x_i          y_i   f(x_i,y_i)\n');
for i = 1:nxeval
for j=1:nyeval
fprintf('%12.2f %12.2f %12.2f\n', xevalm(i),yevalm(j),fevalm(i, j));
end
end
else
fprintf('\nOutputting of the function values on the mesh is disabled\n');
end
```
```

Computing the coefficients of a C^1 spline approximation to Franke's function
Using a 6 by 6 grid
Using an averaged local approximation

Values of computed spline at (x_i,y_i):

x_i          y_i   f(x_i,y_i)
0.00         0.00         0.76

Outputting of the function values on the mesh is disabled

```
```function e02jf_example
xdata = [0; 0.5; 1; 1.5; 2; 2.5; 3; 4; 4.5; 5; 5.5; 6; 7; 7.5; 8];
ydata = [-1.1; -0.372; 0.431; 1.69; 2.11; 3.1; 4.23; 4.35; 4.81; 4.61; 4.79; ...
5.23; 6.35; 7.19; 7.97];
wdata = [1; 1; 1.5; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1];
cstart = 'c';
sfac  = 0.001;
x     = [6.5178; 7.2463; 1.0159; 7.3070; 5.0589; 0.7803; 2.2280; 4.3751; ...
7.6601; 7.7191; 1.2609; 7.7647; 7.6573; 3.8830; 6.4022; 1.1351; ...
3.3741; 7.3259; 6.3377; 7.6759];
nest  = int64(numel(xdata) + 4);
ixloc = zeros(numel(x), 1, 'int64');
wrk   = zeros(4*numel(xdata) + 16*nest + 41, 1);
iwrk1 = zeros(nest, 1, 'int64');
iwrk2 = zeros(3+3*numel(x), 1, 'int64');
lamda = zeros(nest, 1);
xord  = int64(0);
start = int64(0);
deriv = int64(3);

% Generate the data to fit.
[x, y, f, lsminp, lsmaxp, nxcels, nycels] = generate_data();

% Initialize the options arrays and set/get some options.
[iopts, opts] = handle_options();

% Compute the spline coefficients.
[coefs, iopts, opts, ifail] = ...
e02jd(x, y, f, lsminp, lsmaxp, nxcels, nycels, iopts, opts);

% pmin and pmax form the bounding box of the spline. We must not attempt to
% evaluate the spline outside this box.
pmin = [min(x); min(y)];
pmax = [max(x); max(y)];

% Evaluate the approximation at a vector of values.
evaluate_at_vector(coefs, iopts, opts, pmin, pmax);

% Evaluate the approximation on a mesh.
evaluate_on_mesh(coefs, iopts, opts, pmin, pmax);

function [x, y, f, lsminp, lsmaxp, nxcels, nycels] = generate_data()
% Generates an x and a y vector of n pseudorandom uniformly distributed
% values on (0,1]. These are passed to the bivariate function of R. Franke
% to create the data set to fit.  The remaining input data for
% e02jd are set to suitable values for this problem,
% as discussed by Davydov and Zeilfelder.

n = int64(100);

% Initialize the generator to a repeatable sequence
[state, ifail] = g05kf(int64(1), int64(0), int64(32958));

% Generate x values
[state, x, ifail] = g05sa(n, state);

% Generate y values
[state, y, ifail] = g05sa(n, state);

% Ensure that the bounding box stretches all the way to (0,0) and (1,1)
x(1) = 0;
y(1) = 0;
x(n) = 1;
y(n) = 1;

f = 0.75*exp(-((9*x(:)-2).^2+(9*y(:)-2).^2)/4)+0.75*exp(-(9*x(:)+ 1).^2/49 ...
-(9*y(:)+1)/10) + 0.5*exp(-((9*x(:)-7).^2+(9*y(:)- 3).^2)/4) - ...
0.2*exp(-(9*x(:)- 4).^2-(9.*y(:)-7).^2);

% Grid size for the approximation
nxcels = int64(6);
nycels = int64(6);

% Identify the computation.
fprintf('\nComputing the coefficients of a C^1 spline approximation to Franke''s function\n');
fprintf(' Using a %d by %d grid\n', nxcels, nycels);

% Local-approximation control parameters.
lsminp = int64(3);
lsmaxp = int64(100);
function [iopts, opts] = handle_options()
% Initialize the options arrays and demonstrate how to set and get
% optional parameters.
opts  = zeros(100, 1);
iopts = zeros(100, 1, 'int64');

[iopts, opts, ifail] = ...
e02zk('Initialize = e02jd', iopts, opts);

%  Set some non-default parameters for the local approximation method.
optstr = strcat('Minimum Singular Value LPA = ', num2str(1/32));
[iopts, opts, ifail] = e02zk(optstr, iopts, opts);
[iopts, opts, ifail] = ...
e02zk('Polynomial Starting Degree = 3', iopts, opts);

% Set a non-default parameter for the global approximation method.
[iopts, opts, ifail] = e02zk('Averaged Spline = Yes', iopts, opts);

% As an example of how to get the value of an optional parameter,
% display whether averaging of local approximations is in operation.
[~, ~, cvalue, ~, ifail] = e02zl('Averaged Spline', iopts, opts);
if strcmp(cvalue, 'YES')
fprintf(' Using an averaged local approximation\n');
end
function evaluate_at_vector(coefs, iopts, opts, pmin, pmax)
% Evaluates the approximation at a (in this case trivial) vector of values.

xevalv = [0];
yevalv = [0];

% Force the points to be within the bounding box of the spline'
for i = 1:numel(xevalv)
xevalv(i) = max(xevalv(i),pmin(1));
xevalv(i) = min(xevalv(i),pmax(1));
yevalv(i) = max(yevalv(i),pmin(2));
yevalv(i) = min(yevalv(i),pmax(2));
end

[fevalv, ifail] = e02je(xevalv, yevalv, coefs, iopts, opts);

fprintf('\n Values of computed spline at (x_i,y_i):\n\n');
fprintf('         x_i          y_i   f(x_i,y_i)\n');
for i = 1:numel(xevalv)
fprintf('%12.2f %12.2f %12.2f\n', xevalv(i),yevalv(i),fevalv(i));
end
function evaluate_on_mesh(coefs,iopts,opts,pmin,pmax)
% Evaluates the approximation on a mesh of n_x * n_y values.
nxeval = 101;
nyeval = 101;

% Define the mesh by its lower-left and upper-right corners.
ll_corner = [0; 0];
ur_corner = [1; 1];

% Set the mesh spacing and the evaluation points.
% Force the points to be within the bounding box of the spline.
h = [(ur_corner(1)-ll_corner(1))/(nxeval-1); ...
(ur_corner(2)-ll_corner(2))/(nyeval-1)];

xevalm = ll_corner(1) + [0:nxeval-1]*h(1);
yevalm = ll_corner(2) + [0:nyeval-1]*h(2);

% Ensure that the evaluation points are in the bounding box
xevalm = max(pmin(1), xevalm);
xevalm = min(pmax(1), xevalm);
yevalm = max(pmin(2), yevalm);
yevalm = min(pmax(2), yevalm);

% Evaluate
[fevalm, ifail] = e02jf(xevalm, yevalm, coefs, iopts, opts);

print_mesh = false;

if print_mesh
fprintf('\nValues of computed spline at (x_i,y_j):\n\n');
fprintf('         x_i          y_i   f(x_i,y_i)\n');
for i = 1:nxeval
for j=1:nyeval
fprintf('%12.2f %12.2f %12.2f\n', xevalm(i),yevalm(j),fevalm(i, j));
end
end
else
fprintf('\nOutputting of the function values on the mesh is disabled\n');
end
```
```

Computing the coefficients of a C^1 spline approximation to Franke's function
Using a 6 by 6 grid
Using an averaged local approximation

Values of computed spline at (x_i,y_i):

x_i          y_i   f(x_i,y_i)
0.00         0.00         0.76

Outputting of the function values on the mesh is disabled

```