Code covered by the BSD License  

Highlights from
MATLAB en zijn Symbolic Math Toolbox

image thumbnail

MATLAB en zijn Symbolic Math Toolbox

by

Tjibbele Miedema

 

Leergang voor geïnteresseerden in de Computeralgebra - Systemen MATLAB en zijn Symbolic Math Toolbox

ZeroBisect(f, a, b, aTol, rTol, maxIter, trace)
function [xZero, fZero, nIter] = ZeroBisect(f, a, b, aTol, rTol, maxIter, trace)

%
%                 ************
% Funnctie-file : ZeroBisect.m
%                 ************

% Bepaling van het nulpunt van de functie f met behulp van de bisectie methode
%
%  Input parameters : 
% -------------------
%
% --- f          : string die de functie definieert
% --- [a, b]     : het interval waar het nulpunt zich in bevindt
% --- aTol, rTol : respectievelijk absolute- en relatieve tolerantie
% --- maxIter    : maximaal aantal iteraties
% --- trace      : trace-functie : 1 'aan', 0 'af'
%
% Output parameters :
% -------------------
%
% --- xZero      : gevonden nulpunt
% --- fZero      : f(xZero)
% --- nIter      : benodigd aantal iteraties
%

fa  = feval(f, a);
sfa = sign(fa);
fb  = feval(f, b);
sfb = sign(fb);

if sfa == sfb
   error('Nulpuntsbepaling onmogelijk : geen nulpunts-interval')
end

if     sfa == 0
   xZero = a;
   fZero = 0.0;
   nIter = 0;
elseif sfb == 0
   xZero = b;
   fZero = 0.0;
   nIter = 0;
else
   if trace
      disp(' ')
      disp(sprintf ...
         ('**** traceer - optie ****                     a = %15.12e', a))
      disp(sprintf ...
         ('                                              b = %15.12e', b))
      disp(' ')
      disp(sprintf ...
         (' k             a                    b                   f(M)'))
   end
   err = abs(b - a);
   for k = 1 : maxIter
      err = err / 2;
      M   = (a + b) / 2;
      fM  = feval(f, M);
      sfM = sign(fM);
      if trace
         if rem(k, 10) == 1
            disp(' ')
         end % if rem  
         disp(sprintf('%3d  %+15.12e  %+15.12e  %+15.12e', k, a, b, fM))
      end % if trace
      if (sfM == 0) | (err < abs(M) * rTol + aTol)
         xZero = M;
         fZero = fM;
         nIter = k;
         break; % binnen tolerantie
      else
         if sfM == sfa
            a  = M;
            fa = fM;
         else
            b  = M;
            fb = fM;
         end
      end       % buiten tolerantie
   end % for k
   if k == maxIter
      warning('Nulpuntsbepaling niet binnen de vereiste nauwkeurigheid')
   else
      disp(' ')
      disp('Bisectie-methode convergeerde binnen de vereiste nauwkeurigheid')
   end % if k
end %else

%
%                ************
% Functie-file : ZeroBisect.m
%                ************
%

Contact us