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) 

C++ Header Interface
#include <nag.h> extern "C" {
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,\dots ,k$ to the set of data points $\left({x}_{\mathit{r}},{y}_{\mathit{r}}\right)$ with weights ${w}_{\mathit{r}}$, for $\mathit{r}=1,2,\dots ,m$.
The approximation of degree
$i$ has the property that it minimizes
${\sigma}_{i}$ the sum of squares of the weighted residuals
${\epsilon}_{r}$, where
and
${f}_{r}$ is the value of the polynomial of degree
$i$ at the
$r$th data point.
Each polynomial is represented in Chebyshev series form with normalized argument
$\overline{x}$. This argument lies in the range
$1$ to
$+1$ and is related to the original variable
$x$ by the linear transformation
Here
${x}_{\mathrm{max}}$ and
${x}_{\mathrm{min}}$ are respectively the largest and smallest values of
${x}_{r}$. The polynomial approximation of degree
$i$ is represented as
where
${T}_{\mathit{j}}\left(\overline{x}\right)$, for
$\mathit{j}=0,1,\dots ,i$, are the Chebyshev polynomials of the first kind of degree
$j$ with argument
$\left(\overline{x}\right)$.
For $i=0,1,\dots ,k$, the routine produces the values of ${a}_{i+1,\mathit{j}+1}$, for $\mathit{j}=0,1,\dots ,i$, together with the value of the rootmeansquare residual ${s}_{i}=\sqrt{{\sigma}_{i}/\left(mi1\right)}$. In the case $m=i+1$ the routine sets the value of ${s}_{i}$ 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 datafitting package for the nonspecialist 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 nonspecialist 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:
$\mathbf{m}$ – Integer
Input

On entry: the number $m$ of data points.
Constraint:
${\mathbf{m}}\ge \mathit{mdist}\ge 2$, where $\mathit{mdist}$ is the number of distinct $x$ values in the data.

2:
$\mathbf{kplus1}$ – Integer
Input

On entry: $k+1$, where $k$ is the maximum degree required.
Constraint:
$0<{\mathbf{kplus1}}\le \mathit{mdist}$, where $\mathit{mdist}$ is the number of distinct $x$ values in the data.

3:
$\mathbf{lda}$ – Integer
Input

On entry: the first dimension of the array
a as declared in the (sub)program from which
e02adf is called.
Constraint:
${\mathbf{lda}}\ge {\mathbf{kplus1}}$.

4:
$\mathbf{x}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) array
Input

On entry: the values
${x}_{\mathit{r}}$ of the independent variable, for $\mathit{r}=1,2,\dots ,m$.
Constraint:
the values must be supplied in nondecreasing order with ${\mathbf{x}}\left({\mathbf{m}}\right)>{\mathbf{x}}\left(1\right)$.

5:
$\mathbf{y}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) array
Input

On entry: the values
${y}_{\mathit{r}}$ of the dependent variable, for $\mathit{r}=1,2,\dots ,m$.

6:
$\mathbf{w}\left({\mathbf{m}}\right)$ – Real (Kind=nag_wp) array
Input

On entry: the set of weights,
${w}_{\mathit{r}}$, for
$\mathit{r}=1,2,\dots ,m$. For advice on the choice of weights, see
Section 2.1.2 in the
E02 Chapter Introduction.
Constraint:
${\mathbf{w}}\left(\mathit{r}\right)>0.0$, for $\mathit{r}=1,2,\dots ,m$.

7:
$\mathbf{work1}\left(3\times {\mathbf{m}}\right)$ – Real (Kind=nag_wp) array
Workspace

8:
$\mathbf{work2}\left(2\times {\mathbf{kplus1}}\right)$ – Real (Kind=nag_wp) array
Workspace


9:
$\mathbf{a}\left({\mathbf{lda}},{\mathbf{kplus1}}\right)$ – Real (Kind=nag_wp) array
Output

On exit: the coefficients of
${T}_{\mathit{j}}\left(\overline{x}\right)$ in the approximating polynomial of degree $\mathit{i}$. ${\mathbf{a}}\left(\mathit{i}+1,\mathit{j}+1\right)$ contains the coefficient ${a}_{\mathit{i}+1,\mathit{j}+1}$, for $\mathit{i}=0,1,\dots ,k$ and $\mathit{j}=0,1,\dots ,i$.

10:
$\mathbf{s}\left({\mathbf{kplus1}}\right)$ – Real (Kind=nag_wp) array
Output

On exit:
${\mathbf{s}}\left(\mathit{i}+1\right)$ contains the rootmeansquare residual
${s}_{\mathit{i}}$, for
$\mathit{i}=0,1,\dots ,k$, as described in
Section 3. For the interpretation of the values of the
${s}_{\mathit{i}}$ and their use in selecting an appropriate degree, see
Section 3.1 in the
E02 Chapter Introduction.

11:
$\mathbf{ifail}$ – Integer
Input/Output

On entry:
ifail must be set to
$0$,
$1\text{or}1$. If you are unfamiliar with this argument you should refer to
Section 4 in the Introduction to the NAG Library FL Interface for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value
$1\text{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 argument, the recommended value is
$0$.
When the value $\mathbf{1}\text{or}\mathbf{1}$ is used it is essential to test the value of ifail on exit.
On exit:
${\mathbf{ifail}}={\mathbf{0}}$ unless the routine detects an error or a warning has been flagged (see
Section 6).
6
Error Indicators and Warnings
If on entry
${\mathbf{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:
 ${\mathbf{ifail}}=1$

On entry, $\mathit{i}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{w}}\left(\mathit{i}\right)=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{w}}\left(\mathit{i}\right)>0.0$.
 ${\mathbf{ifail}}=2$

On entry, $\mathit{i}=\u2329\mathit{\text{value}}\u232a$, ${\mathbf{x}}\left(\mathit{i}\right)=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{x}}\left(\mathit{i}1\right)=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{x}}\left(\mathit{i}\right)\ge {\mathbf{x}}\left(\mathit{i}1\right)$.
 ${\mathbf{ifail}}=3$

On entry, all ${\mathbf{x}}\left(\mathit{I}\right)$ have the same value: ${\mathbf{x}}\left(1\right)=\u2329\mathit{\text{value}}\u232a$.
 ${\mathbf{ifail}}=4$

On entry, ${\mathbf{kplus1}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{kplus1}}\ge 1$.
On entry, ${\mathbf{kplus1}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{m}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{kplus1}}\le {\mathbf{m}}$.
On entry,
${\mathbf{kplus1}}=\u2329\mathit{\text{value}}\u232a$ and
$\mathit{mdist}=\u2329\mathit{\text{value}}\u232a$.
Constraint:
${\mathbf{kplus1}}\le \mathit{mdist}$, where
mdist is the number of distinct
x values.
 ${\mathbf{ifail}}=5$

On entry, ${\mathbf{lda}}=\u2329\mathit{\text{value}}\u232a$ and ${\mathbf{kplus1}}=\u2329\mathit{\text{value}}\u232a$.
Constraint: ${\mathbf{lda}}\ge {\mathbf{kplus1}}$.
 ${\mathbf{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.
 ${\mathbf{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.
 ${\mathbf{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.
The time taken is approximately proportional to $m\left(k+1\right)\left(k+11\right)$.
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 equallyspaced data, this critical value is about
$2\times \sqrt{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 halfway 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,\dots ,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 selfstarting in that any number of datasets can be supplied.
10.1
Program Text
10.2
Program Data
10.3
Program Results