1 ## Copyright (c) 1998-2011 Andrew V. Knyazev <andrew.knyazev@ucdenver.edu>
2 ## All rights reserved.
4 ## Redistribution and use in source and binary forms, with or without",
5 ## modification, are permitted provided that the following conditions are met:
7 ## 1 Redistributions of source code must retain the above copyright notice,
8 ## this list of conditions and the following disclaimer.
9 ## 2 Redistributions in binary form must reproduce the above copyright
10 ## notice, this list of conditions and the following disclaimer in the
11 ## documentation and/or other materials provided with the distribution.
12 ## 3 Neither the name of the author nor the names of its contributors may be
13 ## used to endorse or promote products derived from this software without
14 ## specific prior written permission.
16 ## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
17 ## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
18 ## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
19 ## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
20 ## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
21 ## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
22 ## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
23 ## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
24 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 % function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
29 % Best polynomial approximation in discrete uniform norm
33 % M : degree of the fitting polynomial
34 % N : number of data points
35 % X(N) : x-coordinates of data points
36 % Y(N) : y-coordinates of data points
37 % K : character of the polynomial:
38 % K = 0 : mixed parity polynomial
39 % K = 1 : odd polynomial ( X(1) must be > 0 )
40 % K = 2 : even polynomial ( X(1) must be >= 0 )
41 % EPSH : tolerance for leveling. A useful value for 24-bit
42 % mantissa is EPSH = 2.0E-7
43 % MAXIT : upper limit for number of exchange steps
44 % REF0(M2): initial alternating set ( N-vector ). This is an
45 % OPTIONAL argument. The length M2 is given by:
46 % M2 = M + 2 , if K = 0
47 % M2 = integer part of (M+3)/2 , if K = 1
48 % M2 = 2 + M/2 (M must be even) , if K = 2
52 % A : polynomial coefficients of the best approximation
53 % in order of increasing powers:
54 % p*(x) = A(1) + A(2)*x + A(3)*x^2 + ...
55 % REF : selected alternating set of points
56 % HMAX : maximum deviation ( uniform norm of p* - f )
57 % H : pointwise approximation errors
58 % R : total number of iterations
59 % EQUAL : success of failure of algorithm
61 % EQUAL=0 : convergence not acheived
62 % EQUAL=-1: input error
63 % EQUAL=-2: algorithm failure
65 % Relies on function EXCH, provided below.
68 % M = 5; N = 10000; K = 0; EPSH = 10^-12; MAXIT = 10;
69 % X = linspace(-1,1,N); % uniformly spaced nodes on [-1,1]
70 % k=1; Y = abs(X).^k; % the function Y to approximate
71 % [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
72 % p = polyval(A,X); plot(X,Y,X,p) % p is the best approximation
74 % Note: using an even value of M, e.g., M=2, in the example above, makes
75 % the algorithm to fail with EQUAL=-2, because of collocation, which
76 % appears because both the appriximating function and the polynomial are
77 % even functions. The way aroung it is to approximate only the right half
78 % of the function, setting K = 2 : even polynomial. For example:
80 % N = 10000; K = 2; EPSH = 10^-12; MAXIT = 10; X = linspace(0,1,N);
82 % k = 2*i-1; Y = abs(X).^k;
85 % [~,~,HMAX] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT);
86 % approxerror(i,j) = HMAX;
89 % disp('Table 3.1 from Approximation theory and methods, M.J.D.POWELL, p. 27');
92 % disp(' '); format short g;
93 % disp([(2.^(1:4))' approxerror']);
97 % Computation of the polynomial that best approximates the data (X,Y)
98 % in the discrete uniform norm, i.e. the polynomial with the minimum
99 % value of max{ | p(x_i) - y_i | , x_i in X } . That polynomial, also
100 % known as minimax polynomial, is obtained by the exchange algorithm,
101 % a finite iterative process requiring, at most,
103 % ( ) iterations ( usually p = M + 2. See also function EXCH ).
105 % since this number can be very large , the routine may not converge
106 % within MAXIT iterations . The other possibility of failure occurs
107 % when there is insufficient floating point precision for the input
110 % CREDITS: This routine was developed and modified as
111 % computer assignments in Approximation Theory courses by
112 % Prof. Andrew Knyazev, University of Colorado Denver, USA.
114 % Team Fall 98 (Revision 1.0):
115 % Chanchai Aniwathananon
120 % Team Spring 11 (Revision 1.1): Manuchehr Aminian
122 % The algorithm and the comments are based on a FORTRAN code written
123 % by Joseph C. Simpson. The code is available on Netlib repository:
124 % http://www.netlib.org/toms/501
125 % See also: Communications of the ACM, V14, pp.355-356(1971)
129 % 1) A may contain the collocation polynomial
130 % 2) If MAXIT is exceeded, REF contains a new reference set
131 % 3) M, EPSH and REF can be altered during the execution
132 % 4) To keep consistency to the original code , EPSH can be
133 % negative. However, the use of REF0 is *NOT* determined by
134 % EPSH< 0, but only by its inclusion as an input parameter.
136 % Some parts of the code can still take advantage of vectorization.
138 % Revision 1.0 from 1998 is a direct human translation of
139 % the FORTRAN code http://www.netlib.org/toms/501
140 % Revision 1.1 is a clean-up and technical update.
141 % Tested on MATLAB Version 7.11.0.584 (R2010b) and
142 % GNU Octave Version 3.2.4
144 % $Revision: 1.1 $ $Date: 2011/08/3 $
146 % ************************************ beginning of POLYFITINF
147 function [A,REF,HMAX,H,R,EQUAL] = polyfitinf(M,N,K,X,Y,EPSH,MAXIT,REF0)
149 % Preassign output variables A,REF,HMAX,H,R,EQUAL in case of error return
150 A = []; REF = []; HMAX = []; H = []; R = 0; EQUAL = -2;
151 %%%% end preassignment
153 % Setting M with respect to K
169 % If the user has input odd M, but wants an even polynomial,
170 % subtract 1 from M to prevent errors later. The outputs should be
171 % mathematically equivalent.
179 warning('polyfitinf:MixedParity','Using mixed parity polynomial...');
189 % Check input data consistency
191 if ( (length(X) ~= N) || (length(Y) ~= N) )
192 error('Input Error: check data lengths');
196 error('Input Error: insufficient data points');
200 error('Input Error: insufficient degree');
203 if ( (K == 2) && (X(1) < 0) ) || ( (K == 1) && (X(1) <= 0) )
204 error('Input Error: X(1) inconsistent with parity');
208 error('Input Error: Abscissae out of order');
223 % Read initial reference set into Z, if available.
228 Z(2:(P+1))= REF0(1:P);
230 % Check if REF is monotonically increasing
231 if ( any(diff(REF0) < 0) || any(REF0 > J) )
232 error('Input Error : Bad initial reference set');
237 % Loads Z with the points closest to the Chebychev abscissas
242 % Setting parity-dependent parameters
251 ITEMP = 2*(M+1) + Q1;
255 % Calculate the J-th Chebyshev abcissa and load Z(J+1)
256 % with the appropriate index from the data abcissas
260 X1 = XA + XE*( cos(Q*(P-J)) );
267 while ( (II <= HIGH) && (FLAG == 1) )
271 % If the Chebyschev abscissa is bracketed by
272 % two input abcissas, get out of the while loop
274 if (X(I)+X(ITEMP) <= X1)
293 % If the lower Chebyshev abcissas are less than X(1),
294 % load the lower elements of Z with the lowest points
296 IND = find(Z(2:end) >= (1:(length(Z)-1)));
298 try TEMP = IND(1); % If IND is empty, do nothing.
299 catch exception % The catch will be that IND is an empty array.
301 if strcmpi(exception.identifier,'MATLAB:badsubscript')
302 % This will be the exception. Do nothing.
307 Z(2:TEMP) = (1:(TEMP-1))';
313 % M1 entry. Initialize variables to prepare for exchange iteration
324 % Load H with the ordinates and XX(I) with the abscissas if the
325 % polynomial is mixed . If it is even or odd , load XX with the
326 % squares of the abscissas.
341 % Iteration entry. R is the iteration index
348 while ( (R < MAXIT) && (FLAG == 1) )
350 R = R + 1; % LABEL 350
353 % Computation of div. differences schemes
357 % If the polynomial is mixed or even:
363 % C(I) = (H(J) + S*T)/Q;
371 C(I) = (H(J) + S*T)./X(J);
377 % If the polynomial is odd:
382 % C(I) = H(ITEMP) + S*T;
389 C(I) = H( Z(ITEMP) ) + S.*T;
405 C(J) = (C(J)-C(ITEMP))/QD;
406 D(J) = (D(J)-D(ITEMP))/QD;
413 % Computation of polynomial coefficients
419 DAA(ITEMP) = C(ITEMP) + DT*D(ITEMP);
425 DAA(LOW:M) = DAA(LOW:M) - QD*DAA(((LOW:M)+1));
429 AA(1:HIGH) = AA(1:HIGH) + DAA(1:HIGH);
431 % Evaluation of the polynomial to get the approximation errors
444 % If the polynomial is odd, multiply SD by X(I)
452 % Load MAXX with the largest magnitude
453 % of the approximation array
458 % Test for alternating signs
464 % This represents a case where the polynomial
465 % exactly predicts a data point
467 warning('polyfitinf:Collocation','Collocation has occured.');
474 warning('polyfitinf:AnotherTry','1 more attempt with middle points');
475 LOW = (N+1)/2 - (P+1)/2 + 1;
477 Z(LOW:HIGH) = ( (LOW:HIGH) -1);
480 disp('Normal Exit.');
494 while ( (I <= P) && (FLAG2 == 1) )
499 warning('polyfitinf:Collocation','Collocation has occured.');
506 warning('polyfitinf:AnotherTry','1 more attempt with middle points');
507 LOW = (N+1)/2 - (P+1)/2 + 1;
509 Z(LOW:HIGH) = ( (LOW:HIGH) -1);
511 disp('Normal Exit.');
524 % Error entry: bad accuracy for calculation
538 % Search for another reference
542 [H,Z,EQUAL] = exch(N, P, EPSH, H, Z);
554 end % end of if over H(ITEMP)
556 end; % end of iteration loop
558 % M2 entry; load output variables and return
562 % Load the coefficients into A array
564 A(Q1 + Q2*(((1:HIGH)-1)) + 1) = AA(1:HIGH);
566 % Load REF with the final reference points
568 REF(1:P) = Z((1:P) + 1);
574 warning('polyfitinf:Collocation','polyfitinf terminates');
578 warning('polyfitinf:NoAlternatingSigns','Alternating signs not observed');
582 warning('polyfitinf:MaxIterationsReached','MAXIT was reached, current ref. set saved in REF.');
585 % Reverse the order of A to make it compatible with MATLAB'S polyval() function.
590 % ****************************************** end of POLYFITINF
592 function [H,Z,EQUAL] = exch(N, P, EPSH, H, Z)
593 % function [H,Z,EQUAL] = exch(N, P, EPSH, H, Z)
595 % EXCH: exchange algorithm
598 % N : number of data points
599 % P : number of reference points
600 % EPSH : tolerance for leveling.
601 % Z : old reference indices
604 % H : pointwise approximation errors
605 % Z : new reference indices
606 % EQUAL : EQUAL=1 : normal exchange
607 % EQUAL=0 : old and new references are equal
609 % CREDITS: This routine was developed and modified as
610 % computer assignments in Approximation Theory courses by
611 % Prof. Andrew Knyazev, University of Colorado Denver, USA.
613 % Team Fall 98 (Revision 1.0):
614 % Chanchai Aniwathananon
619 % Team Spring 11 (Revision 1.1): Manuchehr Aminian
621 % The algorithm and the comments are based on a FORTRAN code written
622 % by Joseph C. Simpson. The code is available on Netlib repository:
623 % http://www.netlib.org/toms/501
624 % See also: Communications of the ACM, V14, pp.355-356(1971)
626 % Revision 1.0 from 1998 is a direct human translation of
627 % the FORTRAN code http://www.netlib.org/toms/501
628 % Revision 1.1 is a clean-up and technical update.
629 % Tested on MATLAB Version 7.11.0.584 (R2010b) and
630 % GNU Octave Version 3.2.4
633 % Copyright 1998-2011 Andrew V. Knyazev
634 % $Revision: 1.1 $ $Date: 2011/05/17 $
636 % ************************************ beginning of exch
642 % SIG is arbitrarily chosen equal to the sign of the input
643 % point. This will be adjusted later if necessary.
651 % The next loop prescans Z to insure it is a proper choice, i.e
652 % resets Z if necessary so that maximum error points are chosen,
653 % given the sign convention mentioned above. In order to work
654 % properly, this section requires Z(1) = 0 and Z(P+2) = N + 1 .
663 % Scan the open point interval containing only the 1th initial
664 % reference point. In the interval pick the point with largest
665 % magnitude and correct sign. Most of the sorting occurs in
666 % this section. SIG contains the sign assumed for H(I).
669 if (SIG*(H(J)-MAXX) > 0)
678 % If the MAX error is significantly greater than the
679 % input point, switch to this point.
681 if (abs( MAXX - H(ITEMP) )/MAXL > EPSH)
695 % Find the error with largest abs value and proper sign
696 % from among the points above the last reference point.
697 % This section is necessary because the set of points
698 % chosen may begin with the wrong sign alternation.
701 if (SIG*(MAXR-H(J)) > 0)
708 % Find the error with largest abs value and proper sign
709 % from among the points below the 1st reference point.
710 % This section is necessary by the same reason as above.
725 if (SIG*(MAXL-H(J)) > 0)
733 % MAXL and MAXR contain the magnitude of the significant
734 % errors outside the reference point set. If either is
735 % zero, the reference point set extends to the end point
736 % on that side of the interval.
745 % L = 0 implies that the previous prescan did not change
746 % any points. If L = 0 and MAXL, MAXR are not significant
747 % if compared with upper and lower reference point errors,
748 % respectively, use the EQUAL exit.
752 if ( (MAXL == 0) || (EPSH >= (MAXL-HZP)/MAXL) )
753 if ( (MAXR == 0) || (EPSH >= (MAXR-HZ1)/MAXR) )
760 if ( (MAXL == 0) && (MAXR == 0) )
764 if ( (MAXL > MAXR) && (MAXL <= HZP) && (MAXR < HZ1) )
768 if ( (MAXL <= MAXR) && (MAXR <= HZ1) && (MAXL < HZP) )
772 % If a point outside the present reference set must be
773 % included, (i.e. the sign of the 1st point previously
774 % assumed is wrong) shift to the side of largest
775 % relative error first.
781 if ( (MAXL > MAXR) && (MAXL > HZP) )
785 if ( (MAXL <= MAXR) && (MAXR <= HZ1) )
791 % SHR entry. This section inserts a point from
792 % above the prescan point set
796 % shift point set down, dropping the lowest point
802 % add the next high point
805 % if MAXL > 0 replace reference points from the left,
806 % stopping the 1st time the candidate for replacement
807 % is greater than in magnitude than the prospective
808 % replacee. Alternation of signs is preserved because
809 % non-replacement immediately terminates the process.
814 while ( (I <= P) && (FLAG3 == 1) )
816 if ( abs(H(INDL)) >= abs(H(ITEMP)) )
830 % SHL entry. This section inserts a point from below the
838 % store lowest point in Z(2)
841 % if MAXR > 0 replace reference points from the right,
842 % stopping the 1st time the candidate for replacement
843 % is greater than in magnitude than the prospective
849 while ( (II <= P) && (FLAG3 == 1) )
853 if ( abs(H(INDR)) >= abs(H(HIGH)) )
870 % ****************************************** end of exch