function x=verenclinthull(A,b,d) % uses vertex % VERENCLINTHULL Verified enclosure of the interval hull of the solution set % of a system of interval linear equations. % % For a square interval matrix A and a matching interval vector b, % x=verenclinthull(A,b,d) % either computes a verified enclosure x of the interval hull of the % solution set of A*X=b, or yields no verified result. "d" is a search % depth parameter; it should be an integer between 0 and 6. % % There are two possible outputs: % % ~isnan(x) : x is a verified enclosure of the interval hull of A*X=b, % isnan(x) : no verified output. % % Generally said, the command % verenclinthull(A,b,d+1) % should give a better enclosure (although it is not theoretically proved) than % verenclinthull(A,b,d), % but it takes approximately twice as much computation time. The initial enclosure % verenclinthull(A,b,0) % is never worse than that computed by % verifylss(A,b). % % EXAMPLE (Barth and Nuding 1974). Computation done for d=0,1,2,3,4. % Observe how the enclosure is successively shrinking down to the optimal one: % % intval A = % [ 2.0000, 4.0000] [ -2.0000, 1.0000] % [ -1.0000, 2.0000] [ 2.0000, 4.0000] % intval b = % [ -2.0000, 2.0000] % [ -2.0000, 2.0000] % % tic, x0=verenclinthull(A,b,0), toc % intval x0 = % [ -14.0001, 14.0001] % [ -14.0001, 14.0001] % Elapsed time is 0.137479 seconds. % % tic, x1=verenclinthull(A,b,1), toc % intval x1 = % [ -14.0001, 14.0001] % [ -14.0001, 14.0001] % Elapsed time is 0.242569 seconds. % % tic, x2=verenclinthull(A,b,2), toc % intval x2 = % [ -12.0953, 12.0953] % [ -12.0953, 12.0953] % Elapsed time is 0.544185 seconds. % % tic, x3=verenclinthull(A,b,3), toc % intval x3 = % [ -4.0001, 4.0001] % [ -4.4445, 4.4445] % Elapsed time is 1.113682 seconds. % % tic, x4=verenclinthull(A,b,4), toc % intval x4 = % [ -4.0001, 4.0001] % [ -4.0001, 4.0001] % Elapsed time is 2.151929 seconds. % COMMENT. Optimal enclosure (interval hull) reached. This, however, is % not guaranteed in the general case. % % See also VERINTERVALHULL, VERIFYLSS. % Copyright 2005-2008 Jiri Rohn % % The branching mechanism (used in the subfunction VERTEX) is due to S. P. Shary. % The terminal enclosures are computed using the Hansen-Bliek-Rohn % enclosure as described in M. Fiedler, J. Nedoma, J. Ramik, J. Rohn % and K. Zimmermann, Linear Optimization Problems with Inexact Data, % Springer-Verlag, New York 2006, pp. 72-73. % % WARRANTY % % Because the program is licensed free of charge, there is % no warranty for the program, to the extent permitted by applicable % law. Except when otherwise stated in writing the copyright holder % and/or other parties provide the program "as is" without warranty % of any kind, either expressed or implied, including, but not % limited to, the implied warranties of merchantability and fitness % for a particular purpose. The entire risk as to the quality and % performance of the program is with you. Should the program prove % defective, you assume the cost of all necessary servicing, repair % or correction. % % History % % 2005-11-03 initial MATLAB version % 2006-02-15 INTLAB version % 2008-12-04 revisited; VERTEX created, embedded % 2008-12-11 final version of EA, p-coded % 2008-12-15 final version (html copy) % gr=getround; setround(0); b=b(:); [m,n]=size(A); x=repmat(infsup(NaN,NaN),n,1); if ~(m==n&&n==length(b)) % wrong data return end if ~(nargin==3&&d==round(d)&&0<=d&&d<=6) % d not assigned or not integer or out of the bounds d=0; % default end if ~isintval(A) A=infsup(A,A); % allows for real input end if ~isintval(b) b=infsup(b,b); end x=vertex(A,b,0,d); % till depth d setround(gr); % end of verenclinthull % %%%%%%%%%%%%%%%%%%%% % Subfunction vertex %%%%%%%%%%%%%%%%%%%% % function x=vertex(A,b,j,d) % recursive procedure; calls itself % j current depth (start with j=0) % d maximal depth allowed % interval hull guaranteed to be achieved for d=2^(n(n+1)) if j<d [m,n]=size(A); x=repmat(intval(NaN),n,1); x0=ea(A,b); Ab=[A b]; if any(any(Ab.inf<Ab.sup)) D=Ab.sup-Ab.inf; mx=max(max(D)); [i1,j1]=find(D==mx,1); if j1<=n A1=A; A1(i1,j1)=A.inf(i1,j1); x1=vertex(A1,b,j+1,d); A2=A; A2(i1,j1)=A.sup(i1,j1); x2=vertex(A2,b,j+1,d); else b1=b; b1(i1)=b.inf(i1); x1=vertex(A,b1,j+1,d); b2=b; b2(i1)=b.sup(i1); x2=vertex(A,b2,j+1,d); end if ~isnan(x1.inf(1))&&~isnan(x2.inf(1)) x=hull(x1,x2); if ~isnan(x0.inf(1)) x=intersect(x,x0); return end else if ~isnan(x0.inf(1)) x=x0; return end end else x=verifylss(A,b); end else % j>=d x=ea(A,b); end % end of vertex %