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