hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_fit_1dcheb_con (e02ag)

Purpose

nag_fit_1dcheb_con (e02ag) 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.

Syntax

[a, s, np1, wrk, ifail] = e02ag(k, xmin, xmax, x, y, w, xf, yf, ip, lwrk, 'm', m, 'mf', mf)
[a, s, np1, wrk, ifail] = nag_fit_1dcheb_con(k, xmin, xmax, x, y, w, xf, yf, ip, lwrk, 'm', m, 'mf', mf)

Description

nag_fit_1dcheb_con (e02ag) determines least squares polynomial approximations of degrees up to kk to the set of data points (xr,yr)(xr,yr) with weights wrwr, for r = 1,2,,mr=1,2,,m. The value of kk, the maximum degree required, is to be prescribed by you. At each of the values xfrxfr, for r = 1,2,,mfr=1,2,,mf, of the independent variable xx, the approximations and their derivatives up to order prpr are constrained to have one of the values yfsyfs, for s = 1,2,,ns=1,2,,n, specified by you, where n = mf + r = 0mfprn=mf+r=0mfpr.
The approximation of degree ii has the property that, subject to the imposed constraints, it minimizes σiσi, the sum of the squares of the weighted residuals εrεr, for r = 1,2,,mr=1,2,,m, where
εr = wr(yrfi(xr))
εr=wr(yr-fi(xr))
and fi(xr)fi(xr) is the value of the polynomial approximation of degree ii at the rrth data point.
Each polynomial is represented in Chebyshev series form with normalized argument xx-. This argument lies in the range 1-1 to + 1+1 and is related to the original variable xx by the linear transformation
x = (2x(xmax + xmin))/((xmaxxmin))
x-=2x-(xmax+xmin) (xmax-xmin)
where xminxmin and xmaxxmax, specified by you, are respectively the lower and upper end points of the interval of xx over which the polynomials are to be defined.
The polynomial approximation of degree ii can be written as
(1/2)ai,0 + ai,1T1(x) + + aijTj(x) + + aiiTi(x)
12ai,0+ai,1T1(x-)++aijTj(x-)++aiiTi(x-)
where Tj(x)Tj(x-) is the Chebyshev polynomial of the first kind of degree jj with argument xx-. For i = n,n + 1,,ki=n,n+1,,k, the function produces the values of the coefficients aijaij, for j = 0,1,,ij=0,1,,i, together with the value of the root mean square residual,
Si = sqrt( ( σi )/((m + ni1)) ),
Si = σ i ( m +n -i -1 ) ,
where mm is the number of data points with nonzero weight.
Values of the approximations may subsequently be computed using nag_fit_1dcheb_eval (e02ae) or nag_fit_1dcheb_eval2 (e02ak).
First nag_fit_1dcheb_con (e02ag) determines a polynomial μ(x)μ(x-), of degree n1n-1, which satisfies the given constraints, and a polynomial ν(x)ν(x-), of degree nn, which has value (or derivative) zero wherever a constrained value (or derivative) is specified. It then fits yrμ(xr)yr-μ(xr), for r = 1,2,,mr=1,2,,m, with polynomials of the required degree in xx- each with factor ν(x)ν(x-). Finally the coefficients of μ(x)μ(x-) are added to the coefficients of these fits to give the coefficients of the constrained polynomial approximations to the data points (xr,yr)(xr,yr), for r = 1,2,,mr=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)).

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

Parameters

Compulsory Input Parameters

1:     k – int64int32nag_int scalar
kk, the maximum degree required.
Constraint: nkm + n1nkm+n-1 where nn is the total number of constraints and mm is the number of data points with nonzero weights and distinct abscissae which do not coincide with any of the xfrxfr.
2:     xmin – double scalar
3:     xmax – double scalar
The lower and upper end points, respectively, of the interval [xmin,xmax][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 xrxr and xfrxfr. This avoids the danger of extrapolation provided there is a constraint point or data point with nonzero weight at each end point.
Constraint: xmax > xminxmax>xmin.
4:     x(m) – double array
m, the dimension of the array, must satisfy the constraint m1m1.
x(r)xr must contain the value xrxr of the independent variable at the rrth data point, for r = 1,2,,mr=1,2,,m.
Constraint: the x(r)xr must be in nondecreasing order and satisfy xminx(r)xmaxxminxrxmax.
5:     y(m) – double array
m, the dimension of the array, must satisfy the constraint m1m1.
y(r)yr must contain yryr, the value of the dependent variable at the rrth data point, for r = 1,2,,mr=1,2,,m.
6:     w(m) – double array
m, the dimension of the array, must satisfy the constraint m1m1.
w(r)wr must contain the weight wrwr to be applied to the data point xrxr, for r = 1,2,,mr=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 xx and yy values both coincide with those of a constraint (otherwise the denominators involved in the root mean square residuals SiSi will be slightly in error).
7:     xf(mf) – double array
mf, the dimension of the array, must satisfy the constraint mf1mf1.
xf(r)xfr must contain xfrxfr, the value of the independent variable at which a constraint is specified, for r = 1,2,,mfr=1,2,,mf.
Constraint: these values need not be ordered but must be distinct and satisfy xminxf(r)xmaxxminxfrxmax.
8:     yf(lyf) – double array
lyf, the dimension of the array, must satisfy the constraint lyfmf + i = 1mf ip(i)lyfmf+ i=1 mf ipi.
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)xfr, yf contains in successive elements the required value of the approximation, its first derivative, second derivative, ,pr,prth derivative, for r = 1,2,,mfr=1,2,,mf. Thus the value, yfsyfs, which the kkth derivative of each approximation (k = 0k=0 referring to the approximation itself) is required to take at the point xf(r)xfr must be contained in yf(s)yfs, where
s = r + k + p1 + p2 + + pr1,
s=r+k+p1+p2++pr-1,
where k = 0,1,,prk=0,1,,pr and r = 1,2,,mfr=1,2,,mf. The derivatives are with respect to the independent variable xx.
9:     ip(mf) – int64int32nag_int array
mf, the dimension of the array, must satisfy the constraint mf1mf1.
ip(r)ipr must contain prpr, the order of the highest-order derivative specified at xf(r)xfr, for r = 1,2,,mfr=1,2,,mf. pr = 0pr=0 implies that the value of the approximation at xf(r)xfr is specified, but not that of any derivative.
Constraint: ip(r)0ipr0, for r = 1,2,,mfr=1,2,,mf.
10:   lwrk – int64int32nag_int scalar
The dimension of the array wrk as declared in the (sub)program from which nag_fit_1dcheb_con (e02ag) is called.
Constraint: lwrkmax (4 × m + 3 × kplus1,8 × n + 5 × ipmax + mf + 10) + 2 × n + 2lwrkmax(4×m+3×kplus1,8×n+5×ipmax+mf+10)+2×n+2, where ipmax = max (ip(r))ipmax=max(ipr), for r = 1,2,,mfr=1,2,,mf.

Optional Input Parameters

1:     m – int64int32nag_int scalar
Default: The dimension of the arrays x, y, w. (An error is raised if these dimensions are not equal.)
mm, the number of data points to be fitted.
Constraint: m1m1.
2:     mf – int64int32nag_int scalar
Default: The dimension of the arrays xf, ip. (An error is raised if these dimensions are not equal.)
mfmf, the number of values of the independent variable at which a constraint is specified.
Constraint: mf1mf1.

Input Parameters Omitted from the MATLAB Interface

kplus1 lda lyf np1 iwrk liwrk

Output Parameters

1:     a(lda,kplus1) – double array
kplus1 = k + 1kplus1=k+1.
ldakplus1ldakplus1.
a(i + 1,j + 1)ai+1j+1 contains the coefficient aijaij in the approximating polynomial of degree ii, for i = n,,ki=n,,k and j = 0,1,,ij=0,1,,i.
2:     s(kplus1) – double array
kplus1 = k + 1kplus1=k+1.
s(i + 1)si+1 contains SiSi, for i = n,,ki=n,,k, the root mean square residual corresponding to the approximating polynomial of degree ii. In the case where the number of data points with nonzero weight is equal to k + 1nk+1-n, SiSi is indeterminate: the function sets it to zero. For the interpretation of the values of SiSi and their use in selecting an appropriate degree, see Section [General] in the E02 Chapter Introduction.
3:     np1 – int64int32nag_int scalar
n + 1n+1, where nn is the total number of constraint conditions imposed: n = mf + p1 + p2 + + pmfn=mf+p1+p2++pmf.
4:     wrk(lwrk) – double array
Contains weighted residuals of the highest degree of fit determined (k)(k). The residual at xrxr is in element 2(n + 1) + 3(m + k + 1) + r2(n+1)+3(m+k+1)+r, for r = 1,2,,mr=1,2,,m. The rest of the array is used as workspace.
5:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry,m < 1m<1,
orkplus1 < n + 1kplus1<n+1,
orlda < kplus1lda<kplus1,
ormf < 1mf<1,
orlyf < nlyf<n,
orlwrk is too small (see Section [Parameters]),
orliwrk < 2 × mf + 2liwrk<2×mf+2.
(Here nn is the total number of constraint conditions.)
  ifail = 2ifail=2
ip(r) < 0ipr<0 for some r = 1,2,,mfr=1,2,,mf.
  ifail = 3ifail=3
xminxmaxxminxmax, or xf(r)xfr is not in the interval xmin to xmax for some r = 1,2,,mfr=1,2,,mf, or the xf(r)xfr are not distinct.
  ifail = 4ifail=4
x(r)xr is not in the interval xmin to xmax for some r = 1,2,,mr=1,2,,m.
  ifail = 5ifail=5
x(r) < x(r1)xr<xr-1 for some r = 2,3,,mr=2,3,,m.
  ifail = 6ifail=6
kplus1 > m + nkplus1>m+n, where mm is the number of data points with nonzero weight and distinct abscissae which do not coincide with any xf(r)xfr. Thus there is no unique solution.
  ifail = 7ifail=7
The polynomials μ(x)μ(x) and/or ν(x)ν(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.

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 function.

Further Comments

The time taken to form the interpolating polynomial is approximately proportional to n3n3, and that to form the approximating polynomials is very approximately proportional to m(k + 1)(k + 1n)m(k+1)(k+1-n).
To carry out a least squares polynomial fit without constraints, use nag_fit_1dcheb_arb (e02ad). To carry out polynomial interpolation only, use nag_interp_1d_cheb (e01ae).

Example

function nag_fit_1dcheb_con_example
k = int64(4);
xmin = 0;
xmax = 4;
x = [0.5;
     1;
     2;
     2.5;
     3];
y = [0.03;
     -0.75;
     -1;
     -0.1;
     1.75];
w = [1;
     1;
     1;
     1;
     1];
xf = [0;
     4];
yf = [1;
     -2;
     9];
ip = [int64(1);0];
lwrk = int64(200);
[a, s, np1, wrk, ifail] = nag_fit_1dcheb_con(k, xmin, xmax, x, y, w, xf, yf, ip, lwrk)
 

a =

         0         0         0         0         0
         0         0         0         0         0
         0         0         0         0         0
    3.9980    3.4995    3.0010    0.5005         0
    3.9980    3.4995    3.0010    0.5005   -0.0000


s =

         0
         0
         0
    0.0025
    0.0029


np1 =

                    4


wrk =

    6.0000
    4.0000
    2.0000
         0
   -1.0000
   -0.2500
    0.5000
    0.2500
   -0.2200
   -0.7500
   -2.0000
   -2.3500
   -2.2500
   -0.7500
    0.8632
   -0.5000
    2.2096
         0
    1.8923
    0.2500
   -0.1262
    0.5000
   -2.3711
   -2.0776
   -1.5269
   -0.8471
   -0.1077
   -0.5131
   -1.5269
   -0.1529
   -0.9462
    1.0000
    1.5269
    1.9443
    2.4006
    2.8146
    2.4006
    3.8531
   -0.0010
    0.0008
    0.0020
   -0.0039
    0.0022
         0
    0.2500
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0


function e02ag_example
k = int64(4);
xmin = 0;
xmax = 4;
x = [0.5;
     1;
     2;
     2.5;
     3];
y = [0.03;
     -0.75;
     -1;
     -0.1;
     1.75];
w = [1;
     1;
     1;
     1;
     1];
xf = [0;
     4];
yf = [1;
     -2;
     9];
ip = [int64(1);0];
lwrk = int64(200);
[a, s, np1, wrk, ifail] = e02ag(k, xmin, xmax, x, y, w, xf, yf, ip, lwrk)
 

a =

         0         0         0         0         0
         0         0         0         0         0
         0         0         0         0         0
    3.9980    3.4995    3.0010    0.5005         0
    3.9980    3.4995    3.0010    0.5005   -0.0000


s =

         0
         0
         0
    0.0025
    0.0029


np1 =

                    4


wrk =

    6.0000
    4.0000
    2.0000
         0
   -1.0000
   -0.2500
    0.5000
    0.2500
   -0.2200
   -0.7500
   -2.0000
   -2.3500
   -2.2500
   -0.7500
    0.8632
   -0.5000
    2.2096
         0
    1.8923
    0.2500
   -0.1262
    0.5000
   -2.3711
   -2.0776
   -1.5269
   -0.8471
   -0.1077
   -0.5131
   -1.5269
   -0.1529
   -0.9462
    1.0000
    1.5269
    1.9443
    2.4006
    2.8146
    2.4006
    3.8531
   -0.0010
    0.0008
    0.0020
   -0.0039
    0.0022
         0
    0.2500
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0



PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013