function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
Best polynomial approximation in discrete uniform norm
INPUT VARIABLES:
M : degree of the fitting polynomial
N : number of data points
X(N) : x-coordinates of data points
Y(N) : y-coordinates of data points
K : character of the polynomial:
K = 0 : mixed parity polynomial
K = 1 : odd polynomial ( X(1) must be > 0 )
K = 2 : even polynomial ( X(1) must be >= 0 )
EPSH : tolerance for leveling. A useful value for 24-bit
mantissa is EPSH = 2.0E-7
MAXIT : upper limit for number of exchange steps
REF0(M2): initial alternating set ( N-vector ). This is an
OPTIONAL argument. The length M2 is given by:
M2 = M + 2 , if K = 0
M2 = integer part of (M+3)/2 , if K = 1
M2 = 2 + M/2 (M must be even) , if K = 2
OUTPUT VARIABLES:
A : polynomial coefficients of the best approximation
in order of increasing powers:
p*(x) = A(1) + A(2)*x + A(3)*x^2 + ...
REF : selected alternating set of points
HMAX : maximum deviation ( uniform norm of p* - f )
H : pointwise approximation errors
R : total number of iterations
EQUAL : success of failure of algorithm
EQUAL=1 : succesful
EQUAL=0 : convergence not acheived
EQUAL=-1: input error
EQUAL=-2: algorithm failure
Relies on function EXCH, provided below.
Example:
M = 5; N = 10000; K = 0; EPSH = 10^-12; MAXIT = 10;
X = linspace(-1,1,N); % uniformly spaced nodes on [-1,1]
k=1; Y = abs(X).^k; % the function Y to approximate
[A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
p = polyval(A,X); plot(X,Y,X,p) % p is the best approximation
Note: using an even value of M, e.g., M=2, in the example above, makes
the algorithm to fail with EQUAL=-2, because of collocation, which
appears because both the appriximating function and the polynomial are
even functions. The way aroung it is to approximate only the right half
of the function, setting K = 2 : even polynomial. For example:
N = 10000; K = 2; EPSH = 10^-12; MAXIT = 10; X = linspace(0,1,N);
for i = 1:2
k = 2*i-1; Y = abs(X).^k;
for j = 1:4
M = 2^j;
[~,~,HMAX] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
approxerror(i,j) = HMAX;
end
end
disp('Table 3.1 from Approximation theory and methods, M.J.D.POWELL, p. 27');
disp(' ');
disp(' n K=1 K=3');
disp(' '); format short g;
disp([(2.^(1:4))' approxerror']);
ALGORITHM:
Computation of the polynomial that best approximates the data (X,Y)
in the discrete uniform norm, i.e. the polynomial with the minimum
value of max{ | p(x_i) - y_i | , x_i in X } . That polynomial, also
known as minimax polynomial, is obtained by the exchange algorithm,
a finite iterative process requiring, at most,
n
( ) iterations ( usually p = M + 2. See also function EXCH ).
p
since this number can be very large , the routine may not converge
within MAXIT iterations . The other possibility of failure occurs
when there is insufficient floating point precision for the input
data chosen.
CREDITS: This routine was developed and modified as
computer assignments in Approximation Theory courses by
Prof. Andrew Knyazev, University of Colorado Denver, USA.
Team Fall 98 (Revision 1.0):
Chanchai Aniwathananon
Crhistopher Mehl
David A. Duran
Saulo P. Oliveira
Team Spring 11 (Revision 1.1): Manuchehr Aminian
The algorithm and the comments are based on a FORTRAN code written
by Joseph C. Simpson. The code is available on Netlib repository:
http://www.netlib.org/toms/501
See also: Communications of the ACM, V14, pp.355-356(1971)
NOTES:
1) A may contain the collocation polynomial
2) If MAXIT is exceeded, REF contains a new reference set
3) M, EPSH and REF can be altered during the execution
4) To keep consistency to the original code , EPSH can be
negative. However, the use of REF0 is *NOT* determined by
EPSH< 0, but only by its inclusion as an input parameter.
Some parts of the code can still take advantage of vectorization.
Revision 1.0 from 1998 is a direct human translation of
the FORTRAN code http://www.netlib.org/toms/501
Revision 1.1 is a clean-up and technical update.
Tested on MATLAB Version 7.11.0.584 (R2010b) and
GNU Octave Version 3.2.4
Package: optim