Home

Gnotide

Gnomesh

Nowcasting

Other Small Sotware

Personal Information

Adaptive Quadrature

Here is a very simple implementation in Matlab of an adaptive quadrature algorithm for approximating the integral of a function f(x) using a simpler simpson 1/3 rule.

function [ integral ] = adaptive( a,b,fx,e )
   s=simpson1_3(fx,a,b);
   s1=simpson1_3(fx,a,(a+b)/2);
   s2=simpson1_3(fx,(a+b)/2,b);
   if abs(s-s1-s2)<15*e
      disp(['accept a: ',num2str(a),' b: ',num2str(b),' h: ',num2str(b-a),' sum: ',num2str(s1+s2)]);
      integral=s1+s2;
   else
      integral=adaptive(a,(a+b)/2,fx,e/2)+adaptive((a+b)/2,b,fx,e/2);
   end
   

function [ integral ] = simpson1_3( fx, a, b )
   f=inline(fx);
   integral=(b-a)*(f(a)+4*f((a+b)/2)+f(b))/6;

Example:

If we are interested in the value of the defined integral of the function f(x)=exp(-x) in the limits [0,100] an analytical solution can be derived easily as F(x) = -exp(-100)-(-exp(-0)) = -(3.72007597602E-44) + 1 = 1 (with normal precision in any standar computer)

Using the above function for the simpson 1/3 rule "simpson1_3('exp(-(x))',0,100)", a rough stimation ban be obtained equal to 16.16667, what is really far away from real value.

Using the adaptive version, a much more acurate value of 1.002 is obtained even using a tolerance of 0.1.

adaptative(0,100,'exp(-(x))',.1)
accept a: 0 b: 3.125 h: 3.125 sum: 0.95791
accept a: 3.125 b: 6.25 h: 3.125 sum: 0.042087
accept a: 6.25 b: 12.5 h: 6.25 sum: 0.0019758
accept a: 12.5 b: 25 h: 12.5 sum: 4.5805e-006
accept a: 25 b: 50 h: 25 sum: 2.9157e-011
accept a: 50 b: 100 h: 50 sum: 8.0366e-022

In the figure, it is possible to see, how the points used for the adaptive method are denser where the slope is high, describing better the shape of the equation, and it has lower density where the function is almost flat, saving computational time.