e02jd:: Curve and Surface Fitting (NAG Toolbox)

function e02jd_example


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

npts  = 15;
xdata = [ 0.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(npts + 4);
ixloc = zeros(numel(x), 1, 'int64');
wrk   = zeros(4*npts + 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 random vectors x, y. 
  % These are used by bivariate function of R. Franke to create data set.
  % The remaining input data 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 and y values
  [state, x, ifail] = g05sa(n, state);
  [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.50*exp(-((9*x(:)-7).^2     + (9*y(:)-3).^2)/4) - ...
      0.20*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

  fig1 = figure;
  meshc(yevalm,xevalm,fevalm);
  title({'Bivariate spline fit from scattered data', ...
         'using two-stage approximation'});
  xlabel('x');
  ylabel('y');
  view(22,28);

$GAUSS$	Gaussian $\exp (- ρ^{2})$
$IMQ$	inverse multiquadric $1 / \sqrt{r^{2} + R^{2}}$
$IMQ2$	inverse multiquadric $1 / (r^{2} + R^{2})$
$IMQ3$	inverse multiquadric $1 / {(r^{2} + R^{2})}^{(3 / 2)}$
$IMQ0_5$	inverse multiquadric $1 / {(r^{2} + R^{2})}^{(1 / 4)}$
$WENDLAND31$	H. Wendland's $C^{2}$ function ${\max (0, 1 - ρ)}^{4} (4 ρ + 1)$
$WENDLAND32$	H. Wendland's $C^{4}$ function ${\max (0, 1 - ρ)}^{6} (35 ρ^{2} + 18 ρ + 3)$
$WENDLAND33$	H. Wendland's $C^{6}$ function ${\max (0, 1 - ρ)}^{8} (32 ρ^{3} + 25 ρ^{2} + 8 ρ + 1)$
$BUHMANNC3$	M. Buhmann's $C^{3}$ function $112 / 45 ρ^{(9 / 2)} + 16 / 3 ρ^{(7 / 2)} - 7 ρ^{4} - 14 / 15 ρ^{2} + 1 / 9$ if $ρ \leq 1$ , $0$ otherwise
$MQ$	multiquadric $\sqrt{r^{2} + R^{2}}$
$MQ1_5$	multiquadric ${(r^{2} + R^{2})}^{(1.5 / 2)}$
$POLYHARMONIC1_5$	polyharmonic spline $ρ^{1.5}$
$POLYHARMONIC1_75$	polyharmonic spline $ρ^{1.75}$

$MQ2$	multiquadric $(r^{2} + R^{2}) \log (r^{2} + R^{2})$
$MQ3$	multiquadric ${(r^{2} + R^{2})}^{(3 / 2)}$
$TPS$	thin plate spline $ρ^{2} \log ρ^{2}$
$POLYHARMONIC3$	polyharmonic spline $ρ^{3}$
$TPS4$	thin plate spline $ρ^{4} \log ρ^{2}$
$POLYHARMONIC5$	polyharmonic spline $ρ^{5}$
$TPS6$	thin plate spline $ρ^{6} \log ρ^{2}$
$POLYHARMONIC7$	polyharmonic spline $ρ^{7}$
$POLYHARMONIC9$	polyharmonic spline $ρ^{9}$

NAG Toolbox: nag_fit_2dspline_ts_sctr (e02jd)

▸▿ Contents

Purpose

Syntax

Description

References

Parameters

Compulsory Input Parameters

Optional Input Parameters

Output Parameters

Error Indicators and Warnings

Accuracy

Further Comments

Example

Optional Parameters

Description of the Optional Parameters