Adaptive Quadrature Algorithm using MATLAB (m file)



MATLAB Program:

% Adaptive quadrature algoritm  
 % Find the integral of y=sin(x) from 0 to pi.

 f = @(x)  sin(x);
 aa = input('Enter lower limit, a:  ');
 bb = input('Enter upper limit, b:  ');
 n = input('Enter maximum no. of levels, n:  ');
 eps = input('Enter tolerance, tol:  ');


 cnt = 0;
 app = 0;
 i = 1;
 tol = zeros(1,n);
 a = zeros(1,n);
 h = zeros(1,n);
 fa = zeros(1,n);
 fc = zeros(1,n);
 fb = zeros(1,n);
 s = zeros(1,n);
 l = zeros(1,n);
 fd = zeros(1,n);
 fe = zeros(1,n);
 v = zeros(1,7);
 tol(i) = 10*eps;
 a(i) = aa;
 h(i) = 0.5*(bb-aa);
 fa(i) = f(aa);
 cnt = cnt+1;
 fc(i) = f((aa+h(i)));
 cnt = cnt+1;
 fb(i) = f(bb);
 cnt = cnt+1;


 s(i) = h(i)*(fa(i)+4*fc(i)+fb(i))/3;
 l(i) = 1;

 while i > 0 
   fd = f((a(i)+0.5*h(i)));
   cnt = cnt+1;
   fe = f((a(i)+1.5*h(i)));
   cnt = cnt+1;
   s1 = h(i)*(fa(i)+4*fd+fc(i))/6;
   s2 = h(i)*(fc(i)+4*fe+fb(i))/6;
   v(1) = a(i);
   v(2) = fa(i);
   v(3) = fc(i);
   v(4) = fb(i);
   v(5) = h(i);
   v(6) = tol(i);
   v(7) = s(i);
   lev = l(i);
   i = i-1;
   if abs(s1+s2-v(7)) < v(6)
     app = app+(s1+s2);
   else
     if lev >= n
       fprintf('Procedure fails');
     else
      i = i+1;
      a(i) = v(1)+v(5);
      fa(i) = v(3);
      fc(i) = fe;
      fb(i) = v(4);
      h(i) = 0.5*v(5);
      tol(i) = 0.5*v(6);
      s(i) = s2;
      l(i) = lev+1;
      i = i+1;
      a(i) = v(1);
      fa(i) = v(2);
      fc(i) = fd;
      fb(i) = v(3);
      h(i) = h(i-1);
      tol(i) = tol(i-1);
      s(i) = s1;
      l(i) = l(i-1);
     end
   end
 end
 fprintf('The approximate integral is: %11.8f \n', app);


OUTPUT: 

>> Adaptive_quadrature_c
Enter lower limit, a:  0
Enter upper limit, b:  pi
Enter maximum no. of levels, n:  6
Enter tolerance, tol:  0.001
The approximate integral is:  2.00026917 
>> 

1 comment:

  1. a=0;
    b=0.8;
    n=2;
    y= @(x) 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;

    choice = menu ( 'choose a methode','trapizoidal','simpson1/3','simpson3/8','Gaussian Quadratute Algorithm');
    switch choice
    case 1
    % trapziodal ;
    h=(b-a)/n;
    s=0.5*(y(a)+y(b));
    for i=1 :n-1
    s=s+y(a+i*h);
    end
    I=h*s
    case 4

    % Gaussian Quadratute Algorithm;
    % Find the integral of f=0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5 from
    % 0 to 0.8.

    a = input('Enter lower limit, a: ');
    b = input('Enter upper limit, b: ');
    n = input('Enter the order, n: ');

    if (n < 2 | n > 8)
    fprintf('This code works only for order from 2 to 8\n');
    else
    a0 = 0.5*(b+a);
    a1 = 0.5*(b-a);
    if (n == 2)
    c(1) = 1.0;
    c(2) = c(1);
    x_table(1) = -0.577350269;
    x_table(2) = -x_table(1);
    elseif (n == 3)
    c(1) = 0.555555556;
    c(2) = 0.888888889;
    c(3) = c(1);
    x_table(1) = -0.774596669;
    x_table(2) = 0.0;
    x_table(3) = -x_table(1);
    elseif (n == 4)
    c(1) = 0.347854845;
    c(2) = 0.652145155;
    c(3) = c(2);
    c(4) = c(1);
    x_table(1) = -0.861136312;
    x_table(2) = -0.339981044;
    x_table(3) = -x_table(2);
    x_table(4) = -x_table(1);
    elseif (n == 5)
    c(1) = 0.236926885;
    c(2) = 0.478628670;
    c(3) = 0.568888889;
    c(4) = c(2);
    c(5) = c(1);
    x_table(1) = -0.906179846;
    x_table(2) = -0.538469310;
    x_table(3) = 0.0;
    x_table(4) = -x_table(2);
    x_table(5) = -x_table(1);
    elseif (n == 6)
    c(1) = 0.171324492;
    c(2) = 0.360761573;
    c(3) = 0.467913935;
    c(4) = c(3);
    c(5) = c(2);
    c(6) = c(1);
    x_table(1) = -0.932469514;
    x_table(2) = -0.661209386;
    x_table(3) = -0.238619186;
    x_table(4) = -x_table(3);
    x_table(5) = -x_table(2);
    x_table(6) = -x_table(1);
    elseif (n == 7)
    c(1) = 0.129484966;
    c(2) = 0.279705391;
    c(3) = 0.381830050;
    c(4) = 0.417959184;
    c(5) = c(3);
    c(6) = c(2);
    c(7) = c(1);
    x_table(1) = -0.949107912;
    x_table(2) = -0.741531186;
    x_table(3) = -0.405845151;
    x_table(4) = 0.0;
    x_table(5) = -x_table(3);
    x_table(6) = -x_table(2);
    x_table(7) = -x_table(1);
    elseif (n == 8)
    c(1) = 0.101228536;
    c(2) = 0.222381034;
    c(3) = 0.313706645;
    c(4) = 0.362683783;
    c(5) = c(4);
    c(6) = c(3);
    c(7) = c(2);
    c(8) = c(1);
    x_table(1) = -0.960289856;
    x_table(2) = -0.796666477;
    x_table(3) = -0.525532409;
    x_table(4) = -0.183434642;
    x_table(5) = -x_table(4);
    x_table(6) = -x_table(3);
    x_table(7) = -x_table(2);
    x_table(8) = -x_table(1);
    end

    integral = 0.0;
    for i = 1:n
    x(i) = a0 + a1*x_table(i);
    f(i) = funkeval(x(i));
    integral = integral + c(i)*f(i);
    end
    integral = integral*a1;
    fprintf('The integral is: %10.8f\n\n',integral);
    end
    FUNCTION f(i) = funkeval(x(i));
    f=0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5
    end

    >> UntitledRAQ
    Enter lower limit, a: 0
    Enter upper limit, b: 0.8
    Enter maximum no. of levels, n: 4
    Enter tolerance, tol: 0.001
    The approximate integral is: 0.30329600
    >> UntitledRAQ
    Enter lower limit, a: 0
    Enter upper limit, b: PI
    Error using input
    Undefined function or variable 'PI'.

    Error in UntitledRAQ (line 8)
    bb = input('Enter upper limit, b: ');

    Enter upper limit, b: CLC
    Error using input
    Undefined function or variable 'CLC'.

    Error in UntitledRAQ (line 8)
    bb = input('Enter upper limit, b: ');

    Enter upper limit, b: PI
    Error using input
    Undefined function or variable 'PI'.

    Error in UntitledRAQ (line 8)
    bb = input('Enter upper limit, b: ');

    Enter upper limit, b: pi
    Enter maximum no. of levels, n: 6
    Enter tolerance, tol: 0.001
    The approximate integral is: 2.00026917
    >>


    ReplyDelete