numapprox[remez] - Remez algorithm for minimax rational approximation
Calling Sequence
remez(w, f, a, b, m, n, crit, ‘maxerror‘)
Parameters
w
-
procedure representing a weight function w(x) > 0 on [a, b]
f
-
procedure representing the function f(x) to be approximated
a, b
-
numeric values specifying the interval [a, b]
m
-
integer specifying the desired degree of the numerator
n
-
integer specifying the desired degree of the denominator
crit
-
Array indexed
1 .. m + n + 2
containing an initial estimate of the critical set (i.e. the points of max/min of the error curve)
maxerror
-
name which will be assigned the minimax norm of
w |f - r|
Description
" This is not usually invoked as a user-level routine. See numapprox[minimax] for the standard user interface to the Remez algorithm.
" This procedure computes the best minimax rational approximation of degree
m, n
for a given real function f(x) on the interval [a, b] with respect to the positive weight function w(x).
" Specifically, it computes the rational expression r(x) such that
max(w(x) |f(x) - r(x)|, `in`(x, [a, b]))
(1)
is minimized over all rational expressions
p(x)
r(x) = ----
q(x)
with numerator of degree m and denominator of degree n.
" The value returned is an operator r such that
r(x)
is the desired approximation as a quotient of polynomials in Horner (nested multiplication) form.
" Note that if f(x) is nonzero on the interval of approximation then the relative error will be minimized by specifying the weight function
1
w(x) = ------
|f(x)|
.
" If
n = 0
then the best minimax polynomial approximation of degree m is computed.
" The last argument ‘maxerror‘ must be a name and upon return, its value will be an estimate of the minimax norm specified by equation (1) above.
" Various levels of user information will be displayed during the computation if infolevel[remez] is assigned values between 1 and 3.
" The command with(numapprox,remez) allows the use of the abbreviated form of this command.
Examples
with(numapprox);
w := proc(x) 1.0 end proc:
f := proc(x) evalf(exp(x)) end proc:
crit := Array(1 .. 7, [0, .10, .25, .50, .75, .90, 1.0]);
remez(w, f, 0, 1, 5, 0, crit, ‘maxerror‘);
x -> 0.9999988705 + (1.000079457 + (0.4990960985
+ (0.1704019738 + (0.03480057115 + 0.01390372811 x) x) x) x) x
maxerror;
0.000001129739955
Digits := 14;
Digits := 14
g := proc(x) if x=0 then 1.0 else evalf(tan(x)/x) end if end proc:
crit := Array(1 .. 8, [0, 0.5e-1, .15, .30, .48, .63, .73, .78]);
remez(w, g, 0, evalf((1/4)*Pi), 3, 3, crit, ‘maxerror‘);
x -> (1.2864938726745 + (-0.50393137136308
