CHEBGRID Name
CHEBGRID
Author
Craig B. Markwardt, NASA/GSFC Code 662, Greenbelt, MD 20770
craigm@lheamail.gsfc.nasa.gov
UPDATED VERSIONs can be found on my WEB PAGE:
http://cow.physics.wisc.edu/~craigm/idl/idl.html
Purpose
Estimate Chebyshev polynomial coefficients of a function on a grid
Major Topics
Curve and Surface Fitting
Calling Sequence
p = CHEBGRID(T, X, [ DXDT, NPOINTS=, NPOLY=, NGRANULE= , $
RMS=, DRMS=, RESIDUALS=, DRESIDUALS= , $
XMATRIX=, DXMATRIX=, RESET=,
DERIV_WEIGHT= ] ) Description
CHEBGRID estimates the coefficients for a finite sum of Chebyshev
polynomials approximating a continuous tabulated function over an
interval. The function (and optionally its derivative) must be
tabulated on a regularly sampled grid. The implementation of this
function is taken from a method described by X. X. Newhall, used
in estimating coefficients for ephemerides in the solar system.
The tabulated function is assumed to be continuous over the entire
interval. A Chebyshev series is fitted to the function over small
segments, called granules. The size of each granule, the number
of points in each granule, and the number of Chebyshev polynomials
are all configurable.
Users may specify either the function alone, or the function and
its first derivative. By also giving the tabulated derivative, a
more accurate Chebyshev polynomial can be developed. Aside from
the constraints mentioned in the next paragraph, the polynomial
that is returned is the best-fit polynomial in a least-squares
sense.
Here is a definition of terms:
GRANULE - a single continuous fitted segment. The length of the
granule, NGRANULE, is specified in units of the tabulated
grid size. Because of the continuity requirements developed
below, granules will always overlap at their endpoints.
Thus, then length of a granule should be a factor of
N_ELEMENTS(X)-1. For simple functions over short intervals,
the granule size can be equal to N_ELEMENTS(X)-1
NUMBER OF POINTS the number of points, NPOINTS, within a
granule to be fitted to the polynomial, not necessarily
equal to the granule size. The greater the number of
points, the more computation time and storage is required.
This number *must* be a factor of NGRANULE. Typically
NPOINTS is a number between 8 and 12. Because of the
single-point overlap between granules (see below), the
actual number of points per fit is NPOINTS+1.
NUMBER OF POLYNOMIALS the number of Chebyshev polynomial terms,
NPOLYNOMIAL, to be fitted per granule. The greater the
number of polynomial terms, the more computation time and
storage is required, but also the greater the approximating
precision of the fit.
The particular set of Chebyshev polynomial coefficients developed
by this function have some special properties. If both the
function and its derivative are specified, then the value and
derivative of the interpolating polynomial at the granule
endpoints will be exactly equal to the tabulated endpoint values.
This feature allows many approximations to be strung together
piecewise, and the function value and first derivative will be
continuous across granule boundaries.
If only the function value is specified, then only the function
value will be continuous at the granule endpoints, and not the
derivative.
An extensive set of statistics are computed to assess the quality
of the Chebyshev polynomial fit. The keywords RESIDUALS and
DRESIDUALS return the residuals of the fit after subtracting the
interpolation. The RMS and DRMS keywords return the root mean
squared deviations between data and model.
If the user does not know how many granules, points, or polynomial
coefficients to use, then he or she should try several
combinations and see which minimizes the r.m.s. value with the
fewest number of coefficients.
If the XMATRIX and DXMATRIX keywords are passed, then CHEBGRID
attempts to avoid recomputing several of the matrices it uses in
estimating the coefficients. If multiple calls to CHEBGRID are to
be made, some compution time savings can be made. In the first
call CHEBGRID the required matrices are computed and returned. In
subsequent calls, CHEBGRID detects the XMATRIX and DXMATRIX
keyword values and uses those values if it can.
The user can also estimate their own coefficients. The matrices
returned are (NPOINTS+1)x(NPOLYNOMIAL). The coefficients from a
NPOINTS+1 tabulation, X, are found by:
PCHEB = XMATRIX ## X + DXMATRIX ## DXDT
if derivative information is known, or
PCHEB = XMATRIX ## X
if no derivative information is known. [ Note: the matrices are
different, depending on whether derivative information is known or
not. ]
Inputs
T - array of regularly sampled *independent* variables. The number
of elements in T should be a multiple of NGRANULE, plus one.
X - array of regularly sampled *dependent* variables. The number
of elements in X should be equal to the number of elements in
T.
DXDT - optionally, a tabulated array of first derivatives of X
with respect to T, at the same grid points.
Keyword Parameters
NGRANULE - size of a "granule", in grid intervals. NGRANULE must
be at least 2, and a factor of N_ELEMENTS(T)-1.
Default: 8
NPOINTS - number of points per granule that are fitted. NPOINTS
must be at least 2, and a factor of NGRANULE.
Default: NGRANULE
NPOLYNOMIAL - number of Chebyshev polynomial terms per fit.
NPOLYNOMIAL must be at least 2 and less than
2*(NPOINTS+1), when derivative information is
specified; or less than NPOINTS+1, when no
derivative information is specified.
Default: 7
RESIDUALS - upon return, an array of size N_ELEMENTS(T), with
residuals of the tabulated function minus the
interpolated function.
DRESIDUALS - same as RESIDUALS, but for the function's first
derivative.
RMS - upon return, the root mean square of the function value
residuals.
DRMS - same as RMS, but for the function's first derivative.
XMATRIX - upon return, the matrix used to compute Chebyshev
polynomial coefficients from the function value.
Upon input, CHEBGRID determines if XMATRIX will apply to
the data, and if so, XMATRIX is reused rather than
computed. If XMATRIX cannot be reused, then it is
computed afresh, and the new value is returned in the
XMATRIX keyword.
The user should not modify the contents of this array.
DXMATRIX - same as XMATRIX, but for the function's first
derivative.
RESET - if set, force a recomputation of XMATRIX and/or DXMATRIX.
DERIV_WEIGHT - amount of weight to give to function derivative,
relative to the function value.
Default: 0.16d Returns
An array of coefficient values. The dimensions of the array are
NPOLYNOMIALxNSEGS, where NSEGS is the number of granules in the
entire interval. Example
;; Estimate Chebyshev coefficients for the function SIN(X), on the
;; interval [-1,+1].
xx = dindgen(9)/4d - 1d ;; Regular grid from -1 to 1 (9 points)
yy = sin(xx) ;; Function values, sin(x), ...
dy = cos(xx) ;; ... and derivatives
;; Estimate coefficients using CHEBGRID (single granule of 8 intervals)
p = chebgrid(xx, yy, dy, npoints=8, ngranule=8, npoly=10)
xxx = dindgen(1001)/500 - 1d ;; New grid for testing
res = sin(xxx) - chebeval(xxx, p)
plot, xxx, res
;; Same as example above, except extended range to [-1, +15],
using eight granules.
xx2 = dindgen(65)/4d - 1
yy2 = sin(xx2)
dy2 = cos(xx2)
p = chebgrid(xx2, yy2, dy2, ngranule=8, npoint=8, npoly=10)
help, p
P DOUBLE = Array[10, 8]
;; (i.e., 10 polynomial coefficients over 8 granules)
References
Abramowitz, M. & Stegun, I., 1965, *Handbook of Mathematical
Functions*, 1965, U.S. Government Printing Office, Washington,
D.C. (Applied Mathematical Series 55)
Newhall, X. X. 1989, Celestial Mechanics, 45, p. 305-310
Modification History
Written, CM, Feb 2002
Documented, CM, 24 Mar 2002
Corrected documentation, CM, 28 Apr 2002
Typo correction, CM, 10 Oct 2002