RATCH

WEIGHT — User-supplied FUNCTION to scale the maximum error. It must be continuous and nonvanishing on the closed interval (A, B). The form is WEIGHT(X), where

X – Independent variable. (Input)
WEIGHT – The function value. (Output)

WEIGHT must be declared EXTERNAL in the calling program.

A — Lower end of the interval on which the approximation is desired. (Input)

B — Upper end of the interval on which the approximation is desired. (Input)

P — Vector of length N + 1 containing the coefficients of the numerator polynomial. (Output)

Q — Vector of length M + 1 containing the coefficients of the denominator polynomial. (Output)

ERROR — Min-max error of approximation. (Output)

Optional Arguments

N — The degree of the numerator. (Input)
Default: N = size (P,1) - 1.

M — The degree of the denominator. (Input)
Default: M = size (Q,1) - 1.

FORTRAN 90 Interface

Generic: CALL RATCH (F, PHI, WEIGHT, A, B, P, Q, ERROR [, …])

Specific: The specific interface names are S_RATCH and D_RATCH.

FORTRAN 77 Interface

Single: CALL RATCH (F, PHI, WEIGHT, A, B, N, M, P, Q, ERROR)

Double: The double precision name is DRATCH.

Description

The routine RATCH is designed to compute the best weighted L∞ (Chebyshev) approximant to a given function. Specifically, given a weight function w = WEIGHT, a monotone function ɸ = PHI, and a function f to be approximated on the interval [a, b], the subroutine RATCH returns the coefficients (in P and Q) for a rational approximation to f on [a, b]. The user must supply the degree of the numerator N and the degree of the denominator M of the rational function

The goal is to produce coefficients which minimize the expression

Notice that setting ɸ(x) = x yields ordinary rational approximation. A typical use of the function ɸ occurs when one wants to approximate an even function on a symmetric interval, say [−a, a] using ordinary rational functions. In this case, it is known that the answer must be an even function. Hence, one can set ɸ(x) = x2, only approximate on [0, a], and decrease by one half the degrees in the numerator and denominator.

The algorithm implemented in this subroutine is designed for fast execution. It assumes that the best approximant has precisely N + M + 2 equi-oscillations. That is, that there exist N + M + 2 points t1 < … < tN+M+2 satisfying

Such points are called alternants. Unfortunately, there are many instances in which the best rational approximant to the given function has either fewer alternants or more alternants. In this case, it is not expected that this subroutine will perform well. For more information on rational Chebyshev approximation, the reader can consult Cheney (1966). The subroutine is based on work of Cody, Fraser, and Hart (1968).

Comments

1. Workspace may be explicitly provided, if desired, by use of R2TCH/DR2TCH. The reference is:

CALL R2TCH (F, PHI, WEIGHT, A, B, N, M, P, Q, ERROR, ITMAX, IWK, WK)

The additional arguments are as follows:

ITMAX — Maximum number of iterations. (Input)
The default value is 20.

IWK — Workspace vector of length (N + M + 2). (Workspace)

WK — Workspace vector of length (N + M + 8) * (N + M + 2). (Workspace)

2. Informational errors

Type	Code	Description
3	1	The maximum number of iterations has been reached. The routine R2TCH may be called directly to set a larger value for ITMAX.
3	2	The error was reduced as far as numerically possible. A good approximation is returned in P and Q, but this does not necessarily give the Chebyshev approximation.
4	3	The linear system that defines P and Q was found to be algorithmically singular. This indicates the possibility of a degenerate approximation.
4	4	A sequence of critical points that was not monotonic generated. This indicates the possibility of a degenerate approximation.
4	5	The value of the error curve at some critical point is too large. This indicates the possibility of poles in the rational function.
4	6	The weight function cannot be zero on the closed interval (A, B).

Example

In this example, we compute the best rational approximation to the gamma function, Γ, on the interval [2, 3] with weight function w = 1 and N = M = 2. We display the maximum error and the coefficients. This problem is taken from the paper of Cody, Fraser, and Hart (1968). We compute in double precision due to the conditioning of this problem.

USE RATCH_INT

USE UMACH_INT

IMPLICIT NONE

INTEGER M, N

PARAMETER (M=2, N=2)

INTEGER NOUT

DOUBLE PRECISION A, B, ERROR, F, P(N+1), PHI, Q(M+1), WEIGHT

EXTERNAL F, PHI, WEIGHT

A = 2.0D0

B = 3.0D0

! Compute double precision rational

! approximation

CALL RATCH (F, PHI, WEIGHT, A, B, P, Q, ERROR)

! Get output unit number

CALL UMACH (2, NOUT)

! Print P, Q and min-max error

WRITE (NOUT,'(1X,A)') 'In double precision we have:'

WRITE (NOUT,99999) 'P = ', P

WRITE (NOUT,99999) 'Q = ', Q

WRITE (NOUT,99999) 'ERROR = ', ERROR

99999 FORMAT (' ', A, 5X, 3F20.12, /)

END

! -----------------------------------------------------------------------

DOUBLE PRECISION FUNCTION F (X)

DOUBLE PRECISION X

DOUBLE PRECISION DGAMMA

EXTERNAL DGAMMA

F = DGAMMA(X)

RETURN

END

! -----------------------------------------------------------------------

DOUBLE PRECISION FUNCTION PHI (X)

DOUBLE PRECISION X

PHI = X

RETURN

END

! -----------------------------------------------------------------------

DOUBLE PRECISION FUNCTION WEIGHT (X)

DOUBLE PRECISION X

DOUBLE PRECISION DGAMMA

EXTERNAL DGAMMA

WEIGHT = DGAMMA(X)

RETURN

END

Output

In double precision we have:

P = 1.265583562487 -0.650585004466 0.197868699191

Q = 1.000000000000 -0.064342721236 -0.028851461855

ERROR = -0.000026934190