Code covered by the BSD License  

Highlights from
MATLAB en zijn Symbolic Math Toolbox

image thumbnail

MATLAB en zijn Symbolic Math Toolbox

by

 

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

IntSimpRichcorr(f, a, b, aTol, rTol, ...
function [Integrl, nInterv ] = IntSimpRichcorr(f, a, b, aTol, rTol, ...
   maxInterv, trace)

%
%                ******************
% Functie-file : IntSimpsRichcorr.m
%                ******************
%
% Bepaling van de integraal van de functie f over het interval [a, b]
% Toegepaste methode : Regel van Simpson, samen met de Richardson-correctie
%
%  Input parameters :
% -------------------
%
% --- f          : string die de te integreren functie definieert
% --- [a, b]     : respectievelijk absolute- en relatieve tolerantie
% --- aTol, rTol : maximaal aantal iteraties
% --- maxInterv  : maximaal aantal intervallen (gehele macht van 2)
% --- trace      : trace-functie : 1 'aan', 0 'af'
%
% Output parameters :
% -------------------

% --- Integrl    : waarde van de bepaalde integraal
% --- nInterv    : benodigd aantal intervallen
%

if trace
   disp(' ')
   disp(sprintf( ...
      '**** traceer - optie ****                      a = %15.12e', a))
   disp(sprintf( ...
      '                                               b = %15.12e', b))
   disp(' ')
   disp(sprintf( ...
      '  n             I2h                   Ih                 rc/Int'))
   disp(' ')
end

n   = 1;
h   = b - a;
fa  = feval(f, a);
fb  = feval(f, b);

s2h = 0;
sh  = feval(f, a + h / 2);
Ih  = (fa + fb + 4 * sh) * (h / 2) / 3; 

while (n * 2) < maxInterv 
   I2h = Ih;
   s2h = sh + s2h;
   sh  = 0;
   
   n   = n * 2;
   h   = h / 2;
   
   x   = a + h / 2;
   
   for i = 1 : n
      sh = sh + feval(f, x);
      x  = x + h;
   end % for i
   
   Ih  = (fa + fb + 4 * sh + 2 * s2h) * (h / 2) / 3;
   Int = (16 * Ih - I2h) / 15;
   
   if trace == 1
      disp(sprintf( ...
         '%4d  %+15.12e  %+15.12e  %+15.12e', n * 2, I2h, Ih, Int - Ih))
      disp(sprintf( ...
         '                                                  %+15.12e', Int))
   end % if trace
   
   if (abs(Int - Ih) < abs(Int) * rTol + aTol) 
      break;
   end % if (abs(Int - Ih)
end % while

Integrl = Int;
nInterv = n * 2;

if n * 2 >= maxInterv
   warning('Geen voldoende nauwkeurigheid : maximaal aantal intervallen')
end

return

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

Contact us