[CSC 335] Gaussian Quadrature

[Andrew J. Pounds]

  So from class today we never really got to play with the 
Gauss-Legendre quadrature and also your book doesn't provide any code -- 
so I have attached two small pieces of code to integrate f(x) = e^x cos(x).

See what you think...

-- 
Andrew J. Pounds, Ph.D.  (pounds at theochem.mercer.edu)
Associate Professor of Chemistry and Computer Science
Mercer University,  Macon, GA 31207   (478) 301-5627

-------------- next part --------------
      program gleg

****  Gauss-Laguerre Qaudrature

****  Program will integrate functon exp(x)*cos(x)

      double precision X(2), W(2), sum, fofx, A, B
      external fofx

          X( 1) =  0.5773502692D0
          X( 2) =  -0.5773502692D0
          W( 1) =  1.0D0
          W( 2) =  1.0D0

      A = -1.0
      B =  1.0 
      sum = 0.0
      do 10 i = 1, 2
         sum = sum + W(i) * fofx(((B-A)*X(i)+(A+B))/(B-A))*(B-A)/2.0
10    continue

      print *, 'The value of the integral is ', sum

      end 
      

      double precision function fofx( x )
      double precision x
      fofx = exp(x)*cos(x)
      return
      end
-------------- next part --------------
      program gleg

****  Gauss-Laguerre Qaudrature

****  Program will integrate functon exp(x)*cos(x)

      double precision X(8), W(8), sum, fofx, A, B
      external fofx

      X(1) = -9.602898564975363E-001  
      X(2) = -7.966664774136267E-001  
      X(3) = -5.255324099163290E-001  
      X(4) = -1.834346424956498E-001  
      X(5) = 1.834346424956498E-001  
      X(6) =  5.255324099163290E-001  
      X(7) =  7.966664774136267E-001  
      X(8) =  9.602898564975363E-001  
      W(1) =  1.012285362903706E-001
      W(2) =  2.223810344533744E-001          
      W(3) =  3.137066458778874E-001         
      W(4) =  3.626837833783621E-001       
      W(5) =  3.626837833783621E-001      
      W(6) =  3.137066458778874E-001
      W(7) =  2.223810344533744E-001     
      W(8) =  1.012285362903706E-001      

      A = -1.0
      B = 1.0
      sum = 0.0
      do 10 i = 1, 8
         sum = sum + W(i) * fofx(((B-A)*X(i)+(A+B))/2.0)*(B-A)/2.0
10    continue

      print *, 'The value of the integral is ', sum

      end 
      

      double precision function fofx( x )
      double precision x
      fofx = exp(x)*cos(x)
      return
      end