NAG FL Interface
e02adf (dim1_​cheb_​arb)

1 Purpose

e02adf computes weighted least squares polynomial approximations to an arbitrary set of data points.

2 Specification

Fortran Interface
Subroutine e02adf ( m, kplus1, lda, x, y, w, work1, work2, a, s, ifail)
Integer, Intent (In) :: m, kplus1, lda
Integer, Intent (Inout) :: ifail
Real (Kind=nag_wp), Intent (In) :: x(m), y(m), w(m)
Real (Kind=nag_wp), Intent (Inout) :: a(lda,kplus1)
Real (Kind=nag_wp), Intent (Out) :: work1(3*m), work2(2*kplus1), s(kplus1)
C Header Interface
#include <nag.h>
void  e02adf_ (const Integer *m, const Integer *kplus1, const Integer *lda, const double x[], const double y[], const double w[], double work1[], double work2[], double a[], double s[], Integer *ifail)
The routine may be called by the names e02adf or nagf_fit_dim1_cheb_arb.

3 Description

e02adf determines least squares polynomial approximations of degrees 0,1,,k to the set of data points xr,yr with weights wr, for r=1,2,,m.
The approximation of degree i has the property that it minimizes σi the sum of squares of the weighted residuals εr, where
εr=wryr-fr  
and fr is the value of the polynomial of degree i at the rth data point.
Each polynomial is represented in Chebyshev series form with normalized argument x¯. This argument lies in the range -1 to +1 and is related to the original variable x by the linear transformation
x¯= 2x-xmax-xmin xmax-xmin .  
Here xmax and xmin are respectively the largest and smallest values of xr. The polynomial approximation of degree i is represented as
12ai+1,1T0x¯+ai+1,2T1x¯+ai+1,3T2x¯++ai+1,i+1Tix¯,  
where Tjx¯, for j=0,1,,i, are the Chebyshev polynomials of the first kind of degree j with argument x¯.
For i=0,1,,k, the routine produces the values of ai+1,j+1, for j=0,1,,i, together with the value of the root-mean-square residual si=σi/m-i-1. In the case m=i+1 the routine sets the value of si to zero.
The method employed is due to Forsythe (1957) and is based on the generation of a set of polynomials orthogonal with respect to summation over the normalized dataset. The extensions due to Clenshaw (1960) to represent these polynomials as well as the approximating polynomials in their Chebyshev series forms are incorporated. The modifications suggested by Reinsch and Gentleman (see Gentleman (1969)) to the method originally employed by Clenshaw for evaluating the orthogonal polynomials from their Chebyshev series representations are used to give greater numerical stability.
For further details of the algorithm and its use see Cox (1974) and Cox and Hayes (1973).
Subsequent evaluation of the Chebyshev series representations of the polynomial approximations should be carried out using e02aef.

4 References

Clenshaw C W (1960) Curve fitting with a digital computer Comput. J. 2 170–173
Cox M G (1974) A data-fitting package for the non-specialist user Software for Numerical Mathematics (ed D J Evans) Academic Press
Cox M G and Hayes J G (1973) Curve fitting: a guide and suite of algorithms for the non-specialist user NPL Report NAC26 National Physical Laboratory
Forsythe G E (1957) Generation and use of orthogonal polynomials for data fitting with a digital computer J. Soc. Indust. Appl. Math. 5 74–88
Gentleman W M (1969) An error analysis of Goertzel's (Watt's) method for computing Fourier coefficients Comput. J. 12 160–165
Hayes J G (ed.) (1970) Numerical Approximation to Functions and Data Athlone Press, London

5 Arguments

1: m Integer Input
On entry: the number m of data points.
Constraint: mmdist2, where mdist is the number of distinct x values in the data.
2: kplus1 Integer Input
On entry: k+1, where k is the maximum degree required.
Constraint: 0<kplus1mdist, where mdist is the number of distinct x values in the data.
3: lda Integer Input
On entry: the first dimension of the array a as declared in the (sub)program from which e02adf is called.
Constraint: ldakplus1.
4: xm Real (Kind=nag_wp) array Input
On entry: the values xr of the independent variable, for r=1,2,,m.
Constraint: the values must be supplied in nondecreasing order with xm>x1.
5: ym Real (Kind=nag_wp) array Input
On entry: the values yr of the dependent variable, for r=1,2,,m.
6: wm Real (Kind=nag_wp) array Input
On entry: the set of weights, wr, for r=1,2,,m. For advice on the choice of weights, see Section 2.1.2 in the E02 Chapter Introduction.
Constraint: wr>0.0, for r=1,2,,m.
7: work13×m Real (Kind=nag_wp) array Workspace
8: work22×kplus1 Real (Kind=nag_wp) array Workspace
9: aldakplus1 Real (Kind=nag_wp) array Output
On exit: the coefficients of Tjx¯ in the approximating polynomial of degree i. ai+1j+1 contains the coefficient ai+1,j+1, for i=0,1,,k and j=0,1,,i.
10: skplus1 Real (Kind=nag_wp) array Output
On exit: si+1 contains the root-mean-square residual si, for i=0,1,,k, as described in Section 3. For the interpretation of the values of the si and their use in selecting an appropriate degree, see Section 3.1 in the E02 Chapter Introduction.
11: ifail Integer Input/Output
On entry: ifail must be set to 0, -1 or 1 to set behaviour on detection of an error; these values have no effect when no error is detected.
A value of 0 causes the printing of an error message and program execution will be halted; otherwise program execution continues. A value of -1 means that an error message is printed while a value of 1 means that it is not.
If halting is not appropriate, the value -1 or 1 is recommended. If message printing is undesirable, then the value 1 is recommended. Otherwise, the value 0 is recommended. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, i=value and wi=value.
Constraint: wi>0.0.
ifail=2
On entry, i=value, xi=value and xi-1=value.
Constraint: xixi-1.
ifail=3
On entry, all xI have the same value: x1=value.
ifail=4
On entry, kplus1=value.
Constraint: kplus11.
On entry, kplus1=value and m=value.
Constraint: kplus1m.
On entry, kplus1=value and mdist=value.
Constraint: kplus1mdist, where mdist is the number of distinct x values.
ifail=5
On entry, lda=value and kplus1=value.
Constraint: ldakplus1.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

7 Accuracy

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

8 Parallelism and Performance

e02adf is not threaded in any implementation.

9 Further Comments

The time taken is approximately proportional to mk+1k+11.
The approximating polynomials may exhibit undesirable oscillations (particularly near the ends of the range) if the maximum degree k exceeds a critical value which depends on the number of data points m and their relative positions. As a rough guide, for equally-spaced data, this critical value is about 2×m. For further details see page 60 of Hayes (1970).

10 Example

Determine weighted least squares polynomial approximations of degrees 0, 1, 2 and 3 to a set of 11 prescribed data points. For the approximation of degree 3, tabulate the data and the corresponding values of the approximating polynomial, together with the residual errors, and also the values of the approximating polynomial at points half-way between each pair of adjacent data points.
The example program supplied is written in a general form that will enable polynomial approximations of degrees 0,1,,k to be obtained to m data points, with arbitrary positive weights, and the approximation of degree k to be tabulated. e02aef is used to evaluate the approximating polynomial. The program is self-starting in that any number of datasets can be supplied.

10.1 Program Text

Program Text (e02adfe.f90)

10.2 Program Data

Program Data (e02adfe.d)

10.3 Program Results

Program Results (e02adfe.r)
GnuplotProduced by GNUPLOT 4.6 patchlevel 3 0 10 20 30 40 0 2 4 6 8 10 −0.1 −0.05 0 0.05 0.1 y |P(xi) -f(xi)| x Example Program Least-squares Cubic Polynomial Approximation to a set of 11 Data Points P(xi) -f(xi) polynomial fit gnuplot_plot_1 points f(xi) gnuplot_plot_2 gnuplot_plot_3