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

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