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

Fortran Interface

Subroutine e02agf (

m, kplus1, lda, xmin, xmax, x, y, w, mf, xf, yf, lyf, ip, a, s, np1, wrk, lwrk, iwrk, liwrk, ifail)

Integer, Intent (In)	::	m, kplus1, lda, mf, lyf, ip(mf), lwrk, liwrk
Integer, Intent (Inout)	::	ifail
Integer, Intent (Out)	::	np1, iwrk(liwrk)
Real (Kind=nag_wp), Intent (In)	::	xmin, xmax, x(m), y(m), w(m), xf(mf), yf(lyf)
Real (Kind=nag_wp), Intent (Inout)	::	a(lda,kplus1)
Real (Kind=nag_wp), Intent (Out)	::	s(kplus1), wrk(lwrk)

C Header Interface

#include <nag.h>

void

e02agf_ (const Integer *m, const Integer *kplus1, const Integer *lda, const double *xmin, const double *xmax, const double x[], const double y[], const double w[], const Integer *mf, const double xf[], const double yf[], const Integer *lyf, const Integer ip[], double a[], double s[], Integer *np1, double wrk[], const Integer *lwrk, Integer iwrk[], const Integer *liwrk, Integer *ifail)

The routine may be called by the names e02agf or nagf_fit_dim1_cheb_con.

3 Description

e02agf determines least squares polynomial approximations of degrees up to

k

to the set of data points

(x_{r}, y_{r})

with weights

w_{r}

, for

r = 1, 2, \dots, m

. The value of

k

, the maximum degree required, is to be prescribed by you. At each of the values

{x f}_{r}

, for

r = 1, 2, \dots, m f

, of the independent variable

x

, the approximations and their derivatives up to order

p_{r}

are constrained to have one of the values

{y f}_{s}

, for

s = 1, 2, \dots, n

, specified by you, where

n = m f + \sum_{r = 0}^{m f} p_{r}

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, \dots, m

, where

ε_{r} = w_{r} (y_{r} - f_{i} (x_{r}))

and

f_{i} (x_{r})

is the value of the polynomial approximation of degree

i

at the

r

th data point.

Each polynomial is represented in Chebyshev series form with normalized argument

\bar{x}

. This argument lies in the range

−1

+ 1

and is related to the original variable

x

by the linear transformation

\bar{x} = \frac{2 x - (x_{\max} + x_{\min})}{(x_{\max} - x_{\min})}

where

x_{\min}

and

x_{\max}

, 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

\frac{1}{2} a_{i, 0} + a_{i, 1} T_{1} (\bar{x}) + \dots + a_{i j} T_{j} (\bar{x}) + \dots + a_{i i} T_{i} (\bar{x})

where

T_{j} (\bar{x})

is the Chebyshev polynomial of the first kind of degree

j

with argument

\bar{x}

. For

i = n, n + 1, \dots, k

, the routine produces the values of the coefficients

a_{i j}

, for

j = 0, 1, \dots, i

, together with the value of the root mean square residual,

S_{i} = \sqrt{\frac{σ_{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

μ (\bar{x})

, of degree

n - 1

, which satisfies the given constraints, and a polynomial

ν (\bar{x})

, of degree

n

, which has value (or derivative) zero wherever a constrained value (or derivative) is specified. It then fits

y_{r} - μ (x_{r})

, for

r = 1, 2, \dots, m

, with polynomials of the required degree in

\bar{x}

each with factor

ν (\bar{x})

. Finally the coefficients of

μ (\bar{x})

are added to the coefficients of these fits to give the coefficients of the constrained polynomial approximations to the data points

(x_{r}, y_{r})

, for

r = 1, 2, \dots, 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 Arguments

1: $m$ – Integer Input

On entry:

m

, the number of data points to be fitted.

Constraint:

m \geq 1

2: $kplus1$ – Integer Input

On entry:

k + 1

, where

k

is the maximum degree required.

Constraint:

n + 1 \leq kplus1 \leq m^{''} + 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

xf (r)

3: $lda$ – Integer Input

On entry: the first dimension of the array a as declared in the (sub)program from which e02agf is called.

Constraint:

lda \geq kplus1

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

[x_{\min}, x_{\max}]

. 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

x_{r}

and

{x f}_{r}

. 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) array Input

On entry:

x (r)

must contain the value

x_{r}

of the independent variable at the

r

th data point, for

r = 1, 2, \dots, m

Constraint: the

x (r)

must be in nondecreasing order and satisfy

xmin \leq x (r) \leq xmax

7: $y (m)$ – Real (Kind=nag_wp) array Input

On entry:

y (r)

must contain

y_{r}

, the value of the dependent variable at the

r

th data point, for

r = 1, 2, \dots, m

8: $w (m)$ – Real (Kind=nag_wp) array Input

On entry:

w (r)

must contain the weight

w_{r}

to be applied to the data point

x_{r}

, for

r = 1, 2, \dots, 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

S_{i}

will be slightly in error).

9: $mf$ – Integer Input

On entry:

m f

, the number of values of the independent variable at which a constraint is specified.

Constraint:

mf \geq 1

10: $xf (mf)$ – Real (Kind=nag_wp) array Input

On entry:

xf (r)

must contain

{x f}_{r}

, the value of the independent variable at which a constraint is specified, for

r = 1, 2, \dots, mf

Constraint: these values need not be ordered but must be distinct and satisfy

xmin \leq xf (r) \leq xmax

11: $yf (lyf)$ – Real (Kind=nag_wp) array Input

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

xf (r)

, yf contains in successive elements the required value of the approximation, its first derivative, second derivative,

\dots, p_{r}

th derivative, for

r = 1, 2, \dots, m f

. Thus the value,

{y f}_{s}

, which the

k

th derivative of each approximation (

k = 0

referring to the approximation itself) is required to take at the point

xf (r)

must be contained in

yf (s)

, where

s = r + k + p_{1} + p_{2} + \dots + p_{r - 1},

where

k = 0, 1, \dots, p_{r}

and

r = 1, 2, \dots, m f

. The derivatives are with respect to the independent variable

x

12: $lyf$ – Integer Input

On entry: the dimension of the array yf as declared in the (sub)program from which e02agf is called.

Constraint:

lyf \geq mf + \sum_{i = 1}^{mf} ip (i)

13: $ip (mf)$ – Integer array Input

On entry:

ip (r)

must contain

p_{r}

, the order of the highest-order derivative specified at

xf (r)

, for

r = 1, 2, \dots, m f

p_{r} = 0

implies that the value of the approximation at

xf (r)

is specified, but not that of any derivative.

Constraint:

ip (r) \geq 0

, for

r = 1, 2, \dots, mf

14: $a (lda, kplus1)$ – Real (Kind=nag_wp) array Output

On exit:

a (i + 1, j + 1)

contains the coefficient

a_{i j}

in the approximating polynomial of degree

i

, for

i = n, \dots, k

and

j = 0, 1, \dots, i

15: $s (kplus1)$ – Real (Kind=nag_wp) array Output

On exit:

s (i + 1)

contains

S_{i}

, for

i = n, \dots, 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

S_{i}

is indeterminate: the routine sets it to zero. For the interpretation of the values of

S_{i}

and their use in selecting an appropriate degree, see Section 3.1 in the E02 Chapter Introduction.

16: $np1$ – Integer Output

On exit:

n + 1

, where

n

is the total number of constraint conditions imposed:

n = m f + p_{1} + p_{2} + \dots + p_{m f}

17: $wrk (lwrk)$ – Real (Kind=nag_wp) array Output

On exit: contains weighted residuals of the highest degree of fit determined

(k)

. The residual at

x_{r}

is in element

2 (n + 1) + 3 (m + k + 1) + r

, for

r = 1, 2, \dots, m

. The rest of the array is used as workspace.

18: $lwrk$ – Integer Input

On entry: the dimension of the array wrk as declared in the (sub)program from which e02agf is called.

Constraint:

lwrk \geq \max (4 \times m + 3 \times kplus1, 8 \times n + 5 \times ipmax + mf + 10) + 2 \times n + 2

, where

ipmax = \max (ip (r))

, for

r = 1, 2, \dots, m f

19: $iwrk (liwrk)$ – Integer array Workspace

20: $liwrk$ – Integer Input

On entry: the dimension of the array iwrk as declared in the (sub)program from which e02agf is called.

Constraint:

liwrk \geq 2 \times mf + 2

21: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

−1

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

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

−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, $kplus1 = ⟨ value ⟩$ and $N = ⟨ value ⟩$ .
Constraint: $kplus1 \geq N + 1$ , where N is the total number of constraints.

On entry, $lda = ⟨ value ⟩$ and $kplus1 = ⟨ value ⟩$ .
Constraint: $lda \geq kplus1$ .

On entry, liwrk is too small. $liwrk = ⟨ value ⟩$ . Minimum possible dimension: $⟨ value ⟩$ .

On entry, lwrk is too small. $lwrk = ⟨ value ⟩$ . Minimum possible dimension: $⟨ value ⟩$ .

On entry, $lyf = ⟨ value ⟩$ and $N = ⟨ value ⟩$ .
Constraint: $lyf \geq N$ , where N is the total number of constraints.

On entry, $m = ⟨ value ⟩$ .
Constraint: $m \geq 1$ .

On entry, $mf = ⟨ value ⟩$ .
Constraint: $mf \geq 1$ .

$ifail = 2$: On entry, $I = ⟨ value ⟩$ and $ip (I) = ⟨ value ⟩$ .
Constraint: $ip (I) \geq 0$ .

$ifail = 3$: On entry, $I = ⟨ value ⟩$ , $xf (I) = ⟨ value ⟩$ , $J = ⟨ value ⟩$ and $xf (J) = ⟨ value ⟩$ .
Constraint: $xf (I) \neq xf (J)$ .

On entry, $xf (I)$ lies outside interval $[xmin, xmax]$ : $I = ⟨ value ⟩$ , $xf (I) = ⟨ value ⟩$ , $xmin = ⟨ value ⟩$ and $xmax = ⟨ value ⟩$ .

On entry, $xmin = ⟨ value ⟩$ and $xmax = ⟨ value ⟩$ .
Constraint: $xmin < xmax$ .

$ifail = 4$: On entry, $x (I)$ lies outside interval $[xmin, xmax]$ : $I = ⟨ value ⟩$ , $x (I) = ⟨ value ⟩$ , $xmin = ⟨ value ⟩$ and $xmax = ⟨ value ⟩$ .

On entry, $x (I)$ lies outside interval $[xmin, xmax]$ for some $I$ .

$ifail = 5$: On entry, $i = ⟨ value ⟩$ , $x (i) = ⟨ value ⟩$ and $x (i - 1) = ⟨ value ⟩$ .
Constraint: $x (i) \geq x (i - 1)$ .

$ifail = 6$: On entry, $kplus1 > mdist + N$ , where mdist is the number of data points with nonzero weight and distinct abscissae different from all xf, and N is the total number of constraints: $kplus1 = ⟨ value ⟩$ , $mdist = ⟨ value ⟩$ and $N = ⟨ value ⟩$ .

$ifail = 7$: The polynomials $μ (x)$ and/or $ν (x)$ cannot be found. The problem 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.

$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 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 Parallelism and Performance

Background information to multithreading can be found in the Multithreading documentation.

e02agf is not threaded in any implementation.

9 Further Comments

The time taken to form the interpolating polynomial is approximately proportional to

n^{3}

, and that to form the approximating polynomials is very approximately proportional to

m (k + 1) (k + 1 - n)

To carry out a least squares polynomial fit without constraints, use e02adf. To carry out polynomial interpolation only, use e01aef.

10 Example

This example reads data in the following order, using the notation of the argument list above:

mf
$ip (i)$ , $xf (i)$ , Y-value and derivative values (if any) at $xf (i)$ , for $i = 1, 2, \dots, mf$
m
$x (i)$ , $y (i)$ , $w (i)$ , for $i = 1, 2, \dots, m$
$k$ , xmin, xmax

The output is:

the root mean square residual for each degree from $n$ to $k$ ;
the Chebyshev coefficients for the fit of degree $k$ ;
the data points, and the fitted values and residuals for the fit of degree $k$ .

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:

$p (0.0) = 1.0, p^{'} (0.0) = - 2.0, p (4.0) = 9.0 .$

e02ag: FL CL CPP AD PY MB

NAG FL Interfacee02agf (dim1_​cheb_​con)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

NAG FL Interface
e02agf (dim1_cheb_con)