driver_prob83.m

Integraton routines in matlab

%
% Driver for problem 8.3
%

for prob = 1:5

   if( prob == 1 )
      % error function
      f = inline('x./(1+x.^2)', 'x');
      a = 0; b = 3;
      integ_ex = log(10)/2;
      fprintf(1, ' f(x) = x/(1+x^2), [a,b] = [0,3] \n')
   end
   if( prob == 2 )
      % error function
      f = inline('1./(1-x)', 'x');
      a = 0; b = 0.95;
      integ_ex = log(20);
      fprintf(1, ' f(x) = 1/(1-x), [a,b] = [0,0.95] \n')
   end
   if( prob == 3 )
      % error function
      f = inline('1./sqrt(1-0.5*sin(x).^2)', 'x');
      a = 0; b = pi/2;
      integ_ex = 1;
      fprintf(1, ' f(x) = 1/sqrt(1-0.5*sin(x)^2), [a,b] = [0,pi/2] \n')
   end
   if( prob == 4 )
      % error function
      f = inline('1./sqrt(1-0.8*sin(x).^2)', 'x');
      a = 0; b = pi/2;
      integ_ex = 1;
      fprintf(1, ' f(x) = 1/sqrt(1-0.8*sin(x)^2), [a,b] = [0,pi/2] \n')
   end
   if( prob == 5 )
      % error function
      f = inline('1./sqrt(1-0.95*sin(x).^2)', 'x');
      a = 0; b = pi/2;
      integ_ex = 1;
      fprintf(1, ' f(x) = 1/sqrt(1-0.95*sin(x)^2), [a,b] = [0,pi/2] \n')
   end

   tol = 1.e-6;
   
   [integ, T, nfv] = trapcomp(f, a, b, tol);
   fprintf(1, ' Composite trapezoidal rule \n')
   fprintf(1, ' Approximate value of the integral = %12.6e \n', integ)
   fprintf(1, ' Absolute error                    = %12.6e \n', ...
              abs(integ-integ_ex) )
   fprintf(1, ' Relative error                    = %12.6e \n', ...
              abs(integ-integ_ex)/abs(integ_ex) )
   fprintf(1, ' Number of f evaluations required  = %6d \n\n', nfv )
   
   [integ, S, nfv] = simpcomp(f, a, b, tol);
   fprintf(1, ' Composite Simpson rule \n')
   fprintf(1, ' Approximate value of the integral = %12.6e \n', integ)
   fprintf(1, ' Absolute error                    = %12.6e \n', ...
              abs(integ-integ_ex) )
   fprintf(1, ' Relative error                    = %12.6e \n', ...
              abs(integ-integ_ex)/abs(integ_ex) )
   fprintf(1, ' Number of f evaluations required  = %6d \n\n', nfv )

   [integ, T, nfv] = romberg(f, a, b, tol);
   fprintf(1, ' Romberg Integration \n')
   fprintf(1, ' Approximate value of the integral = %12.6e \n', integ)
   fprintf(1, ' Absolute error                    = %12.6e \n', ...
              abs(integ-integ_ex) )
   fprintf(1, ' Relative error                    = %12.6e \n', ...
              abs(integ-integ_ex)/abs(integ_ex) )
   fprintf(1, ' Number of f evaluations required  = %6d \n\n', nfv )
end