E02AGF (PDF version)
E02 Chapter Contents
E02 Chapter Introduction
NAG Library Manual

NAG Library Routine Document

E02AGF

Note:  before using this routine, please read the Users' Note for your implementation to check the interpretation of bold italicised terms and other implementation-dependent details.

+ Contents

    1  Purpose
    7  Accuracy

1  Purpose

E02AGF computes constrained weighted least squares polynomial approximations in Chebyshev series form to an arbitrary set of data points. The values of the approximations and any number of their derivatives can be specified at selected points.

2  Specification

SUBROUTINE E02AGF ( M, KPLUS1, LDA, XMIN, XMAX, X, Y, W, MF, XF, YF, LYF, IP, A, S, NP1, WRK, LWRK, IWRK, LIWRK, IFAIL)
INTEGER  M, KPLUS1, LDA, MF, LYF, IP(MF), NP1, LWRK, IWRK(LIWRK), LIWRK, IFAIL
REAL (KIND=nag_wp)  XMIN, XMAX, X(M), Y(M), W(M), XF(MF), YF(LYF), A(LDA,KPLUS1), S(KPLUS1), WRK(LWRK)

3  Description

E02AGF determines least squares polynomial approximations of degrees up to k to the set of data points xr,yr with weights wr, for r=1,2,,m. The value of k, the maximum degree required, is to be prescribed by you. At each of the values xfr, for r=1,2,,mf, of the independent variable x, the approximations and their derivatives up to order pr are constrained to have one of the values yfs, for s=1,2,,n, specified by you, where n=mf+r=0mfpr.
The approximation of degree i has the property that, subject to the imposed constraints, it minimizes σi, the sum of the squares of the weighted residuals εr, for r=1,2,,m, where
εr=wryr-fixr
and fixr is the value of the polynomial approximation 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
where xmin and xmax, specified by you, are respectively the lower and upper end points of the interval of x over which the polynomials are to be defined.
The polynomial approximation of degree i can be written as
12ai,0+ai,1T1x-++aijTjx-++aiiTix-
where Tjx- is the Chebyshev polynomial of the first kind of degree j with argument x-. For i=n,n+1,,k, the routine produces the values of the coefficients aij, for j=0,1,,i, together with the value of the root mean square residual,
Si = σ i m +n -i -1 ,
where m is the number of data points with nonzero weight.
Values of the approximations may subsequently be computed using E02AEF or E02AKF.
First E02AGF determines a polynomial μx-, of degree n-1, which satisfies the given constraints, and a polynomial νx-, of degree n, which has value (or derivative) zero wherever a constrained value (or derivative) is specified. It then fits yr-μxr, for r=1,2,,m, with polynomials of the required degree in x- each with factor νx-. Finally the coefficients of μx- are added to the coefficients of these fits to give the coefficients of the constrained polynomial approximations to the data points xr,yr, for r=1,2,,m. The method employed is given in Hayes (1970): it is an extension of Forsythe's orthogonal polynomials method (see Forsythe (1957)) as modified by Clenshaw (see Clenshaw (1960)).

4  References

Clenshaw C W (1960) Curve fitting with a digital computer Comput. J. 2 170–173
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
Hayes J G (ed.) (1970) Numerical Approximation to Functions and Data Athlone Press, London

5  Parameters

1:     M – INTEGERInput
On entry: m, the number of data points to be fitted.
Constraint: M1.
2:     KPLUS1 – INTEGERInput
On entry: k+1, where k is the maximum degree required.
Constraint: n+1KPLUS1m+n is the total number of constraints and m is the number of data points with nonzero weights and distinct abscissae which do not coincide with any of the XFr.
3:     LDA – INTEGERInput
On entry: the first dimension of the array A as declared in the (sub)program from which E02AGF is called.
Constraint: LDAKPLUS1.
4:     XMIN – REAL (KIND=nag_wp)Input
5:     XMAX – REAL (KIND=nag_wp)Input
On entry: the lower and upper end points, respectively, of the interval xmin,xmax. Unless there are specific reasons to the contrary, it is recommended that XMIN and XMAX be set respectively to the lowest and highest value among the xr and xfr. This avoids the danger of extrapolation provided there is a constraint point or data point with nonzero weight at each end point.
Constraint: XMAX>XMIN.
6:     X(M) – REAL (KIND=nag_wp) arrayInput
On entry: Xr must contain the value xr of the independent variable at the rth data point, for r=1,2,,m.
Constraint: the Xr must be in nondecreasing order and satisfy XMINXrXMAX.
7:     Y(M) – REAL (KIND=nag_wp) arrayInput
On entry: Yr must contain yr, the value of the dependent variable at the rth data point, for r=1,2,,m.
8:     W(M) – REAL (KIND=nag_wp) arrayInput
On entry: Wr must contain the weight wr to be applied to the data point xr, for r=1,2,,m. For advice on the choice of weights see the E02 Chapter Introduction. Negative weights are treated as positive. A zero weight causes the corresponding data point to be ignored. Zero weight should be given to any data point whose x and y values both coincide with those of a constraint (otherwise the denominators involved in the root mean square residuals Si will be slightly in error).
9:     MF – INTEGERInput
On entry: mf, the number of values of the independent variable at which a constraint is specified.
Constraint: MF1.
10:   XF(MF) – REAL (KIND=nag_wp) arrayInput
On entry: XFr must contain xfr, the value of the independent variable at which a constraint is specified, for r=1,2,,MF.
Constraint: these values need not be ordered but must be distinct and satisfy XMINXFrXMAX.
11:   YF(LYF) – REAL (KIND=nag_wp) arrayInput
On entry: the values which the approximating polynomials and their derivatives are required to take at the points specified in XF. For each value of XFr, YF contains in successive elements the required value of the approximation, its first derivative, second derivative, ,prth derivative, for r=1,2,,mf. Thus the value, yfs, which the kth derivative of each approximation (k=0 referring to the approximation itself) is required to take at the point XFr must be contained in YFs, where
s=r+k+p1+p2++pr-1,
where k=0,1,,pr and r=1,2,,mf. The derivatives are with respect to the independent variable x.
12:   LYF – INTEGERInput
On entry: the dimension of the array YF as declared in the (sub)program from which E02AGF is called.
Constraint: LYFMF+ i=1 MF IPi.
13:   IP(MF) – INTEGER arrayInput
On entry: IPr must contain pr, the order of the highest-order derivative specified at XFr, for r=1,2,,mf. pr=0 implies that the value of the approximation at XFr is specified, but not that of any derivative.
Constraint: IPr0, for r=1,2,,MF.
14:   A(LDA,KPLUS1) – REAL (KIND=nag_wp) arrayOutput
On exit: Ai+1j+1 contains the coefficient aij in the approximating polynomial of degree i, for i=n,,k and j=0,1,,i.
15:   S(KPLUS1) – REAL (KIND=nag_wp) arrayOutput
On exit: Si+1 contains Si, for i=n,,k, the root mean square residual corresponding to the approximating polynomial of degree i. In the case where the number of data points with nonzero weight is equal to k+1-n, Si is indeterminate: the routine sets it to zero. For the interpretation of the values of Si and their use in selecting an appropriate degree, see Section 3.1 in the E02 Chapter Introduction.
16:   NP1 – INTEGEROutput
On exit: n+1, where n is the total number of constraint conditions imposed: n=mf+p1+p2++pmf.
17:   WRK(LWRK) – REAL (KIND=nag_wp) arrayOutput
On exit: contains weighted residuals of the highest degree of fit determined k. The residual at xr is in element 2n+1+3m+k+1+r, for r=1,2,,m. The rest of the array is used as workspace.
18:   LWRK – INTEGERInput
On entry: the dimension of the array WRK as declared in the (sub)program from which E02AGF is called.
Constraint: LWRKmax4×M+3×KPLUS1,8×n+5×ipmax+MF+10+2×n+2, where ipmax=maxIPr, for r=1,2,,mf.
19:   IWRK(LIWRK) – INTEGER arrayWorkspace
20:   LIWRK – INTEGERInput
On entry: the dimension of the array IWRK as declared in the (sub)program from which E02AGF is called.
Constraint: LIWRK2×MF+2.
21:   IFAIL – INTEGERInput/Output
On entry: IFAIL must be set to 0, -1​ or ​1. 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 -1​ or ​1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this parameter, the recommended value is 0. 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,M<1,
orKPLUS1<n+1,
orLDA<KPLUS1,
orMF<1,
orLYF<n,
orLWRK is too small (see Section 5),
orLIWRK<2×MF+2.
(Here n is the total number of constraint conditions.)
IFAIL=2
IPr<0 for some r=1,2,,MF.
IFAIL=3
XMINXMAX, or XFr is not in the interval XMIN to XMAX for some r=1,2,,MF, or the XFr are not distinct.
IFAIL=4
Xr is not in the interval XMIN to XMAX for some r=1,2,,M.
IFAIL=5
Xr<Xr-1 for some r=2,3,,M.
IFAIL=6
KPLUS1>m+n, where m is the number of data points with nonzero weight and distinct abscissae which do not coincide with any XFr. Thus there is no unique solution.
IFAIL=7
The polynomials μx and/or νx cannot be determined. The problem supplied is too ill-conditioned. This may occur when the constraint points are very close together, or large in number, or when an attempt is made to constrain high-order derivatives.

7  Accuracy

No complete error analysis exists for either the interpolating algorithm or the approximating algorithm. However, considerable experience with the approximating algorithm shows that it is generally extremely satisfactory. Also the moderate number of constraints, of low-order, which are typical of data fitting applications, are unlikely to cause difficulty with the interpolating routine.

8  Further Comments

The time taken to form the interpolating polynomial is approximately proportional to n3, and that to form the approximating polynomials is very approximately proportional to mk+1k+1-n.
To carry out a least squares polynomial fit without constraints, use E02ADF. To carry out polynomial interpolation only, use E01AEF.

9  Example

This example reads data in the following order, using the notation of the parameter list above:
The output is:
The program is written in a generalized form which will read any number of datasets.
The dataset supplied specifies 5 data points in the interval 0.0,4.0 with unit weights, to which are to be fitted polynomials, p, of degrees up to 4, subject to the 3 constraints:

9.1  Program Text

Program Text (e02agfe.f90)

9.2  Program Data

Program Data (e02agfe.d)

9.3  Program Results

Program Results (e02agfe.r)

Produced by GNUPLOT 4.4 patchlevel 0 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0.5 1 1.5 2 2.5 3 -0.004 -0.002 0 0.002 0.004 Polynomial Fit P(x) Residual P(xi)-yi x Example Program Constrained Least-squares Polynomial Approximation residual polynomial fit data points

E02AGF (PDF version)
E02 Chapter Contents
E02 Chapter Introduction
NAG Library Manual

© The Numerical Algorithms Group Ltd, Oxford, UK. 2012