E02CAF determines a bivariate polynomial approximation of degree in and in to the set of data points , with weights , for and . That is, the data points are on lines , but the values may be different on each line. The values of and are prescribed by you (for guidance on their choice, see Section 8). The subroutine 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 and . The arguments lie in the range to and are related to the original variables and by the transformations
Here and are set by the subroutine to, respectively, the largest and smallest value of , but and are functions of prescribed by you (see Section 8). For this subroutine, only their values and at each are required. For each , must not be less than the largest on the line , and, similarly, must not be greater than the smallest .
The double Chebyshev series can be written as
where is the Chebyshev polynomial of the first kind of degree with argument , and is similarly defined. However, the standard convention, followed in this subroutine, is that coefficients in the above expression which have either or zero are written as , instead of simply , and the coefficient with both and equal to zero is written as . The series with coefficients output by the subroutine should be summed using this convention. E02CBF is available to compute values of the fitted function from these coefficients.
The subroutine first obtains Chebyshev series coefficients , for , of the weighted least squares polynomial curve fit of degree in to the data on each line , for , in turn, using an auxiliary subroutine. The same subroutine is then called times to fit , for , by a polynomial of degree in , for each . The resulting coefficients are the required .
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 , , and or indeed along any lines constant or constant (see Section 8 of Clenshaw and Hayes (1965)).
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
On entry: the weights to be assigned to the data points, in the same sequence as the values. These weights should be calculated from estimates of the absolute accuracies of the , expressed as standard deviations, probable errors or some other measure which is of the same dimensions as . Specifically, each should be inversely proportional to the accuracy estimate of . Often weights all equal to unity will be satisfactory. If a particular weight is zero, the corresponding data point is omitted from the fit.
9: MTOT – INTEGERInput
On entry: the dimension of the arrays X, F and W as declared in the (sub)program from which E02CAF is called.
On exit: contains the Chebyshev coefficients of the fit. is the coefficient of Section 3 defined according to the standard convention. These coefficients are used by E02CBF to calculate values of the fitted function.
11: NA – INTEGERInput
On entry: the dimension of the array A as declared in the (sub)program from which E02CAF is called.
, the total number of coefficients in the fit.
On entry: must contain the coefficient of the Chebyshev polynomial of degree in , in the Chebyshev series representation of the polynomial factor in which you require the fit to contain, for . These coefficients are defined according to the standard convention of Section 3.
must be nonzero, unless , in which case NUX is ignored.
15: INUXP1 – INTEGERInput
On entry: , where is the degree of a polynomial factor in which you require the fit to contain. (See Section 3, last paragraph.)
If this option is not required, INUXP1 should be set equal to .
On entry: must contain the coefficient of the Chebyshev polynomial of degree in , in the Chebyshev series representation of the polynomial factor which you require the fit to contain, for . These coefficients are defined according to the standard convention of Section 3.
must be nonzero, unless , in which case NUY is ignored.
17: INUYP1 – INTEGERInput
On entry: , where is the degree of a polynomial factor in which you require the fit to contain. (See Section 3, last paragraph.) If this option is not required, INUYP1 should be set equal to .
18: WORK(NWORK) – REAL (KIND=nag_wp) arrayWorkspace
19: NWORK – INTEGERInput
On entry: the dimension of the array WORK as declared in the (sub)program from which E02CAF is called.
20: IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to , . If you are unfamiliar with this parameter you should refer to Section 3.3 in the Essential Introduction for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value is recommended. If the output of error messages is undesirable, then the value is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is . When the value is used it is essential to test the value of IFAIL on exit.
On exit: unless the routine detects an error or a warning has been flagged (see Section 6).
6 Error Indicators and Warnings
If on entry or , explanatory error messages are output on the current error message unit (as defined by X04AAF).
and do not span the data X values on for some , possibly because .
The data X values on are not nondecreasing for some , or the themselves are not strictly increasing.
The number of distinct X values with nonzero weight on is less than for some .
No error analysis for this method has been published. Practical experience with the method, however, is generally extremely satisfactory.
8 Further Comments
The time taken is approximately proportional to .
The reason for allowing and (which are used to normalize the range of ) to vary with is that unsatisfactory fits can result if the highest (or lowest) data values of the normalized on each line 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 , while not being approximately constant, do lie close to some smooth curve in the plane. Using values from this curve as the values of , different in general on each line, causes the lowest transformed data values to be approximately constant. Sometimes, appropriate curves for and 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 and to the corresponding highest data values of , low degree polynomials in , using routine E02ADF, 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 subroutine, only their values at each ; 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 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 and are taken greater than certain values, generally unknown but depending on the total number of coefficients should be significantly smaller than, say not more than half, the total number of data points. Similarly, should be significantly smaller than most (preferably all) the , and significantly smaller than . Closer spacing of the data near the ends of the and ranges is an advantage. In particular, if , for and , for , (thus for all ), then the values and (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 replaced by and by (more precisely, if for every the largest interval between consecutive values of , for , is not greater than , and similarly for the ). 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 and : 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 subroutine should produce a good approximation to this surface. It is not, however, the true least squares surface fit nor even a polynomial in and , 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 are the same on each line , (i.e., the data points lie on a rectangular mesh in the plane), the weights of the data points are all equal, and and are both constants (in this case they should be set to the largest and smallest data values of , respectively).
If the dataset is such that it can be satisfactorily approximated by a polynomial of degrees and , say, then if higher values are used for and in the subroutine, all the coefficients for or 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 and , you should first be seeking (within the limitations discussed above) values for and which are large enough to exhibit the behaviour described. Values for and should then be chosen as the smallest which do not exclude any coefficients significantly larger than the random ones. A polynomial of degrees and should then be fitted to the data.
If the option to force the fit to contain a given polynomial factor in is used and if zeros of the chosen factor coincide with data 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 -direction. No account is taken of this by the subroutine when testing that the degrees and have not been chosen too large.
This example reads data in the following order, using the notation of the parameter list for E02CAF above:
The data points are fitted using E02CAF, and then the fitting polynomial is evaluated at the data points using E02CBF.