hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_fit_2dcheb_lines (e02ca)

Purpose

nag_fit_2dcheb_lines (e02ca) forms an approximation to the weighted, least squares Chebyshev series surface fit to data arbitrarily distributed on lines parallel to one independent coordinate axis.

Syntax

[a, ifail] = e02ca(m, k, l, x, y, f, w, xmin, xmax, nux, nuy, 'n', n, 'inuxp1', inuxp1, 'inuyp1', inuyp1)
[a, ifail] = nag_fit_2dcheb_lines(m, k, l, x, y, f, w, xmin, xmax, nux, nuy, 'n', n, 'inuxp1', inuxp1, 'inuyp1', inuyp1)

Description

nag_fit_2dcheb_lines (e02ca) determines a bivariate polynomial approximation of degree kk in xx and ll in yy to the set of data points (xr,s,ys,fr,s)(xr,s,ys,fr,s), with weights wr,swr,s, for s = 1,2,,ns=1,2,,n and r = 1,2,,msr=1,2,,ms. That is, the data points are on lines y = ysy=ys, but the xx values may be different on each line. The values of kk and ll are prescribed by you (for guidance on their choice, see Section [Further Comments]). The function is based on the method described in Sections 5 and 6 of Clenshaw and Hayes (1965).
The polynomial is represented in double Chebyshev series form with arguments xx- and yy-. The arguments lie in the range 1-1 to + 1+1 and are related to the original variables xx and yy by the transformations
x = (2x(xmax + xmin))/((xmaxxmin))  and  y = (2y(ymax + ymin))/((ymaxymin)).
x-=2x-(xmax+xmin) (xmax-xmin)   and  y-=2y-(ymax+ymin) (ymax-ymin) .
Here ymaxymax and yminymin are set by the function to, respectively, the largest and smallest value of ysys, but xmaxxmax and xminxmin are functions of yy prescribed by you (see Section [Further Comments]). For this function, only their values xmax(s) x max (s)  and xmin(s) x min (s)  at each y = ysy=ys are required. For each s = 1,2,,ns=1,2,,n, xmax(s) x max (s)  must not be less than the largest xr,sxr,s on the line y = ysy=ys, and, similarly, xmin(s) x min (s)  must not be greater than the smallest xr,sxr,s.
The double Chebyshev series can be written as
kl
aijTi(x)Tj(y)
i = 0j = 0
i=0kj=0laijTi(x-)Tj(y-)
where Ti(x)Ti(x-) is the Chebyshev polynomial of the first kind of degree ii with argument xx-, and Tj(y)Tj(y) is similarly defined. However, the standard convention, followed in this function, is that coefficients in the above expression which have either ii or jj zero are written as (1/2)aij12aij, instead of simply aijaij, and the coefficient with both ii and jj equal to zero is written as (1/4)a0,014a0,0. The series with coefficients output by the function should be summed using this convention. nag_fit_2dcheb_eval (e02cb) is available to compute values of the fitted function from these coefficients.
The function first obtains Chebyshev series coefficients cs,ics,i, for i = 0,1,,ki=0,1,,k, of the weighted least squares polynomial curve fit of degree kk in xx- to the data on each line y = ysy=ys, for s = 1,2,,ns=1,2,,n, in turn, using an auxiliary function. The same function is then called k + 1k+1 times to fit cs,ics,i, for s = 1,2,,ns=1,2,,n, by a polynomial of degree ll in yy-, for each i = 0,1,,ki=0,1,,k. The resulting coefficients are the required aijaij.
You can force the fit to contain a given polynomial factor. This allows for the surface fit to be constrained to have specified values and derivatives along the boundaries x = xminx=xmin, x = xmaxx=xmax, y = yminy=ymin and y = ymaxy=ymax or indeed along any lines x = x-= constant or y = y-= constant (see Section 8 of Clenshaw and Hayes (1965)).

References

Clenshaw C W and Hayes J G (1965) Curve and surface fitting J. Inst. Math. Appl. 1 164–183
Hayes J G (ed.) (1970) Numerical Approximation to Functions and Data Athlone Press, London

Parameters

Compulsory Input Parameters

1:     m(n) – int64int32nag_int array
n, the dimension of the array, must satisfy the constraint n > 0n>0.
m(s)ms must be set to msms, the number of data xx values on the line y = ysy=ys, for s = 1,2,,ns=1,2,,n.
Constraint: m(s) > 0ms>0, for s = 1,2,,ns=1,2,,n.
2:     k – int64int32nag_int scalar
kk, the required degree of xx in the fit.
Constraint: for s = 1,2,,ns=1,2,,n, inuxp11k < mdist(s) + inuxp11inuxp1-1k<mdist(s)+inuxp1-1, where mdist(s)mdist(s) is the number of distinct xx values with nonzero weight on the line y = ysy=ys. See Section [Further Comments].
3:     l – int64int32nag_int scalar
ll, the required degree of yy in the fit.
Constraints:
  • l0l0;
  • inuyp11l < n + inuyp11inuyp1-1l<n+inuyp1-1.
4:     x(mtot) – double array
mtot, the dimension of the array, must satisfy the constraint mtot s = 1n m(s) mtot s=1 n ms .
The xx values of the data points. The sequence must be
  • all points on y = y1y=y1, followed by
  • all points on y = y2y=y2, followed by
  • all points on y = yny=yn.
Constraint: for each ysys, the xx values must be in nondecreasing order.
5:     y(n) – double array
n, the dimension of the array, must satisfy the constraint n > 0n>0.
y(s)ys must contain the yy value of line y = ysy=ys, for s = 1,2,,ns=1,2,,n, on which data is given.
Constraint: the ysys values must be in strictly increasing order.
6:     f(mtot) – double array
mtot, the dimension of the array, must satisfy the constraint mtot s = 1n m(s) mtot s=1 n ms .
ff, the data values of the dependent variable in the same sequence as the xx values.
7:     w(mtot) – double array
mtot, the dimension of the array, must satisfy the constraint mtot s = 1n m(s) mtot s=1 n ms .
The weights to be assigned to the data points, in the same sequence as the xx values. These weights should be calculated from estimates of the absolute accuracies of the frfr, expressed as standard deviations, probable errors or some other measure which is of the same dimensions as frfr. Specifically, each wrwr should be inversely proportional to the accuracy estimate of frfr. Often weights all equal to unity will be satisfactory. If a particular weight is zero, the corresponding data point is omitted from the fit.
8:     xmin(n) – double array
n, the dimension of the array, must satisfy the constraint n > 0n>0.
xmin(s)xmins must contain xmin(s)xmin (s) , the lower end of the range of xx on the line y = ysy=ys, for s = 1,2,,ns=1,2,,n. It must not be greater than the lowest data value of xx on the line. Each xmin(s)xmin (s)  is scaled to 1.0-1.0 in the fit. (See also Section [Further Comments].)
9:     xmax(n) – double array
n, the dimension of the array, must satisfy the constraint n > 0n>0.
xmax(s)xmaxs must contain xmax(s) x max (s) , the upper end of the range of xx on the line y = ysy=ys, for s = 1,2,,ns=1,2,,n. It must not be less than the highest data value of xx on the line. Each xmax(s)xmax (s)  is scaled to + 1.0+1.0 in the fit. (See also Section [Further Comments].)
Constraint: xmax(s) > xmin(s)xmaxs>xmins.
10:   nux(inuxp1) – double array
inuxp1, the dimension of the array, must satisfy the constraint 1inuxp1k + 11inuxp1k+1.
nux(i)nuxi must contain the coefficient of the Chebyshev polynomial of degree (i1)(i-1) in xx-, in the Chebyshev series representation of the polynomial factor in xx- which you require the fit to contain, for i = 1,2,,inuxp1i=1,2,,inuxp1. These coefficients are defined according to the standard convention of Section [Description].
Constraint: nux(inuxp1)nuxinuxp1 must be nonzero, unless inuxp1 = 1inuxp1=1, in which case nux is ignored.
11:   nuy(inuyp1) – double array
nuy(i)nuyi must contain the coefficient of the Chebyshev polynomial of degree (i1)(i-1) in yy-, in the Chebyshev series representation of the polynomial factor which you require the fit to contain, for i = 1,2,,inuyp1i=1,2,,inuyp1. These coefficients are defined according to the standard convention of Section [Description].
Constraint: nuy(inuyp1)nuyinuyp1 must be nonzero, unless inuyp1 = 1inuyp1=1, in which case nuy is ignored.

Optional Input Parameters

1:     n – int64int32nag_int scalar
Default: The dimension of the arrays m, y, xmin, xmax. (An error is raised if these dimensions are not equal.)
The number of lines y = y= constant on which data points are given.
Constraint: n > 0n>0.
2:     inuxp1 – int64int32nag_int scalar
Default: The dimension of the array nux.
inux + 1inux+1, where inuxinux is the degree of a polynomial factor in xx- which you require the fit to contain. (See Section [Description], last paragraph.)
If this option is not required, inuxp1 should be set equal to 11.
Constraint: 1inuxp1k + 11inuxp1k+1.
3:     inuyp1 – int64int32nag_int scalar
Default: The dimension of the array nuy.
inuy + 1inuy+1, where inuyinuy is the degree of a polynomial factor in yy- which you require the fit to contain. (See Section [Description], last paragraph.) If this option is not required, inuyp1 should be set equal to 11.

Input Parameters Omitted from the MATLAB Interface

mtot na work nwork

Output Parameters

1:     a(na) – double array
na(k + 1) × (l + 1)na(k+1)×(l+1), the total number of coefficients in the fit.
Contains the Chebyshev coefficients of the fit. a(i × (l + 1) + j)ai×(l+1)+j is the coefficient aijaij of Section [Description] defined according to the standard convention. These coefficients are used by nag_fit_2dcheb_eval (e02cb) to calculate values of the fitted function.
2:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,k or l < 0l<0,
orinuxp1 or inuyp1 < 1inuyp1<1,
orinuxp1 > k + 1inuxp1>k+1,
orinuyp1 > l + 1inuyp1>l+1,
orm(i) < kinuxp1 + 2mi<k-inuxp1+2 for some i = 1,2,,ni=1,2,,n,
orn < linuyp1 + 2n<l-inuyp1+2,
orna is too small,
ornwork is too small,
ormtot is too small.
  ifail = 2ifail=2
xmin(i)xmini and xmax(i)xmaxi do not span the data x values on y = y(i)y=yi for some i = 1,2,,ni=1,2,,n, possibly because xmin(i)xmax(i)xminixmaxi.
  ifail = 3ifail=3
The data x values on y = y(i)y=yi are not nondecreasing for some i = 1,2,,ni=1,2,,n, or the y(i)yi themselves are not strictly increasing.
  ifail = 4ifail=4
The number of distinct x values with nonzero weight on y = y(i)y=yi is less than kinuxp1 + 2k-inuxp1+2 for some i = 1,2,,ni=1,2,,n.
  ifail = 5ifail=5
On entry,nux(inuxp1) = 0.0nuxinuxp1=0.0 and inuxp11inuxp11,
ornuy(inuyp1) = 0.0nuyinuyp1=0.0 and inuyp11inuyp11.

Accuracy

No error analysis for this method has been published. Practical experience with the method, however, is generally extremely satisfactory.

Further Comments

The time taken is approximately proportional to k × (k × mtot + n × l2)k×(k×mtot+n×l2).
The reason for allowing xmaxxmax and xminxmin (which are used to normalize the range of xx) to vary with yy is that unsatisfactory fits can result if the highest (or lowest) data values of the normalized xx on each line y = ysy=ys are not approximately the same. (For an explanation of this phenomenon, see page 176 of Clenshaw and Hayes (1965).) Commonly in practice, the lowest (for example) data values x1,sx1,s, while not being approximately constant, do lie close to some smooth curve in the (x,y)(x,y) plane. Using values from this curve as the values of xminxmin, different in general on each line, causes the lowest transformed data values x1,sx-1,s to be approximately constant. Sometimes, appropriate curves for xmaxxmax and xminxmin will be clear from the context of the problem (they need not be polynomials). If this is not the case, suitable curves can often be obtained by fitting to the lowest data values x1,sx1,s and to the corresponding highest data values of xx, low degree polynomials in yy, using function nag_fit_1dcheb_arb (e02ad), and then shifting the two curves outwards by a small amount so that they just contain all the data between them. The complete curves are not in fact supplied to the present function, only their values at each ysys; and the values simply need to lie on smooth curves. More values on the complete curves will be required subsequently, when computing values of the fitted surface at arbitrary yy values.
Naturally, a satisfactory approximation to the surface underlying the data cannot be expected if the character of the surface is not adequately represented by the data. Also, as always with polynomials, the approximating function may exhibit unwanted oscillations (particularly near the ends of the ranges) if the degrees kk and ll are taken greater than certain values, generally unknown but depending on the total number of coefficients (k + 1) × (l + 1)(k+1)×(l+1) should be significantly smaller than, say not more than half, the total number of data points. Similarly, k + 1k+1 should be significantly smaller than most (preferably all) the msms, and l + 1l+1 significantly smaller than nn. Closer spacing of the data near the ends of the xx and yy ranges is an advantage. In particular, if ys = cos(π(s1) / (n1)) y-s = - cos( π(s-1)/ (n-1) ) , for s = 1,2,,ns=1,2,,n and xr,s = cos(π(r1) / (m1)) x-r,s = - cos( π(r-1)/ (m-1) ) , for r = 1,2,,mr=1,2,,m, (thus ms = mms=m for all ss), then the values k = m1k=m-1 and l = n1l=n-1 (so that the polynomial passes exactly through all the data points) should not give unwanted oscillations. Other datasets should be similarly satisfactory if they are everywhere at least as closely spaced as the above cosine values with mm replaced by k + 1k+1 and nn by l + 1l+1 (more precisely, if for every ss the largest interval between consecutive values of arccosxr,sarccosx-r,s, for r = 1,2,,mr=1,2,,m, is not greater than π / kπ/k, and similarly for the ysy-s). The polynomial obtained should always be examined graphically before acceptance. Note that, for this purpose it is not sufficient to plot the polynomial only at the data values of xx and yy: intermediate values should also be plotted, preferably via a graphics facility.
Provided the data are adequate, and the surface underlying the data is of a form that can be represented by a polynomial of the chosen degrees, the function should produce a good approximation to this surface. It is not, however, the true least squares surface fit nor even a polynomial in xx and yy, the original variables (see Section 6 of Clenshaw and Hayes (1965), ), except in certain special cases. The most important of these is where the data values of xx are the same on each line y = ysy=ys, (i.e., the data points lie on a rectangular mesh in the (x,y)(x,y) plane), the weights of the data points are all equal, and xmaxxmax and xminxmin are both constants (in this case they should be set to the largest and smallest data values of xx, respectively).
If the dataset is such that it can be satisfactorily approximated by a polynomial of degrees kk and ll, say, then if higher values are used for kk and ll in the function, all the coefficients aijaij for i > ki>k or j > lj>l will take apparently random values within a range bounded by the size of the data errors, or rather less. (This behaviour of the Chebyshev coefficients, most readily observed if they are set out in a rectangular array, closely parallels that in curve-fitting, examples of which are given in Section 8 of Hayes (1970).) In practice, therefore, to establish suitable values of kk and ll, you should first be seeking (within the limitations discussed above) values for kk and ll which are large enough to exhibit the behaviour described. Values for kk and ll should then be chosen as the smallest which do not exclude any coefficients significantly larger than the random ones. A polynomial of degrees kk and ll should then be fitted to the data.
If the option to force the fit to contain a given polynomial factor in xx is used and if zeros of the chosen factor coincide with data xx values on any line, then the effective number of data points on that line is reduced by the number of such coincidences. A similar consideration applies when forcing the yy-direction. No account is taken of this by the function when testing that the degrees kk and ll have not been chosen too large.

Example

function nag_fit_2dcheb_lines_example
m = [int64(8);7;7;6];
k = int64(3);
l = int64(2);
x = [0.1;
     1;
     1.6;
     2.1;
     3.3;
     3.9;
     4.2;
     4.9;
     0.1;
     1.1;
     1.9;
     2.7;
     3.2;
     4.1;
     4.5;
     0.5;
     1.1;
     1.3;
     2.2;
     2.9;
     3.5;
     3.9;
     1.7;
     2;
     2.4;
     2.7;
     3.1;
     3.5 ];
y = [0;
     1;
     2;
     4];
f = [1.01005;
     1.10517;
     1.17351;
     1.23368;
     1.39097;
     1.47698;
     1.52196;
     1.63232;
     2.0201;
     2.23256;
     2.4185;
     2.61993;
     2.75426;
     3.01364;
     3.13662;
     3.15381;
     3.34883;
     3.41649;
     3.73823;
     4.00928;
     4.2572;
     4.43094;
     5.92652;
     6.10701;
     6.35625;
     6.54982;
     6.81713;
     7.09534 ];
w =  ones(28, 1);
xmin = [0;
     0.1;
     0.4;
     1.6];
xmax = [5;
     4.5;
     4;
     3.5];
nux = [0];
nuy = [0];
[a, ifail] = nag_fit_2dcheb_lines(m, k, l, x, y, f, w, xmin, xmax, nux, nuy)
 

a =

   15.3482
    5.1507
    0.1014
    1.1472
    0.1442
   -0.1046
    0.0490
   -0.0031
   -0.0070
    0.0015
   -0.0003
   -0.0002


ifail =

                    0


function e02ca_example
m = [int64(8);7;7;6];
k = int64(3);
l = int64(2);
x = [0.1;
     1;
     1.6;
     2.1;
     3.3;
     3.9;
     4.2;
     4.9;
     0.1;
     1.1;
     1.9;
     2.7;
     3.2;
     4.1;
     4.5;
     0.5;
     1.1;
     1.3;
     2.2;
     2.9;
     3.5;
     3.9;
     1.7;
     2;
     2.4;
     2.7;
     3.1;
     3.5 ];
y = [0;
     1;
     2;
     4];
f = [1.01005;
     1.10517;
     1.17351;
     1.23368;
     1.39097;
     1.47698;
     1.52196;
     1.63232;
     2.0201;
     2.23256;
     2.4185;
     2.61993;
     2.75426;
     3.01364;
     3.13662;
     3.15381;
     3.34883;
     3.41649;
     3.73823;
     4.00928;
     4.2572;
     4.43094;
     5.92652;
     6.10701;
     6.35625;
     6.54982;
     6.81713;
     7.09534 ];
w =  ones(28, 1);
xmin = [0;
     0.1;
     0.4;
     1.6];
xmax = [5;
     4.5;
     4;
     3.5];
nux = [0];
nuy = [0];
[a, ifail] = e02ca(m, k, l, x, y, f, w, xmin, xmax, nux, nuy)
 

a =

   15.3482
    5.1507
    0.1014
    1.1472
    0.1442
   -0.1046
    0.0490
   -0.0031
   -0.0070
    0.0015
   -0.0003
   -0.0002


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013