Two Incomplete Matlab Codes for Numerical Integration

 Related to a Computer Project
 


Note:  You need to modify the codes to get them work. 




% integral1.m: comparison of numerical integration rules 

   a = 0; b = pi/2; M = 10;

   for n = 1:M
     h = (b-a)/n; 

     rl = 0; rm = 0;
     for i = 1:n
       rl = rl + sin(a+(i-1)*h);         % left endpont rectangle rule
       rm = rm + sin( ? );               % midpoint rectangle rule
     end 
     rl = h*rl; erl(n) = log(abs(rl-1));
     rm = h*rm; erm(n) = log(abs(rm-1)); 

     t = sin(a) + sin(b);                % trapezoid rule
     for i = 1:n-1
       t = t + ?;                   
     end 
     t = 0.5*h*t; et(n) = log(abs(t-1));

     s = sin(a) + sin(b);                % Simpson's rule
     for i = 1:n-1
       s = s + 2*sin(a+i*h);
     end  
     for i = 1:n
       s = s + ?*sin( ? );
     end
     s = s*h/6; es(n) = log(abs(s-1));

     g2 = 0;                             % Two-point Gaussian quadrature
     for i = 1:n
       g2 = g2 + sin(a+(i-0.5)*h-0.5*h/sqrt(3)) ...
            + sin(a+(i-0.5)*h+0.5*h/sqrt(3));
     end 
     g2 = 0.5*h*g2; eg2(n) = log(abs(g2-1));

     disp(sprintf('%3.0f %12.8f %12.8f %12.8f %12.8f %12.8f',n,rl,rm,t,s,g2))
   end

   x = [];
   for n = 1:M
     x(n) = log(n);
   end 
   plot(x, erl, 'o', x, erm, '--', x, et, '-', x, es, '-.', x, eg2, ':')
   print error1.ps

% integral2.m : comparison of the trapezoid rule, the extrapolated 
%               trapezoid rule, and the Simpson rule. 

   a = 0; b = pi/2; M = 8;

   for n = 1:M+1
     k = 2^(n-1);
     h = (b-a)/k;

     t = sin(a) + sin(b);                     % trapezoid rule
     for i = 1:k-1
       t = t + 2*sin(a+i*h);
     end 
     t = 0.5*h*t; tra(n) = t;

     s = sin(a) + sin(b);                     % Simpson's rule
     for i = 1:k-1
       s = s + 2*sin(a+i*h);
     end
     for i = 1:k
       s = s + 4*sin(a+(i-0.5)*h);
     end
     s = s*h/6; sim(n) = s;
   end

   for n = 1:M                                % Richardson extrapolation
     etr(n) = ? ;             
   end 

   for n = 1:M
     k = 2^(n-1);
     disp(sprintf('%3.0f %20.16f %20.16f %20.16f',k,tra(n),etr(n),sim(n)))
   end
   disp(sprintf('%3.0f %20.16f %41.16f',256,tra(9),sim(9)))
   disp(sprintf(''))

   for n = 1:M
     k = 2^(n-1);
     eetr(n) = log(abs(etr(n)-1));
     etra(n) = log(abs(tra(n)-1));
     esim(n) = log(abs(sim(n)-1));
     disp(sprintf('%3.0f %20.16f %20.16f %20.16f',k,etra(n),eetr(n),esim(n)))
   end

   x = [];
   for n = 1:M
     x(n) = log(2^(n-1));
   end 
   plot(x, etra, ':', x, eetr, '-', x, esim, 'o')
   print error2.ps