Adaptive numerical integration

Contents

Example for integral

We want to approximate the integral

with an approximation Q such that Q-I<=1e-2. We want to use the composite trapezoid rule with a subdivision which is adapted to the behavior of the function f: large intervals where |f''| is small, small intervals where |f''| is large.

We use 63 function evaluations. The actual error is less than 1e-3.

function adapt_test

f = @(x) (x.^3-x)./(1+x.^4);
ezplot(f,[0,6]); hold on

Nev = 2;                       % initial evaluations at a,b
Q = adaptint(f,0,6,1e-2)
Nev                            % number of function evaluations

title('subdivision for adaptive quadrature')
hold off;

I = log(1297)/4 - atan(36)/2;  % exact value for integral
error = abs(Q-I)               % Quadrature error

Q =
           1.0214243535841
Nev =
    63
error =
      0.000984902605729143

Recursive function adaptint(f,a,b,Tol)

We use the composite trapezoid rule. On each subinterval we estimate the error |QT-I| by |QT-QS| where QS is the Simpson rule value. We subdivide the intervals until the sum of estimated errors is <= Tol.

function Q = adaptint(f,a,b,Tol)
  fa = f(a); fb = f(b);
  QT = (b-a)/2*(fa+fb);        % Trapezoid rule

  c = (a+b)/2; fc = f(c);      % evaluate f in midpoint
  Nev = Nev + 1;               % count function evaluations
  plot([c c],[-1,1]*1e-2,'b'); % plot new point

  QS = (b-a)/6*(fa+4*fc+fb);   % Simpson rule
  % for small intervals QS is much closer to I than QT
  % hence we can approximate error QT-I by QT-QS

  if abs(QT-QS)<=Tol           % if estimated error is <= Tol
    Q = QT;                    %   accept trapezoid rule value
  else
    Q = adaptint(f,a,c,Tol/2) + adaptint(f,c,b,Tol/2);
                               %   use algorithm for [a,c] and [c,b]
  end
end

end

Published with MATLAB® R2016b