3.3 Krawczyk の定理に基づく解の検証アリゴリズム

function [ X , xs , k ] = verifynlss(f,xs,slopes,see,varargin)
%VERIFYNLSS   Verified solution of nonlinear system
%
%   [ X , xs , k ] = verifynlss(f,xs,slopes,see,P1,P2,...)
%
% f is name of function, to be called by f(xs), xs is approximation
%
% optional input    slopes  'g'  use gradient, proves uniqueness
%                           's'  use slopes, better, but w/o uniqueness
%                   see     see intermediate results
%                   P1,...  extra parameters for function evaluation
%                           f(x,P1,P2,...)
%
% optional output   xs      improved approximation (column vector)
%                   k       interval iteration steps
%
%
%Simple, one-dimensional nonlinear functions which can be written in one
%formula string, can be entered directly. The unknown must be 'x'. E.g.,
%    X = verifynlss('x*exp(x)-1',.6)
%evaluates an inclusion of the zero of x*exp(x)-1=0 near xs=.6.
%
%Nonlinear system solver based on the Krawcyzk operator, see
%  R. Krawczyk: Newton-Algorithmen zur Bestimmung von Nullstellen mit
%    Fehlerschranken, Computing 4, 187-201, 1969.
%  R.E. Moore: A Test for Existence of Solutions for Non-Linear Systems,
%    SIAM J. Numer. Anal. 4, 611-615, 1977.
%with modifications for enclosing the error with respect to an approximate
%solution, an iteration scheme, epsilon-inflation and an improved interval
%iteration as in
%  S.M. Rump: Solving Algebraic Systems with High Accuracy, in "A New
%    Approach to Scientific Computation", eds. U. Kulisch and W. Miranker,
%    Academic Press, 51-120, 1983.
%  S.M. Rump: Verification methods for dense and sparse systems of equations, 
%    in : J. Herzberger (ed.), Topics in Validated Computations - Studies in 
%    Computational Mathematics, Elsevier, Amsterdam, 63-136, 1994.
%
%Using gradient verifies existence and uniqueness of a zero of the nonlinear
%function within the inclusion interval. This also implies multiplicity 1 of
%the zero, and nonsingularity of the Jacobian at the zero.
%Using slopes implies existence but not uniqueness of a zero. This allows
%inclusion of zero clusters and multiple zeros. For details, see
%  S.M. Rump: Inclusion of Zeros of Nowhere Differentiable n-dimensional
%    Functions, Reliable Computing, 3:5-16 (1997).
%

% written  10/16/98     S.M. Rump
% modified 10/12/99     S.M. Rump  output NaN-vector in case of failure,
%                                  interval iteration stops if NaN occurred
% modified 06/26/02     S.M. Rump  output always of type intval, also for NaN
%

  xs = xs(:);
  if ( nargin<3 ) | isempty(slopes)
    slopes = 0;
  else
    if ~ischar(slopes)
      error('third parameter of verifynlss must be char')
    end
    slopes = isequal(lower(slopes(1)),'s');
  end

  if ( nargin<4 ) | isempty(see)
    see = 0;
  end

% Convert to inline function as needed
  strfun = fcnchk(f,length(varargin));

% floating point Newton iteration
  n = length(xs);
  xsold = xs;
  k = 0;
  while ( norm(xs-xsold)>1e-10*norm(xs) & k<10 ) | k<1
    k = k+1;                     % at most 10, at least 1 iteration performed
    xsold = xs;
    y = feval(strfun,gradientinit(xs),varargin{:});
    if see
      disp(['residual norm(f(xs_k)), floating point iteration ' sprintf('%d',k)])
      norm(y.x)
    end
    xs = xs - y.dx\y.x;
  end

% interval iteration
  R = inv(y.dx);
  Z = - R * feval(strfun,intval(xs),varargin{:});
  X = Z;
  E = 0.1*rad(X)*hull(-1,1) + midrad(0,realmin);
  ready = 0; k = 0; kmax = 10;
  while ( ~ready ) & ( k<kmax ) & ( ~any(isnan(X)) )
    k = k+1;
    if see
      disp(['interval iteration ' sprintf('%d',k)])
    end
    Y = hull( X + E , 0 );       % epsilon inflation
    Yold = Y;
    if slopes                    % use slopes
      x = slopeinit(xs,xs+Y);
      y = feval(strfun,x,varargin{:});
      C = eye(n) - R * y.s;      % automatic slopes
    else                         % use gradients
      x = gradientinit(xs+Y);
      y = feval(strfun,x,varargin{:});
      C = eye(n) - R * y.dx;     % automatic gradients
    end
    i=0;
    while ( ~ready ) & ( i<2 )   % improved interval iteration
      i = i+1;
      X = Z + C * Y;
      ready = all(all(in0(X,Y)));
      Y = intersect(X,Yold);
    end
  end
  if ready
    X = xs+Y;                    % verified inclusion
  else
    X = intval(repmat(NaN,n,1)); % inclusion failed
  end

option+a å

option+u,a ä