gausquad.m

matlab算法集 matlab算法集

function y = gausquad (a,b,n,f)
%-------------------------------------------------------------------- 
% Usage:       y = gausquad (a,b,n,f);
%
% Description: Use the n-point Gauss-Legendre formula to numerically
%              integrate the function f(x) over the interval [a,b].
%
% Inputs:      a = lower limit of integration
%              b = upper limit of integration
%              n = number of points (1 <= n <= 10)
%              f = string containing name of function to be integrated.   
%                  The form of f is:
%
%                       function y = f(x)
%
%                  When f is called with scalar xit must return the 
%                  value f(x).
%
% Outputs:     y = estimate of integral
%-------------------------------------------------------------------- 

% Initialize 

   y = 0;
   n = args (n,1,10,3,'gausquad');
   x = zeros (n,1);
   c = zeros (n,1);

% Compute parameters 

   switch (n) 

      case 1; x(1) = 0.0;
              c(1) = 2.0;

      case 2; x(1) = 0.5773503; x(2) = -x(1);
              c(1) = 1; c(2) = c(1);

      case 3; x(1) = 0.0;
              x(2) = 0.7745967; x(3) = -x(2);
              c(1) = 0.8888889;         
              c(2) = 0.5555556; c(3) = c(2);   

      case 4; x(1) = 0.3399810; x(2) = -x(1);
              x(3) = 0.8611363; x(4) = -x(3);
              c(1) = 0.6521452; c(2) = c(1);         
              c(3) = 0.3478548; c(4) = c(3);   

      case 5; x(1) = 0.0;
              x(2) = 0.5384693; x(3) = -x(2);
              x(4) = 0.9061798; x(5) = -x(4);
              c(1) = 0.5688880;
              c(2) = 0.4786287; c(3) = c(2);   
              c(4) = 0.2369269; c(5) = c(4);   

      case 6; x(1) = 0.2386192; x(2) = -x(1);
              x(3) = 0.6612094; x(4) = -x(3);
              x(5) = 0.9324695; x(6) = -x(5);
              c(1) = 0.4679139; c(2) = c(1);
              c(3) = 0.3607616; c(4) = c(3);  
              c(5) = 0.1713245; c(6) = c(5);  

      case 7; x(1) = 0.0;
              x(2) = 0.4058452; x(3) = -x(2);
              x(4) = 0.7415312; x(5) = -x(4);
              x(6) = 0.9491079; x(7) = -x(6);
              c(1) = 0.4179592;
              c(2) = 0.3818301; c(3) = c(2);   
              c(4) = 0.2797054; c(5) = c(4);   
              c(6) = 0.1294850; c(7) = c(6);   

      case 8; x(1) = 0.1834346; x(2) = -x(1);
              x(3) = 0.5255324; x(4) = -x(3);
              x(5) = 0.7966665; x(6) = -x(5);
              x(7) = 0.9620899; x(8) = -x(7);
              c(1) = 0.3626838; c(2) = c(1);
              c(3) = 0.3137066; c(4) = c(3);   
              c(5) = 0.2223810; c(6) = c(5);  
              c(7) = 0.1012285; c(8) = c(7);   

      case 9; x(1) = 0.0;
              x(2) = 0.3242534; x(3) = -x(2);
              x(4) = 0.6133714; x(5) = -x(4);
              x(6) = 0.8360311; x(7) = -x(6);
              x(8) = 0.9681602; x(9) = -x(8);
              c(1) = 0.3302394;
              c(2) = 0.3123471; c(3) = c(2);   
              c(4) = 0.2606107; c(5) = c(4);   
              c(6) = 0.1806482; c(7) = c(6);  
              c(8) = 0.0812744; c(9) = c(8);   

      case 10; x(1) = 0.1488743; x(2) = -x(1);
               x(3) = 0.4333954; x(4) = -x(3);
               x(5) = 0.6794096; x(6) = -x(5);
               x(7) = 0.8650634; x(8) = -x(7);
               x(9) = 0.9739065; x(10) = -x(9);
               c(1) = 0.2955242; c(2) = c(1);
               c(3) = 0.2692602; c(4) = c(3);   
               c(5) = 0.2190864; c(6) = c(5);   
               c(7) = 0.1494513; c(8) = c(7);   
               c(9) = 0.0666713; c(10) = c(9);  
      end
              
% Estimate integral 

   alpha = (b - a)/2;
   beta  = (b + a)/2;
   for i = 1 : n
      y = y + c(i)*feval(f,alpha*x(i)+beta);
   end
   y = alpha*y;