program main c*********************************************************************72 c cc MAIN is the main program for EXACTNESS_PRB. c c Discussion: c c EXACTNESS_PRB tests the EXACTNESS library. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 27 May 2014 c c Author: c c John Burkardt c implicit none call timestamp ( ) write ( *, '(a)' ) '' write ( *, '(a)' ) 'EXACTNESS_PRB' write ( *, '(a)' ) ' FORTRAN77 version' write ( *, '(a)' ) ' Test the EXACTNESS library.' call test01 ( ) call test015 ( ) call test02 ( ) call test03 ( ) call test04 ( ) call test05 ( ) call test06 ( ) call test07 ( ) call test08 ( ) call test085 ( ) call test09 ( ) c c Terminate. c write ( *, '(a)' ) '' write ( *, '(a)' ) 'EXACTNESS_PRB' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) '' call timestamp ( ) stop end subroutine test01 ( ) c*********************************************************************72 c cc TEST01 tests Gauss-Legendre rules for the Legendre integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST01' write ( *, '(a)' ) & ' Gauss-Legendre rules for the Legendre integral.' write ( *, '(a)' ) ' Density function rho(x) = 1.' write ( *, '(a)' ) ' Region: -1 <= x <= +1.' write ( *, '(a)' ) ' Exactness: 2*N-1.' do n = 1, 5 call legendre_set ( n, x, w ) p_max = 2 * n call legendre_exactness ( n, x, w, p_max ) end do return end subroutine test015 ( ) c*********************************************************************72 c cc TEST015 tests Fejer Type 1 rules for the Legendre integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 27 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST015' write ( *, '(a)' ) & ' Fejer Type 1 rules for the Legendre integral.' write ( *, '(a)' ) ' Density function rho(x) = 1.' write ( *, '(a)' ) ' Region: -1 <= x <= +1.' write ( *, '(a)' ) ' Exactness: N for N odd,' write ( *, '(a)' ) ' N-1 for N even.' do n = 1, 5 call fejer1_set ( n, x, w ) if ( mod ( n, 2 ) .eq. 1 ) then p_max = n + 1 else p_max = n end if call legendre_exactness ( n, x, w, p_max ) end do return end subroutine test02 ( ) c*********************************************************************72 c cc TEST02 tests Fejer Type 2 rules for the Legendre integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST02' write ( *, '(a)' ) & ' Fejer Type 2 rules for the Legendre integral.' write ( *, '(a)' ) ' Density function rho(x) = 1.' write ( *, '(a)' ) ' Region: -1 <= x <= +1.' write ( *, '(a)' ) ' Exactness: N for N odd,' write ( *, '(a)' ) ' N-1 for N even.' do n = 1, 5 call fejer2_set ( n, x, w ) if ( mod ( n, 2 ) .eq. 1 ) then p_max = n + 1 else p_max = n end if call legendre_exactness ( n, x, w, p_max ) end do return end subroutine test03 ( ) c*********************************************************************72 c cc TEST03 tests Clenshaw-Curtis rules for the Legendre integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 22 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST03' write ( *, '(a)' ) & ' Clenshaw-Curtis rules for the Legendre integral.' write ( *, '(a)' ) ' Density function rho(x) = 1.' write ( *, '(a)' ) ' Region: -1 <= x <= +1.' write ( *, '(a)' ) ' Exactness: N for N odd,' write ( *, '(a)' ) ' N-1 for N even.' do n = 1, 5 call clenshaw_curtis_set ( n, x, w ) if ( mod ( n, 2 ) .eq. 1 ) then p_max = n + 1 else p_max = n end if call legendre_exactness ( n, x, w, p_max ) end do return end subroutine test04 ( ) c*********************************************************************72 c cc TEST04 tests Gauss-Laguerre rules for the Laguerre integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST04' write ( *, '(a)' ) & ' Gauss-Laguerre rules for the Laguerre integral.' write ( *, '(a)' ) ' Density function rho(x) = exp(-x).' write ( *, '(a)' ) ' Region: 0 <= x .lt. +oo.' write ( *, '(a)' ) ' Exactness: 2N-1.' do n = 1, 5 call laguerre_set ( n, x, w ) p_max = 2 * n call laguerre_exactness ( n, x, w, p_max ) end do return end subroutine test05 ( ) c*********************************************************************72 c cc TEST05 tests Gauss-Laguerre rules for the Laguerre integral with rho=1. c c Discussion: c c Instead of the usual density rho(x)=exp(-x), these rules apply to c approximating the integral: c I(f) = integral ( 0 <= x .lt. +oo ) f(x) dx c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer i integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST05' write ( *, '(a)' ) & ' Gauss-Laguerre quadrature rules for the Laguerre integral.' write ( *, '(a)' ) ' Density function rho(x) = 1.' write ( *, '(a)' ) ' Region: 0 <= x .lt. +oo.' write ( *, '(a)' ) ' Exactness: 2N-1.' do n = 1, 5 call laguerre_1_set ( n, x, w ) c c Standardize the rule by multiplying every weight w(i) by exp(-x(i)). c do i = 1, n w(i) = exp ( - x(i) ) * w(i) end do c c Now test the rule in standard form. c p_max = 2 * n call laguerre_exactness ( n, x, w, p_max ) end do return end subroutine test06 ( ) c*********************************************************************72 c cc TEST06 tests Gauss-Hermite rules for the Hermite integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST06' write ( *, '(a)' ) & ' Gauss-Hermite rules for the Hermite integral.' write ( *, '(a)' ) ' Density function rho(x) = exp(-x^2).' write ( *, '(a)' ) ' Region: -oo .lt. x .lt. +oo.' write ( *, '(a)' ) ' Exactness: 2N-1.' do n = 1, 5 call hermite_set ( n, x, w ) p_max = 2 * n call hermite_exactness ( n, x, w, p_max ) end do return end subroutine test07 ( ) c*********************************************************************72 c cc TEST07 tests Gauss-Hermite rules for the Hermite integral. c c Discussion: c c Instead of the usual density rho(x)=exp(-x*x), these rules apply to c approximating the integral: c I(f) = integral ( -oo .lt. x .lt. +oo ) f(x) dx c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer i integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST07' write ( *, '(a)' ) & ' Gauss-Hermite rules for the Hermite integral.' write ( *, '(a)' ) ' Density function rho(x) = 1.' write ( *, '(a)' ) ' Region: -oo .lt. x .lt. +oo.' write ( *, '(a)' ) ' Exactness: 2N-1.' do n = 1, 5 call hermite_1_set ( n, x, w ) c c Standardize the rule by multiplying every weight w(i) by exp(-x(i)^2). c do i = 1, n w(i) = exp ( - x(i) * x(i) ) * w(i) end do c c Now test the rule in standard form. c p_max = 2 * n call hermite_exactness ( n, x, w, p_max ) end do return end subroutine test08 ( ) c*********************************************************************72 c cc TEST08 tests Gauss-Chebyshev1 rules for the Chebyshev1 integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 27 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer i integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST08' write ( *, '(a)' ) & ' Gauss-Chebyshev1 rules for the Chebyshev1 integral.' write ( *, '(a)' ) ' Density function rho(x) = 1/sqrt(1-x^2).' write ( *, '(a)' ) ' Region: -1 .lt. x .lt. +1.' write ( *, '(a)' ) ' Exactness: 2N-1.' do n = 1, 5 call chebyshev1_set ( n, x, w ) p_max = 2 * n call chebyshev1_exactness ( n, x, w, p_max ) end do return end subroutine test085 ( ) c*********************************************************************72 c cc TEST085 tests Gauss-Chebyshev3 rules for the Chebyshev1 integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer i integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST085' write ( *, '(a)' ) & ' Gauss-Chebyshev3 rules for the Chebyshev1 integral.' write ( *, '(a)' ) ' Density function rho(x) = 1/sqrt(1-x^2).' write ( *, '(a)' ) ' Region: -1 .lt. x .lt. +1.' write ( *, '(a)' ) ' Exactness: 2N-3.' do n = 1, 5 call chebyshev3_set ( n, x, w ) if ( n .eq. 1 ) then p_max = 2 * n else p_max = 2 * n - 2 end if call chebyshev1_exactness ( n, x, w, p_max ) end do return end subroutine test09 ( ) c*********************************************************************72 c cc TEST09 tests Gauss-Chebyshev2 rules for the Chebyshev2 integral. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 May 2014 c c Author: c c John Burkardt c implicit none integer n_max parameter ( n_max = 10 ) integer i integer n integer p_max double precision w(n_max) double precision x(n_max) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TEST09' write ( *, '(a)' ) & ' Gauss-Chebyshev1 rules for the Chebyshev2 integral.' write ( *, '(a)' ) ' Density function rho(x) = sqrt(1-x^2).' write ( *, '(a)' ) ' Region: -1 .lt. x .lt. +1.' write ( *, '(a)' ) ' Exactness: 2N-1.' do n = 1, 5 call chebyshev2_set ( n, x, w ) p_max = 2 * n call chebyshev2_exactness ( n, x, w, p_max ) end do return end subroutine chebyshev1_set ( n, x, w ) c*********************************************************************72 c cc CHEBYSHEV1_SET sets a Chebyshev Type 1 quadrature rule. c c Discussion: c c The integral: c c integral ( -1 <= x <= 1 ) f(x) / sqrt ( 1 - x * x ) dx c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 27 May 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n == 1 ) then x( 1) = 0.0D+00 w( 1) = 3.141592653589793D+00 else if ( n == 2 ) then x( 1) = -0.7071067811865475D+00 x( 2) = 0.7071067811865476D+00 w( 1) = 1.570796326794897D+00 w( 2) = 1.570796326794897D+00 else if ( n == 3 ) then x( 1) = -0.8660254037844387D+00 x( 2) = 0.0D+00 x( 3) = 0.8660254037844387D+00 w( 1) = 1.047197551196598D+00 w( 2) = 1.047197551196598D+00 w( 3) = 1.047197551196598D+00 else if ( n == 4 ) then x( 1) = -0.9238795325112867D+00 x( 2) = -0.3826834323650897D+00 x( 3) = 0.3826834323650898D+00 x( 4) = 0.9238795325112867D+00 w( 1) = 0.7853981633974483D+00 w( 2) = 0.7853981633974483D+00 w( 3) = 0.7853981633974483D+00 w( 4) = 0.7853981633974483D+00 else if ( n == 5 ) then x( 1) = -0.9510565162951535D+00 x( 2) = -0.5877852522924730D+00 x( 3) = 0.0D+00 x( 4) = 0.5877852522924731D+00 x( 5) = 0.9510565162951535D+00 w( 1) = 0.6283185307179586D+00 w( 2) = 0.6283185307179586D+00 w( 3) = 0.6283185307179586D+00 w( 4) = 0.6283185307179586D+00 w( 5) = 0.6283185307179586D+00 else if ( n == 6 ) then x( 1) = -0.9659258262890682D+00 x( 2) = -0.7071067811865475D+00 x( 3) = -0.2588190451025206D+00 x( 4) = 0.2588190451025207D+00 x( 5) = 0.7071067811865476D+00 x( 6) = 0.9659258262890683D+00 w( 1) = 0.5235987755982988D+00 w( 2) = 0.5235987755982988D+00 w( 3) = 0.5235987755982988D+00 w( 4) = 0.5235987755982988D+00 w( 5) = 0.5235987755982988D+00 w( 6) = 0.5235987755982988D+00 else if ( n == 7 ) then x( 1) = -0.9749279121818237D+00 x( 2) = -0.7818314824680295D+00 x( 3) = -0.4338837391175581D+00 x( 4) = 0.0D+00 x( 5) = 0.4338837391175582D+00 x( 6) = 0.7818314824680298D+00 x( 7) = 0.9749279121818236D+00 w( 1) = 0.4487989505128276D+00 w( 2) = 0.4487989505128276D+00 w( 3) = 0.4487989505128276D+00 w( 4) = 0.4487989505128276D+00 w( 5) = 0.4487989505128276D+00 w( 6) = 0.4487989505128276D+00 w( 7) = 0.4487989505128276D+00 else if ( n == 8 ) then x( 1) = -0.9807852804032304D+00 x( 2) = -0.8314696123025453D+00 x( 3) = -0.5555702330196020D+00 x( 4) = -0.1950903220161282D+00 x( 5) = 0.1950903220161283D+00 x( 6) = 0.5555702330196023D+00 x( 7) = 0.8314696123025452D+00 x( 8) = 0.9807852804032304D+00 w( 1) = 0.3926990816987241D+00 w( 2) = 0.3926990816987241D+00 w( 3) = 0.3926990816987241D+00 w( 4) = 0.3926990816987241D+00 w( 5) = 0.3926990816987241D+00 w( 6) = 0.3926990816987241D+00 w( 7) = 0.3926990816987241D+00 w( 8) = 0.3926990816987241D+00 else if ( n == 9 ) then x( 1) = -0.9848077530122080D+00 x( 2) = -0.8660254037844385D+00 x( 3) = -0.6427876096865394D+00 x( 4) = -0.3420201433256685D+00 x( 5) = 0.0D+00 x( 6) = 0.3420201433256688D+00 x( 7) = 0.6427876096865394D+00 x( 8) = 0.8660254037844387D+00 x( 9) = 0.9848077530122080D+00 w( 1) = 0.3490658503988659D+00 w( 2) = 0.3490658503988659D+00 w( 3) = 0.3490658503988659D+00 w( 4) = 0.3490658503988659D+00 w( 5) = 0.3490658503988659D+00 w( 6) = 0.3490658503988659D+00 w( 7) = 0.3490658503988659D+00 w( 8) = 0.3490658503988659D+00 w( 9) = 0.3490658503988659D+00 else if ( n == 10 ) then x( 1) = -0.9876883405951377D+00 x( 2) = -0.8910065241883678D+00 x( 3) = -0.7071067811865475D+00 x( 4) = -0.4539904997395467D+00 x( 5) = -0.1564344650402306D+00 x( 6) = 0.1564344650402309D+00 x( 7) = 0.4539904997395468D+00 x( 8) = 0.7071067811865476D+00 x( 9) = 0.8910065241883679D+00 x(10) = 0.9876883405951378D+00 w( 1) = 0.3141592653589793D+00 w( 2) = 0.3141592653589793D+00 w( 3) = 0.3141592653589793D+00 w( 4) = 0.3141592653589793D+00 w( 5) = 0.3141592653589793D+00 w( 6) = 0.3141592653589793D+00 w( 7) = 0.3141592653589793D+00 w( 8) = 0.3141592653589793D+00 w( 9) = 0.3141592653589793D+00 w(10) = 0.3141592653589793D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHEBYSHEV1_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 through 10.' stop 1 end if return end subroutine chebyshev2_set ( n, x, w ) c*********************************************************************72 c cc CHEBYSHEV2_SET sets a Chebyshev Type 2 quadrature rule. c c Discussion: c c The integral: c c integral ( -1 <= x <= 1 ) f(x) * sqrt ( 1 - x * x ) dx c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 28 May 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x( 1) = 0.0D+00 w( 1) = 1.570796326794897D+00 else if ( n .eq. 2 ) then x( 1) = -0.5000000000000000D+00 x( 2) = 0.5000000000000000D+00 w( 1) = 0.7853981633974484D+00 w( 2) = 0.7853981633974481D+00 else if ( n .eq. 3 ) then x( 1) = -0.7071067811865475D+00 x( 2) = 0.0D+00 x( 3) = 0.7071067811865476D+00 w( 1) = 0.3926990816987243D+00 w( 2) = 0.7853981633974483D+00 w( 3) = 0.3926990816987240D+00 else if ( n .eq. 4 ) then x( 1) = -0.8090169943749473D+00 x( 2) = -0.3090169943749473D+00 x( 3) = 0.3090169943749475D+00 x( 4) = 0.8090169943749475D+00 w( 1) = 0.2170787134227061D+00 w( 2) = 0.5683194499747424D+00 w( 3) = 0.5683194499747423D+00 w( 4) = 0.2170787134227060D+00 else if ( n .eq. 5 ) then x( 1) = -0.8660254037844387D+00 x( 2) = -0.5000000000000000D+00 x( 3) = 0.0D+00 x( 4) = 0.5000000000000000D+00 x( 5) = 0.8660254037844387D+00 w( 1) = 0.1308996938995747D+00 w( 2) = 0.3926990816987242D+00 w( 3) = 0.5235987755982988D+00 w( 4) = 0.3926990816987240D+00 w( 5) = 0.1308996938995747D+00 else if ( n .eq. 6 ) then x( 1) = -0.9009688679024190D+00 x( 2) = -0.6234898018587335D+00 x( 3) = -0.2225209339563143D+00 x( 4) = 0.2225209339563144D+00 x( 5) = 0.6234898018587336D+00 x( 6) = 0.9009688679024191D+00 w( 1) = 0.08448869089158863D+00 w( 2) = 0.2743330560697779D+00 w( 3) = 0.4265764164360819D+00 w( 4) = 0.4265764164360819D+00 w( 5) = 0.2743330560697778D+00 w( 6) = 0.08448869089158857D+00 else if ( n .eq. 7 ) then x( 1) = -0.9238795325112867D+00 x( 2) = -0.7071067811865475D+00 x( 3) = -0.3826834323650897D+00 x( 4) = 0.0D+00 x( 5) = 0.3826834323650898D+00 x( 6) = 0.7071067811865476D+00 x( 7) = 0.9238795325112867D+00 w( 1) = 0.05750944903191316D+00 w( 2) = 0.1963495408493621D+00 w( 3) = 0.3351896326668110D+00 w( 4) = 0.3926990816987241D+00 w( 5) = 0.3351896326668110D+00 w( 6) = 0.1963495408493620D+00 w( 7) = 0.05750944903191313D+00 else if ( n .eq. 8 ) then x( 1) = -0.9396926207859083D+00 x( 2) = -0.7660444431189779D+00 x( 3) = -0.5000000000000000D+00 x( 4) = -0.1736481776669303D+00 x( 5) = 0.1736481776669304D+00 x( 6) = 0.5000000000000000D+00 x( 7) = 0.7660444431189780D+00 x( 8) = 0.9396926207859084D+00 w( 1) = 0.04083294770910712D+00 w( 2) = 0.1442256007956728D+00 w( 3) = 0.2617993877991495D+00 w( 4) = 0.3385402270935190D+00 w( 5) = 0.3385402270935190D+00 w( 6) = 0.2617993877991494D+00 w( 7) = 0.1442256007956727D+00 w( 8) = 0.04083294770910708D+00 else if ( n .eq. 9 ) then x( 1) = -0.9510565162951535D+00 x( 2) = -0.8090169943749473D+00 x( 3) = -0.5877852522924730D+00 x( 4) = -0.3090169943749473D+00 x( 5) = 0.0D+00 x( 6) = 0.3090169943749475D+00 x( 7) = 0.5877852522924731D+00 x( 8) = 0.8090169943749475D+00 x( 9) = 0.9510565162951535D+00 w( 1) = 0.02999954037160818D+00 w( 2) = 0.1085393567113530D+00 w( 3) = 0.2056199086476263D+00 w( 4) = 0.2841597249873712D+00 w( 5) = 0.3141592653589793D+00 w( 6) = 0.2841597249873711D+00 w( 7) = 0.2056199086476263D+00 w( 8) = 0.1085393567113530D+00 w( 9) = 0.02999954037160816D+00 else if ( n .eq. 10 ) then x( 1) = -0.9594929736144974D+00 x( 2) = -0.8412535328311811D+00 x( 3) = -0.6548607339452850D+00 x( 4) = -0.4154150130018863D+00 x( 5) = -0.1423148382732850D+00 x( 6) = 0.1423148382732851D+00 x( 7) = 0.4154150130018864D+00 x( 8) = 0.6548607339452851D+00 x( 9) = 0.8412535328311812D+00 x(10) = 0.9594929736144974D+00 w( 1) = 0.02266894250185884D+00 w( 2) = 0.08347854093418908D+00 w( 3) = 0.1631221774548166D+00 w( 4) = 0.2363135602034873D+00 w( 5) = 0.2798149423030966D+00 w( 6) = 0.2798149423030965D+00 w( 7) = 0.2363135602034873D+00 w( 8) = 0.1631221774548166D+00 w( 9) = 0.08347854093418902D+00 w(10) = 0.02266894250185884D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHEBYSHEV2_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 through 10.' stop 1 end if return end subroutine chebyshev3_set ( n, x, w ) !*********************************************************************72 ! !! CHEBYSHEV3_SET sets a Chebyshev Type 3 quadrature rule. ! ! Discussion: ! ! The integral: ! ! integral ( -1 <= x <= 1 ) f(x) / sqrt ( 1 - x * x ) dx ! ! The quadrature rule: ! ! sum ( 1 <= i <= n ) w(i) * f ( x(i) ) ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 28 May 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the order. ! N must be between 1 and 10. ! ! Output, double precision X(N), the abscissas. ! ! Output, double precision W(N), the weights. ! implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x( 1) = 0.000000000000000D+00 w( 1) = 3.141592653589793D+00 else if ( n .eq. 2 ) then x( 1) = -1.000000000000000D+00 x( 2) = 1.000000000000000D+00 w( 1) = 1.570796326794897D+00 w( 2) = 1.570796326794897D+00 else if ( n .eq. 3 ) then x( 1) = -1.000000000000000D+00 x( 2) = 0.0D+00 x( 3) = 1.000000000000000D+00 w( 1) = 0.7853981633974483D+00 w( 2) = 1.570796326794897D+00 w( 3) = 0.7853981633974483D+00 else if ( n .eq. 4 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.5000000000000000D+00 x( 3) = 0.5000000000000000D+00 x( 4) = 1.000000000000000D+00 w( 1) = 0.5235987755982988D+00 w( 2) = 1.047197551196598D+00 w( 3) = 1.047197551196598D+00 w( 4) = 0.5235987755982988D+00 else if ( n .eq. 5 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.7071067811865475D+00 x( 3) = 0.0D+00 x( 4) = 0.7071067811865476D+00 x( 5) = 1.000000000000000D+00 w( 1) = 0.3926990816987241D+00 w( 2) = 0.7853981633974483D+00 w( 3) = 0.7853981633974483D+00 w( 4) = 0.7853981633974483D+00 w( 5) = 0.3926990816987241D+00 else if ( n .eq. 6 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.8090169943749473D+00 x( 3) = -0.3090169943749473D+00 x( 4) = 0.3090169943749475D+00 x( 5) = 0.8090169943749475D+00 x( 6) = 1.000000000000000D+00 w( 1) = 0.3141592653589793D+00 w( 2) = 0.6283185307179586D+00 w( 3) = 0.6283185307179586D+00 w( 4) = 0.6283185307179586D+00 w( 5) = 0.6283185307179586D+00 w( 6) = 0.3141592653589793D+00 else if ( n .eq. 7 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.8660254037844387D+00 x( 3) = -0.5000000000000000D+00 x( 4) = 0.0D+00 x( 5) = 0.5000000000000001D+00 x( 6) = 0.8660254037844387D+00 x( 7) = 1.000000000000000D+00 w( 1) = 0.2617993877991494D+00 w( 2) = 0.5235987755982988D+00 w( 3) = 0.5235987755982988D+00 w( 4) = 0.5235987755982988D+00 w( 5) = 0.5235987755982988D+00 w( 6) = 0.5235987755982988D+00 w( 7) = 0.2617993877991494D+00 else if ( n .eq. 8 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.9009688679024190D+00 x( 3) = -0.6234898018587335D+00 x( 4) = -0.2225209339563143D+00 x( 5) = 0.2225209339563144D+00 x( 6) = 0.6234898018587336D+00 x( 7) = 0.9009688679024191D+00 x( 8) = 1.000000000000000D+00 w( 1) = 0.2243994752564138D+00 w( 2) = 0.4487989505128276D+00 w( 3) = 0.4487989505128276D+00 w( 4) = 0.4487989505128276D+00 w( 5) = 0.4487989505128276D+00 w( 6) = 0.4487989505128276D+00 w( 7) = 0.4487989505128276D+00 w( 8) = 0.2243994752564138D+00 else if ( n .eq. 9 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.9238795325112867D+00 x( 3) = -0.7071067811865475D+00 x( 4) = -0.3826834323650897D+00 x( 5) = 0.0D+00 x( 6) = 0.3826834323650898D+00 x( 7) = 0.7071067811865476D+00 x( 8) = 0.9238795325112867D+00 x( 9) = 1.000000000000000D+00 w( 1) = 0.1963495408493621D+00 w( 2) = 0.3926990816987241D+00 w( 3) = 0.3926990816987241D+00 w( 4) = 0.3926990816987241D+00 w( 5) = 0.3926990816987241D+00 w( 6) = 0.3926990816987241D+00 w( 7) = 0.3926990816987241D+00 w( 8) = 0.3926990816987241D+00 w( 9) = 0.1963495408493621D+00 else if ( n .eq. 10 ) then x( 1) = -1.000000000000000D+00 x( 2) = -0.9396926207859083D+00 x( 3) = -0.7660444431189779D+00 x( 4) = -0.5000000000000000D+00 x( 5) = -0.1736481776669303D+00 x( 6) = 0.1736481776669304D+00 x( 7) = 0.5000000000000001D+00 x( 8) = 0.7660444431189780D+00 x( 9) = 0.9396926207859084D+00 x(10) = 1.000000000000000D+00 w( 1) = 0.1745329251994329D+00 w( 2) = 0.3490658503988659D+00 w( 3) = 0.3490658503988659D+00 w( 4) = 0.3490658503988659D+00 w( 5) = 0.3490658503988659D+00 w( 6) = 0.3490658503988659D+00 w( 7) = 0.3490658503988659D+00 w( 8) = 0.3490658503988659D+00 w( 9) = 0.3490658503988659D+00 w(10) = 0.1745329251994329D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CHEBYSHEV3_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 through 10.' stop 1 end if return end subroutine clenshaw_curtis_set ( n, x, w ) c*********************************************************************72 c cc CLENSHAW_CURTIS_SET sets a Clenshaw-Curtis quadrature rule. c c Discussion: c c The integral: c c integral ( -1 <= x <= 1 ) f(x) dx c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c The abscissas for the rule of order N can be regarded c as the cosines of equally spaced angles between 180 and 0 degrees: c c X(I) = cos ( ( I - 1 ) * PI / ( N - 1 ) ) c c except for the basic case N = 1, when c c X(1) = 0. c c A Clenshaw-Curtis rule that uses N points will integrate c exactly all polynomials of degrees 0 through N-1. If N c is odd, then by symmetry the polynomial of degree N will c also be integrated exactly. c c If the value of N is increased in a sensible way, then c the new set of abscissas will include the old ones. One such c sequence would be N(K) = 2*K+1 for K = 0, 1, 2, ... c Thus, in the table below, the abscissas for order 9 include c those for order 5. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 19 April 2011 c c Author: c c John Burkardt c c Reference: c c Charles Clenshaw, Alan Curtis, c A Method for Numerical Integration on an Automatic Computer, c Numerische Mathematik, c Volume 2, Number 1, December 1960, pages 197-205. c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 0.00000000000000000000D+00 w(1) = 2.00000000000000000000D+00 else if ( n .eq. 2 ) then x(1) = -1.00000000000000000000D+00 x(2) = 1.00000000000000000000D+00 w(1) = 1.00000000000000000000D+00 w(2) = 1.00000000000000000000D+00 else if ( n .eq. 3 ) then x(1) = -1.00000000000000000000D+00 x(2) = 0.00000000000000000000D+00 x(3) = 1.00000000000000000000D+00 w(1) = 0.33333333333333333333D+00 w(2) = 1.33333333333333333333D+00 w(3) = 0.33333333333333333333D+00 else if ( n .eq. 4 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.50000000000000000000D+00 x(3) = 0.50000000000000000000D+00 x(4) = 1.00000000000000000000D+00 w(1) = 0.11111111111111111111D+00 w(2) = 0.88888888888888888889D+00 w(3) = 0.88888888888888888889D+00 w(4) = 0.11111111111111111111D+00 else if ( n .eq. 5 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.70710678118654752440D+00 x(3) = 0.00000000000000000000D+00 x(4) = 0.70710678118654752440D+00 x(5) = 1.00000000000000000000D+00 w(1) = 0.06666666666666666667D+00 w(2) = 0.53333333333333333333D+00 w(3) = 0.80000000000000000000D+00 w(4) = 0.53333333333333333333D+00 w(5) = 0.06666666666666666667D+00 else if ( n .eq. 6 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.80901699437494742410D+00 x(3) = -0.30901699437494742410D+00 x(4) = 0.30901699437494742410D+00 x(5) = 0.80901699437493732410D+00 x(6) = 1.00000000000000000000D+00 w(1) = 0.04000000000000000000D+00 w(2) = 0.36074304120001121619D+00 w(3) = 0.59925695879998878381D+00 w(4) = 0.59925695879998878381D+00 w(5) = 0.36074304120001121619D+00 w(6) = 0.04000000000000000000D+00 else if ( n .eq. 7 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.86602540378443864676D+00 x(3) = -0.50000000000000000000D+00 x(4) = 0.00000000000000000000D+00 x(5) = 0.50000000000000000000D+00 x(6) = 0.86602540378443864676D+00 x(7) = 1.00000000000000000000D+00 w(1) = 0.02857142857142857143D+00 w(2) = 0.25396825396825396825D+00 w(3) = 0.45714285714285714286D+00 w(4) = 0.52063492063492063492D+00 w(5) = 0.45714285714285714286D+00 w(6) = 0.25396825396825396825D+00 w(7) = 0.02857142857142857143D+00 else if ( n .eq. 8 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.90096886790241912624D+00 x(3) = -0.62348980185873353053D+00 x(4) = -0.22252093395631440429D+00 x(5) = 0.22252093395631440429D+00 x(6) = 0.62348980185873353053D+00 x(7) = 0.90096886790241910624D+00 x(8) = 1.00000000000000000000D+00 w(1) = 0.02040816326530612245D+00 w(2) = 0.19014100721820835178D+00 w(3) = 0.35224242371815911533D+00 w(4) = 0.43720840579832641044D+00 w(5) = 0.43720840579832641044D+00 w(6) = 0.35224242371815911533D+00 w(7) = 0.19014100721820835178D+00 w(8) = 0.02040816326530612245D+00 else if ( n .eq. 9 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.92387953251128675613D+00 x(3) = -0.70710678118654752440D+00 x(4) = -0.38268343236508977173D+00 x(5) = 0.00000000000000000000D+00 x(6) = 0.38268343236508977173D+00 x(7) = 0.70710678118654752440D+00 x(8) = 0.92387953251128675613D+00 x(9) = 1.00000000000000000000D+00 w(1) = 0.01587301587301587302D+00 w(2) = 0.14621864921601815501D+00 w(3) = 0.27936507936507936508D+00 w(4) = 0.36171785872048978150D+00 w(5) = 0.39365079365079365079D+00 w(6) = 0.36171785872048978150D+00 w(7) = 0.27936507936507936508D+00 w(8) = 0.14621864921601815501D+00 w(9) = 0.01587301587301587302D+00 else if ( n .eq. 10 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.93969262078590838405D+00 x(3) = -0.76604444311897903520D+00 x(4) = -0.50000000000000000000D+00 x(5) = -0.17364817766693034885D+00 x(6) = 0.17364817766693034885D+00 x(7) = 0.50000000000000000000D+00 x(8) = 0.76604444311897903520D+00 x(9) = 0.93969262078590838405D+00 x(10) = 1.00000000000000000000D+00 w(1) = 0.01234567901234567901D+00 w(2) = 0.11656745657203712296D+00 w(3) = 0.22528432333810440813D+00 w(4) = 0.30194003527336860670D+00 w(5) = 0.34386250580414418320D+00 w(6) = 0.34386250580414418320D+00 w(7) = 0.30194003527336860670D+00 w(8) = 0.22528432333810440813D+00 w(9) = 0.11656745657203712296D+00 w(10) = 0.01234567901234567901D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'CLENSHAW_CURTIS_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 to 10.' stop 1 end if return end subroutine fejer1_set ( n, x, w ) c*********************************************************************72 c cc FEJER1_SET sets abscissas and weights for Fejer type 1 quadrature. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 19 April 2011 c c Author: c c John Burkardt c c Reference: c c Philip Davis, Philip Rabinowitz, c Methods of Numerical Integration, c Second Edition, c Dover, 2007, c ISBN: 0486453391, c LC: QA299.3.D28. c c Walter Gautschi, c Numerical Quadrature in the Presence of a Singularity, c SIAM Journal on Numerical Analysis, c Volume 4, Number 3, 1967, pages 357-362. c c Joerg Waldvogel, c Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules, c BIT Numerical Mathematics, c Volume 43, Number 1, 2003, pages 1-18. c c Parameters: c c Input, integer N, the order. c N should be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 0.000000000000000D+00 w(1) = 2.000000000000000D+00 else if ( n .eq. 2 ) then x(1) = -0.7071067811865475D+00 x(2) = 0.7071067811865475D+00 w(1) = 1.000000000000000D+00 w(2) = 1.000000000000000D+00 else if ( n .eq. 3 ) then x( 1) = -0.8660254037844387D+00 x( 2) = 0.0000000000000000D+00 x( 3) = 0.8660254037844387D+00 w(1) = 0.4444444444444444D+00 w(2) = 1.1111111111111111D+00 w(3) = 0.4444444444444444D+00 else if ( n .eq. 4 ) then x( 1) = -0.9238795325112867D+00 x( 2) = -0.3826834323650897D+00 x( 3) = 0.3826834323650898D+00 x( 4) = 0.9238795325112867D+00 w( 1) = 0.2642977396044841D+00 w( 2) = 0.7357022603955158D+00 w( 3) = 0.7357022603955158D+00 w( 4) = 0.2642977396044841D+00 else if ( n .eq. 5 ) then x( 1) = -0.9510565162951535D+00 x( 2) = -0.5877852522924730D+00 x( 3) = 0.0000000000000000D+00 x( 4) = 0.5877852522924731D+00 x( 5) = 0.9510565162951535D+00 w( 1) = 0.1677812284666835D+00 w( 2) = 0.5255521048666498D+00 w( 3) = 0.6133333333333333D+00 w( 4) = 0.5255521048666498D+00 w( 5) = 0.1677812284666835D+00 else if ( n .eq. 6 ) then x( 1) = -0.9659258262890682D+00 x( 2) = -0.7071067811865475D+00 x( 3) = -0.2588190451025206D+00 x( 4) = 0.2588190451025207D+00 x( 5) = 0.7071067811865476D+00 x( 6) = 0.9659258262890683D+00 w( 1) = 0.1186610213812358D+00 w( 2) = 0.3777777777777778D+00 w( 3) = 0.5035612008409863D+00 w( 4) = 0.5035612008409863D+00 w( 5) = 0.3777777777777778D+00 w( 6) = 0.1186610213812358D+00 else if ( n .eq. 7 ) then x( 1) = -0.9749279121818237D+00 x( 2) = -0.7818314824680295D+00 x( 3) = -0.4338837391175581D+00 x( 4) = 0.0000000000000000D+00 x( 5) = 0.4338837391175582D+00 x( 6) = 0.7818314824680298D+00 x( 7) = 0.9749279121818236D+00 w( 1) = 0.08671618072672234D+00 w( 2) = 0.2878313947886919D+00 w( 3) = 0.3982415401308441D+00 w( 4) = 0.4544217687074830D+00 w( 5) = 0.3982415401308441D+00 w( 6) = 0.2878313947886919D+00 w( 7) = 0.08671618072672234D+00 else if ( n .eq. 8 ) then x( 1) = -0.9807852804032304D+00 x( 2) = -0.8314696123025453D+00 x( 3) = -0.5555702330196020D+00 x( 4) = -0.1950903220161282D+00 x( 5) = 0.1950903220161283D+00 x( 6) = 0.5555702330196023D+00 x( 7) = 0.8314696123025452D+00 x( 8) = 0.9807852804032304D+00 w( 1) = 0.06698294569858981D+00 w( 2) = 0.2229879330145788D+00 w( 3) = 0.3241525190645244D+00 w( 4) = 0.3858766022223071D+00 w( 5) = 0.3858766022223071D+00 w( 6) = 0.3241525190645244D+00 w( 7) = 0.2229879330145788D+00 w( 8) = 0.06698294569858981D+00 else if ( n .eq. 9 ) then x( 1) = -0.9848077530122080D+00 x( 2) = -0.8660254037844385D+00 x( 3) = -0.6427876096865394D+00 x( 4) = -0.3420201433256685D+00 x( 5) = 0.0000000000000000D+00 x( 6) = 0.3420201433256688D+00 x( 7) = 0.6427876096865394D+00 x( 8) = 0.8660254037844387D+00 x( 9) = 0.9848077530122080D+00 w( 1) = 0.05273664990990676D+00 w( 2) = 0.1791887125220458D+00 w( 3) = 0.2640372225410044D+00 w( 4) = 0.3308451751681364D+00 w( 5) = 0.3463844797178130D+00 w( 6) = 0.3308451751681364D+00 w( 7) = 0.2640372225410044D+00 w( 8) = 0.1791887125220458D+00 w( 9) = 0.05273664990990676D+00 else if ( n .eq. 10 ) then x( 1) = -0.9876883405951377D+00 x( 2) = -0.8910065241883678D+00 x( 3) = -0.7071067811865475D+00 x( 4) = -0.4539904997395467D+00 x( 5) = -0.1564344650402306D+00 x( 6) = 0.1564344650402309D+00 x( 7) = 0.4539904997395468D+00 x( 8) = 0.7071067811865476D+00 x( 9) = 0.8910065241883679D+00 x(10) = 0.9876883405951378D+00 w( 1) = 0.04293911957413078D+00 w( 2) = 0.1458749193773909D+00 w( 3) = 0.2203174603174603D+00 w( 4) = 0.2808792186638755D+00 w( 5) = 0.3099892820671425D+00 w( 6) = 0.3099892820671425D+00 w( 7) = 0.2808792186638755D+00 w( 8) = 0.2203174603174603D+00 w( 9) = 0.1458749193773909D+00 w(10) = 0.04293911957413078D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEJER1_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 through 10.' stop 1 end if return end subroutine fejer2_set ( n, x, w ) c*********************************************************************72 c cc FEJER2_SET sets abscissas and weights for Fejer type 2 quadrature. c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 19 April 2011 c c Author: c c John Burkardt c c Reference: c c Philip Davis, Philip Rabinowitz, c Methods of Numerical Integration, c Second Edition, c Dover, 2007, c ISBN: 0486453391, c LC: QA299.3.D28. c c Walter Gautschi, c Numerical Quadrature in the Presence of a Singularity, c SIAM Journal on Numerical Analysis, c Volume 4, Number 3, 1967, pages 357-362. c c Joerg Waldvogel, c Fast Construction of the Fejer and Clenshaw-Curtis Quadrature Rules, c BIT Numerical Mathematics, c Volume 43, Number 1, 2003, pages 1-18. c c Parameters: c c Input, integer N, the order. c N should be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 0.000000000000000D+00 w(1) = 2.000000000000000D+00 else if ( n .eq. 2 ) then x(1) = -0.5000000000000000D+00 x(2) = 0.5000000000000000D+00 w(1) = 1.0000000000000000D+00 w(2) = 1.0000000000000000D+00 else if ( n .eq. 3 ) then x( 1) = -0.7071067811865476D+00 x( 2) = 0.0000000000000000D+00 x( 3) = 0.7071067811865476D+00 w(1) = 0.6666666666666666D+00 w(2) = 0.6666666666666666D+00 w(3) = 0.6666666666666666D+00 else if ( n .eq. 4 ) then x( 1) = -0.8090169943749475D+00 x( 2) = -0.3090169943749475D+00 x( 3) = 0.3090169943749475D+00 x( 4) = 0.8090169943749475D+00 w( 1) = 0.4254644007500070D+00 w( 2) = 0.5745355992499930D+00 w( 3) = 0.5745355992499930D+00 w( 4) = 0.4254644007500070D+00 else if ( n .eq. 5 ) then x( 1) = -0.8660254037844387D+00 x( 2) = -0.5000000000000000D+00 x( 3) = 0.0000000000000000D+00 x( 4) = 0.5000000000000000D+00 x( 5) = 0.8660254037844387D+00 w( 1) = 0.3111111111111111D+00 w( 2) = 0.4000000000000000D+00 w( 3) = 0.5777777777777777D+00 w( 4) = 0.4000000000000000D+00 w( 5) = 0.3111111111111111D+00 else if ( n .eq. 6 ) then x( 1) = -0.9009688679024191D+00 x( 2) = -0.6234898018587336D+00 x( 3) = -0.2225209339563144D+00 x( 4) = 0.2225209339563144D+00 x( 5) = 0.6234898018587336D+00 x( 6) = 0.9009688679024191D+00 w( 1) = 0.2269152467244296D+00 w( 2) = 0.3267938603769863D+00 w( 3) = 0.4462908928985841D+00 w( 4) = 0.4462908928985841D+00 w( 5) = 0.3267938603769863D+00 w( 6) = 0.2269152467244296D+00 else if ( n .eq. 7 ) then x( 1) = -0.9238795325112867D+00 x( 2) = -0.7071067811865476D+00 x( 3) = -0.3826834323650898D+00 x( 4) = 0.0000000000000000D+00 x( 5) = 0.3826834323650898D+00 x( 6) = 0.7071067811865476D+00 x( 7) = 0.9238795325112867D+00 w( 1) = 0.1779646809620499D+00 w( 2) = 0.2476190476190476D+00 w( 3) = 0.3934638904665215D+00 w( 4) = 0.3619047619047619D+00 w( 5) = 0.3934638904665215D+00 w( 6) = 0.2476190476190476D+00 w( 7) = 0.1779646809620499D+00 else if ( n .eq. 8 ) then x( 1) = -0.9396926207859084D+00 x( 2) = -0.7660444431189780D+00 x( 3) = -0.5000000000000000D+00 x( 4) = -0.1736481776669304D+00 x( 5) = 0.1736481776669304D+00 x( 6) = 0.5000000000000000D+00 x( 7) = 0.7660444431189780D+00 x( 8) = 0.9396926207859084D+00 w( 1) = 0.1397697435050225D+00 w( 2) = 0.2063696457302284D+00 w( 3) = 0.3142857142857143D+00 w( 4) = 0.3395748964790348D+00 w( 5) = 0.3395748964790348D+00 w( 6) = 0.3142857142857143D+00 w( 7) = 0.2063696457302284D+00 w( 8) = 0.1397697435050225D+00 else if ( n .eq. 9 ) then x( 1) = -0.9510565162951535D+00 x( 2) = -0.8090169943749475D+00 x( 3) = -0.5877852522924731D+00 x( 4) = -0.3090169943749475D+00 x( 5) = 0.0000000000000000D+00 x( 6) = 0.3090169943749475D+00 x( 7) = 0.5877852522924731D+00 x( 8) = 0.8090169943749475D+00 x( 9) = 0.9510565162951535D+00 w( 1) = 0.1147810750857217D+00 w( 2) = 0.1654331942222276D+00 w( 3) = 0.2737903534857068D+00 w( 4) = 0.2790112502222170D+00 w( 5) = 0.3339682539682539D+00 w( 6) = 0.2790112502222170D+00 w( 7) = 0.2737903534857068D+00 w( 8) = 0.1654331942222276D+00 w( 9) = 0.1147810750857217D+00 else if ( n .eq. 10 ) then x( 1) = -0.9594929736144974D+00 x( 2) = -0.8412535328311812D+00 x( 3) = -0.6548607339452851D+00 x( 4) = -0.4154150130018864D+00 x( 5) = -0.1423148382732851D+00 x( 6) = 0.1423148382732851D+00 x( 7) = 0.4154150130018864D+00 x( 8) = 0.6548607339452851D+00 x( 9) = 0.8412535328311812D+00 x(10) = 0.9594929736144974D+00 w( 1) = 0.09441954173982806D+00 w( 2) = 0.1411354380109716D+00 w( 3) = 0.2263866903636005D+00 w( 4) = 0.2530509772156453D+00 w( 5) = 0.2850073526699544D+00 w( 6) = 0.2850073526699544D+00 w( 7) = 0.2530509772156453D+00 w( 8) = 0.2263866903636005D+00 w( 9) = 0.1411354380109716D+00 w(10) = 0.09441954173982806D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEJER2_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 through 10.' stop 1 end if return end subroutine hermite_set ( n, x, w ) c*********************************************************************72 c cc HERMITE_SET sets abscissas and weights for Hermite quadrature. c c Discussion: c c The integral: c c integral ( -oo .lt. x .lt. +oo ) f(x) * rho(x) dx c c The weight function: c c rho(x) = exp ( - x * x ) c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c Mathematica can numerically estimate the abscissas of the rule c of order N to P digits by the command: c c NSolve [ HermiteH [ n, x ] .eq. 0, x, p ] c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 22 April 2010 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Vladimir Krylov, c Approximate Calculation of Integrals, c Dover, 2006, c ISBN: 0486445798, c LC: QA311.K713. c c Arthur Stroud, Don Secrest, c Gaussian Quadrature Formulas, c Prentice Hall, 1966, c LC: QA299.4G3S7. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Daniel Zwillinger, editor, c CRC Standard Mathematical Tables and Formulae, c 30th Edition, c CRC Press, 1996, c ISBN: 0-8493-2479-3, c LC: QA47.M315. c c Parameters: c c Input, integer N, the order. c N must be between 1 and 20, 31/32/33, 63/64/65, 127/128/129. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 0.0D+00 w(1) = 1.77245385090551602729816748334D+00 else if ( n .eq. 2 ) then x(1) = - 0.707106781186547524400844362105D+00 x(2) = 0.707106781186547524400844362105D+00 w(1) = 0.886226925452758013649083741671D+00 w(2) = 0.886226925452758013649083741671D+00 else if ( n .eq. 3 ) then x(1) = - 0.122474487139158904909864203735D+01 x(2) = 0.0D+00 x(3) = 0.122474487139158904909864203735D+01 w(1) = 0.295408975150919337883027913890D+00 w(2) = 0.118163590060367735153211165556D+01 w(3) = 0.295408975150919337883027913890D+00 else if ( n .eq. 4 ) then x(1) = - 0.165068012388578455588334111112D+01 x(2) = - 0.524647623275290317884060253835D+00 x(3) = 0.524647623275290317884060253835D+00 x(4) = 0.165068012388578455588334111112D+01 w(1) = 0.813128354472451771430345571899D-01 w(2) = 0.804914090005512836506049184481D+00 w(3) = 0.804914090005512836506049184481D+00 w(4) = 0.813128354472451771430345571899D-01 else if ( n .eq. 5 ) then x(1) = - 0.202018287045608563292872408814D+01 x(2) = - 0.958572464613818507112770593893D+00 x(3) = 0.0D+00 x(4) = 0.958572464613818507112770593893D+00 x(5) = 0.202018287045608563292872408814D+01 w(1) = 0.199532420590459132077434585942D-01 w(2) = 0.393619323152241159828495620852D+00 w(3) = 0.945308720482941881225689324449D+00 w(4) = 0.393619323152241159828495620852D+00 w(5) = 0.199532420590459132077434585942D-01 else if ( n .eq. 6 ) then x(1) = - 0.235060497367449222283392198706D+01 x(2) = - 0.133584907401369694971489528297D+01 x(3) = - 0.436077411927616508679215948251D+00 x(4) = 0.436077411927616508679215948251D+00 x(5) = 0.133584907401369694971489528297D+01 x(6) = 0.235060497367449222283392198706D+01 w(1) = 0.453000990550884564085747256463D-02 w(2) = 0.157067320322856643916311563508D+00 w(3) = 0.724629595224392524091914705598D+00 w(4) = 0.724629595224392524091914705598D+00 w(5) = 0.157067320322856643916311563508D+00 w(6) = 0.453000990550884564085747256463D-02 else if ( n .eq. 7 ) then x(1) = - 0.265196135683523349244708200652D+01 x(2) = - 0.167355162876747144503180139830D+01 x(3) = - 0.816287882858964663038710959027D+00 x(4) = 0.0D+00 x(5) = 0.816287882858964663038710959027D+00 x(6) = 0.167355162876747144503180139830D+01 x(7) = 0.265196135683523349244708200652D+01 w(1) = 0.971781245099519154149424255939D-03 w(2) = 0.545155828191270305921785688417D-01 w(3) = 0.425607252610127800520317466666D+00 w(4) = 0.810264617556807326764876563813D+00 w(5) = 0.425607252610127800520317466666D+00 w(6) = 0.545155828191270305921785688417D-01 w(7) = 0.971781245099519154149424255939D-03 else if ( n .eq. 8 ) then x(1) = - 0.293063742025724401922350270524D+01 x(2) = - 0.198165675669584292585463063977D+01 x(3) = - 0.115719371244678019472076577906D+01 x(4) = - 0.381186990207322116854718885584D+00 x(5) = 0.381186990207322116854718885584D+00 x(6) = 0.115719371244678019472076577906D+01 x(7) = 0.198165675669584292585463063977D+01 x(8) = 0.293063742025724401922350270524D+01 w(1) = 0.199604072211367619206090452544D-03 w(2) = 0.170779830074134754562030564364D-01 w(3) = 0.207802325814891879543258620286D+00 w(4) = 0.661147012558241291030415974496D+00 w(5) = 0.661147012558241291030415974496D+00 w(6) = 0.207802325814891879543258620286D+00 w(7) = 0.170779830074134754562030564364D-01 w(8) = 0.199604072211367619206090452544D-03 else if ( n .eq. 9 ) then x(1) = - 0.319099320178152760723004779538D+01 x(2) = - 0.226658058453184311180209693284D+01 x(3) = - 0.146855328921666793166701573925D+01 x(4) = - 0.723551018752837573322639864579D+00 x(5) = 0.0D+00 x(6) = 0.723551018752837573322639864579D+00 x(7) = 0.146855328921666793166701573925D+01 x(8) = 0.226658058453184311180209693284D+01 x(9) = 0.319099320178152760723004779538D+01 w(1) = 0.396069772632643819045862946425D-04 w(2) = 0.494362427553694721722456597763D-02 w(3) = 0.884745273943765732879751147476D-01 w(4) = 0.432651559002555750199812112956D+00 w(5) = 0.720235215606050957124334723389D+00 w(6) = 0.432651559002555750199812112956D+00 w(7) = 0.884745273943765732879751147476D-01 w(8) = 0.494362427553694721722456597763D-02 w(9) = 0.396069772632643819045862946425D-04 else if ( n .eq. 10 ) then x(1) = - 0.343615911883773760332672549432D+01 x(2) = - 0.253273167423278979640896079775D+01 x(3) = - 0.175668364929988177345140122011D+01 x(4) = - 0.103661082978951365417749191676D+01 x(5) = - 0.342901327223704608789165025557D+00 x(6) = 0.342901327223704608789165025557D+00 x(7) = 0.103661082978951365417749191676D+01 x(8) = 0.175668364929988177345140122011D+01 x(9) = 0.253273167423278979640896079775D+01 x(10) = 0.343615911883773760332672549432D+01 w(1) = 0.764043285523262062915936785960D-05 w(2) = 0.134364574678123269220156558585D-02 w(3) = 0.338743944554810631361647312776D-01 w(4) = 0.240138611082314686416523295006D+00 w(5) = 0.610862633735325798783564990433D+00 w(6) = 0.610862633735325798783564990433D+00 w(7) = 0.240138611082314686416523295006D+00 w(8) = 0.338743944554810631361647312776D-01 w(9) = 0.134364574678123269220156558585D-02 w(10) = 0.764043285523262062915936785960D-05 else write ( *, '(a)' ) '' write ( *, '(a)' ) 'HERMITE_SET - Fatal error!' write ( *, '(a)' ) ' Illegal value of N.' stop 1 end if return end subroutine hermite_1_set ( n, x, w ) c*********************************************************************72 c cc HERMITE_1_SET sets abscissas and weights for Hermite quadrature. c c Discussion: c c This routine is for the case with unit density: c integral ( -oo .lt. x .lt. +oo ) f(x) dx c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 0.0D+00 w(1) = 1.7724538509055161D+00 else if ( n .eq. 2 ) then x(1) = - 0.707106781186547524400844362105D+00 x(2) = 0.707106781186547524400844362105D+00 w(1) = 1.4611411826611391D+00 w(2) = 1.4611411826611391D+00 else if ( n .eq. 3 ) then x(1) = - 0.122474487139158904909864203735D+01 x(2) = 0.0D+00 x(3) = 0.122474487139158904909864203735D+01 w(1) = 1.3239311752136438D+00 w(2) = 1.1816359006036774D+00 w(3) = 1.3239311752136438D+00 else if ( n .eq. 4 ) then x(1) = - 0.165068012388578455588334111112D+01 x(2) = - 0.524647623275290317884060253835D+00 x(3) = 0.524647623275290317884060253835D+00 x(4) = 0.165068012388578455588334111112D+01 w(1) = 1.2402258176958150D+00 w(2) = 1.0599644828949693D+00 w(3) = 1.0599644828949693D+00 w(4) = 1.2402258176958150D+00 else if ( n .eq. 5 ) then x(1) = - 0.202018287045608563292872408814D+01 x(2) = - 0.958572464613818507112770593893D+00 x(3) = 0.0D+00 x(4) = 0.958572464613818507112770593893D+00 x(5) = 0.202018287045608563292872408814D+01 w(1) = 1.1814886255359869D+00 w(2) = 0.98658099675142830D+00 w(3) = 0.94530872048294190D+00 w(4) = 0.98658099675142830D+00 w(5) = 1.1814886255359869D+00 else if ( n .eq. 6 ) then x(1) = - 0.235060497367449222283392198706D+01 x(2) = - 0.133584907401369694971489528297D+01 x(3) = - 0.436077411927616508679215948251D+00 x(4) = 0.436077411927616508679215948251D+00 x(5) = 0.133584907401369694971489528297D+01 x(6) = 0.235060497367449222283392198706D+01 w(1) = 1.1369083326745253D+00 w(2) = 0.93558055763118075D+00 w(3) = 0.87640133443623058D+00 w(4) = 0.87640133443623058D+00 w(5) = 0.93558055763118075D+00 w(6) = 1.1369083326745253D+00 else if ( n .eq. 7 ) then x(1) = - 0.265196135683523349244708200652D+01 x(2) = - 0.167355162876747144503180139830D+01 x(3) = - 0.816287882858964663038710959027D+00 x(4) = 0.0D+00 x(5) = 0.816287882858964663038710959027D+00 x(6) = 0.167355162876747144503180139830D+01 x(7) = 0.265196135683523349244708200652D+01 w(1) = 1.1013307296103216D+00 w(2) = 0.89718460022518409D+00 w(3) = 0.82868730328363926D+00 w(4) = 0.81026461755680734D+00 w(5) = 0.82868730328363926D+00 w(6) = 0.89718460022518409D+00 w(7) = 1.1013307296103216D+00 else if ( n .eq. 8 ) then x(1) = - 0.293063742025724401922350270524D+01 x(2) = - 0.198165675669584292585463063977D+01 x(3) = - 0.115719371244678019472076577906D+01 x(4) = - 0.381186990207322116854718885584D+00 x(5) = 0.381186990207322116854718885584D+00 x(6) = 0.115719371244678019472076577906D+01 x(7) = 0.198165675669584292585463063977D+01 x(8) = 0.293063742025724401922350270524D+01 w(1) = 1.0719301442479805D+00 w(2) = 0.86675260656338138D+00 w(3) = 0.79289004838640131D+00 w(4) = 0.76454412865172916D+00 w(5) = 0.76454412865172916D+00 w(6) = 0.79289004838640131D+00 w(7) = 0.86675260656338138D+00 w(8) = 1.0719301442479805D+00 else if ( n .eq. 9 ) then x(1) = - 0.319099320178152760723004779538D+01 x(2) = - 0.226658058453184311180209693284D+01 x(3) = - 0.146855328921666793166701573925D+01 x(4) = - 0.723551018752837573322639864579D+00 x(5) = 0.0D+00 x(6) = 0.723551018752837573322639864579D+00 x(7) = 0.146855328921666793166701573925D+01 x(8) = 0.226658058453184311180209693284D+01 x(9) = 0.319099320178152760723004779538D+01 w(1) = 1.0470035809766838D+00 w(2) = 0.84175270147867043D+00 w(3) = 0.76460812509455023D+00 w(4) = 0.73030245274509220D+00 w(5) = 0.72023521560605097D+00 w(6) = 0.73030245274509220D+00 w(7) = 0.76460812509455023D+00 w(8) = 0.84175270147867043D+00 w(9) = 1.0470035809766838D+00 else if ( n .eq. 10 ) then x(1) = - 0.343615911883773760332672549432D+01 x(2) = - 0.253273167423278979640896079775D+01 x(3) = - 0.175668364929988177345140122011D+01 x(4) = - 0.103661082978951365417749191676D+01 x(5) = - 0.342901327223704608789165025557D+00 x(6) = 0.342901327223704608789165025557D+00 x(7) = 0.103661082978951365417749191676D+01 x(8) = 0.175668364929988177345140122011D+01 x(9) = 0.253273167423278979640896079775D+01 x(10) = 0.343615911883773760332672549432D+01 w(1) = 1.0254516913657352D+00 w(2) = 0.82066612640481640D+00 w(3) = 0.74144193194356511D+00 w(4) = 0.70329632310490608D+00 w(5) = 0.68708185395127341D+00 w(6) = 0.68708185395127341D+00 w(7) = 0.70329632310490608D+00 w(8) = 0.74144193194356511D+00 w(9) = 0.82066612640481640D+00 w(10) = 1.0254516913657352D+00 else write ( *, '(a)' ) '' write ( *, '(a)' ) 'HERMITE_1_SET - Fatal error!' write ( *, '(a)' ) ' Illegal value of N.' stop 1 end if return end subroutine laguerre_set ( n, x, w ) c*********************************************************************72 c cc LAGUERRE_SET sets abscissas and weights for Laguerre quadrature. c c Discussion: c c The abscissas are the zeroes of the Laguerre polynomial L(N)(X). c c The integral: c c integral ( 0 <= x .lt. +oo ) exp ( -x ) * f(x) dx c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c The integral: c c integral ( 0 <= x .lt. +oo ) f(x) dx c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * exp ( x(i) ) * f ( x(i) ) c c Mathematica can numerically estimate the abscissas for the c n-th order polynomial to p digits of precision by the command: c c NSolve [ LaguerreL[n,x] .eq. 0, x, p ] c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 22 April 2010 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Vladimir Krylov, c Approximate Calculation of Integrals, c Dover, 2006, c ISBN: 0486445798, c LC: QA311.K713. c c Arthur Stroud, Don Secrest, c Gaussian Quadrature Formulas, c Prentice Hall, 1966, c LC: QA299.4G3S7. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Daniel Zwillinger, editor, c CRC Standard Mathematical Tables and Formulae, c 30th Edition, c CRC Press, 1996, c ISBN: 0-8493-2479-3, c LC: QA47.M315. c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 1.00000000000000000000000000000D+00 w(1) = 1.00000000000000000000000000000D+00 else if ( n .eq. 2 ) then x(1) = 0.585786437626904951198311275790D+00 x(2) = 3.41421356237309504880168872421D+00 w(1) = 0.85355339059327376220042218105D+00 w(2) = 0.146446609406726237799577818948D+00 else if ( n .eq. 3 ) then x(1) = 0.415774556783479083311533873128D+00 x(2) = 2.29428036027904171982205036136D+00 x(3) = 6.28994508293747919686641576551D+00 w(1) = 0.71109300992917301544959019114D+00 w(2) = 0.27851773356924084880144488846D+00 w(3) = 0.010389256501586135748964920401D+00 else if ( n .eq. 4 ) then x(1) = 0.322547689619392311800361459104D+00 x(2) = 1.74576110115834657568681671252D+00 x(3) = 4.53662029692112798327928538496D+00 x(4) = 9.39507091230113312923353644342D+00 w(1) = 0.60315410434163360163596602382D+00 w(2) = 0.35741869243779968664149201746D+00 w(3) = 0.03888790851500538427243816816D+00 w(4) = 0.0005392947055613274501037905676D+00 else if ( n .eq. 5 ) then x(1) = 0.263560319718140910203061943361D+00 x(2) = 1.41340305910651679221840798019D+00 x(3) = 3.59642577104072208122318658878D+00 x(4) = 7.08581000585883755692212418111D+00 x(5) = 12.6408008442757826594332193066D+00 w(1) = 0.52175561058280865247586092879D+00 w(2) = 0.3986668110831759274541333481D+00 w(3) = 0.0759424496817075953876533114D+00 w(4) = 0.00361175867992204845446126257D+00 w(5) = 0.00002336997238577622789114908455D+00 else if ( n .eq. 6 ) then x(1) = 0.222846604179260689464354826787D+00 x(2) = 1.18893210167262303074315092194D+00 x(3) = 2.99273632605931407769132528451D+00 x(4) = 5.77514356910451050183983036943D+00 x(5) = 9.83746741838258991771554702994D+00 x(6) = 15.9828739806017017825457915674D+00 w(1) = 0.45896467394996359356828487771D+00 w(2) = 0.4170008307721209941133775662D+00 w(3) = 0.1133733820740449757387061851D+00 w(4) = 0.01039919745314907489891330285D+00 w(5) = 0.000261017202814932059479242860D+00 w(6) = 8.98547906429621238825292053D-07 else if ( n .eq. 7 ) then x(1) = 0.193043676560362413838247885004D+00 x(2) = 1.02666489533919195034519944317D+00 x(3) = 2.56787674495074620690778622666D+00 x(4) = 4.90035308452648456810171437810D+00 x(5) = 8.18215344456286079108182755123D+00 x(6) = 12.7341802917978137580126424582D+00 x(7) = 19.3957278622625403117125820576D+00 w(1) = 0.40931895170127390213043288002D+00 w(2) = 0.4218312778617197799292810054D+00 w(3) = 0.1471263486575052783953741846D+00 w(4) = 0.0206335144687169398657056150D+00 w(5) = 0.00107401014328074552213195963D+00 w(6) = 0.0000158654643485642012687326223D+00 w(7) = 3.17031547899558056227132215D-08 else if ( n .eq. 8 ) then x(1) = 0.170279632305100999788861856608D+00 x(2) = 0.903701776799379912186020223555D+00 x(3) = 2.25108662986613068930711836697D+00 x(4) = 4.26670017028765879364942182690D+00 x(5) = 7.04590540239346569727932548212D+00 x(6) = 10.7585160101809952240599567880D+00 x(7) = 15.7406786412780045780287611584D+00 x(8) = 22.8631317368892641057005342974D+00 w(1) = 0.36918858934163752992058283938D+00 w(2) = 0.4187867808143429560769785813D+00 w(3) = 0.175794986637171805699659867D+00 w(4) = 0.033343492261215651522132535D+00 w(5) = 0.0027945362352256725249389241D+00 w(6) = 0.00009076508773358213104238501D+00 w(7) = 8.4857467162725315448680183D-07 w(8) = 1.04800117487151038161508854D-09 else if ( n .eq. 9 ) then x(1) = 0.152322227731808247428107073127D+00 x(2) = 0.807220022742255847741419210952D+00 x(3) = 2.00513515561934712298303324701D+00 x(4) = 3.78347397333123299167540609364D+00 x(5) = 6.20495677787661260697353521006D+00 x(6) = 9.37298525168757620180971073215D+00 x(7) = 13.4662369110920935710978818397D+00 x(8) = 18.8335977889916966141498992996D+00 x(9) = 26.3740718909273767961410072937D+00 w(1) = 0.336126421797962519673467717606D+00 w(2) = 0.411213980423984387309146942793D+00 w(3) = 0.199287525370885580860575607212D+00 w(4) = 0.0474605627656515992621163600479D+00 w(5) = 0.00559962661079458317700419900556D+00 w(6) = 0.000305249767093210566305412824291D+00 w(7) = 6.59212302607535239225572284875D-06 w(8) = 4.1107693303495484429024104033D-08 w(9) = 3.29087403035070757646681380323D-11 else if ( n .eq. 10 ) then x(1) = 0.137793470540492430830772505653D+00 x(2) = 0.729454549503170498160373121676D+00 x(3) = 1.80834290174031604823292007575D+00 x(4) = 3.40143369785489951448253222141D+00 x(5) = 5.55249614006380363241755848687D+00 x(6) = 8.33015274676449670023876719727D+00 x(7) = 11.8437858379000655649185389191D+00 x(8) = 16.2792578313781020995326539358D+00 x(9) = 21.9965858119807619512770901956D+00 x(10) = 29.9206970122738915599087933408D+00 w(1) = 0.30844111576502014154747083468D+00 w(2) = 0.4011199291552735515157803099D+00 w(3) = 0.218068287611809421588648523D+00 w(4) = 0.062087456098677747392902129D+00 w(5) = 0.009501516975181100553839072D+00 w(6) = 0.0007530083885875387754559644D+00 w(7) = 0.00002825923349599565567422564D+00 w(8) = 4.249313984962686372586577D-07 w(9) = 1.839564823979630780921535D-09 w(10) = 9.911827219609008558377547D-13 else write ( *, '(a)' ) '' write ( *, '(a)' ) 'LAGUERRE_SET - Fatal error!' write ( *, '(a)' ) ' Illegal value of N.' stop 1 end if return end subroutine laguerre_1_set ( n, x, w ) c*********************************************************************72 c cc LAGUERRE_1_SET sets abscissas and weights for Laguerre quadrature. c c Discussion: c c This routine is specialized for the case where the density function is 1. c c The integral is: c I(f) = integral ( 0 <= x .lt. +oo ) f(x) dx c The quadrature rule is c Q(f) = sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 16 May 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 1.00000000000000000000000000000D+00 w(1) = 2.7182818284590451D+00 else if ( n .eq. 2 ) then x(1) = 0.585786437626904951198311275790D+00 x(2) = 3.41421356237309504880168872421D+00 w(1) = 1.5333260331194167D+00 w(2) = 4.4509573350545928D+00 else if ( n .eq. 3 ) then x(1) = 0.415774556783479083311533873128D+00 x(2) = 2.29428036027904171982205036136D+00 x(3) = 6.28994508293747919686641576551D+00 w(1) = 1.0776928592709207D+00 w(2) = 2.7621429619015876D+00 w(3) = 5.6010946254344267D+00 else if ( n .eq. 4 ) then x(1) = 0.322547689619392311800361459104D+00 x(2) = 1.74576110115834657568681671252D+00 x(3) = 4.53662029692112798327928538496D+00 x(4) = 9.39507091230113312923353644342D+00 w(1) = 0.83273912383788917D+00 w(2) = 2.0481024384542965D+00 w(3) = 3.6311463058215168D+00 w(4) = 6.4871450844076604D+00 else if ( n .eq. 5 ) then x(1) = 0.263560319718140910203061943361D+00 x(2) = 1.41340305910651679221840798019D+00 x(3) = 3.59642577104072208122318658878D+00 x(4) = 7.08581000585883755692212418111D+00 x(5) = 12.6408008442757826594332193066D+00 w(1) = 0.67909404220775038D+00 w(2) = 1.6384878736027471D+00 w(3) = 2.7694432423708375D+00 w(4) = 4.3156569009208940D+00 w(5) = 7.2191863543544450D+00 else if ( n .eq. 6 ) then x(1) = 0.222846604179260689464354826787D+00 x(2) = 1.18893210167262303074315092194D+00 x(3) = 2.99273632605931407769132528451D+00 x(4) = 5.77514356910451050183983036943D+00 x(5) = 9.83746741838258991771554702994D+00 x(6) = 15.9828739806017017825457915674D+00 w(1) = 0.57353550742273818D+00 w(2) = 1.3692525907123045D+00 w(3) = 2.2606845933826722D+00 w(4) = 3.3505245823554555D+00 w(5) = 4.8868268002108213D+00 w(6) = 7.8490159455958279D+00 else if ( n .eq. 7 ) then x(1) = 0.193043676560362413838247885004D+00 x(2) = 1.02666489533919195034519944317D+00 x(3) = 2.56787674495074620690778622666D+00 x(4) = 4.90035308452648456810171437810D+00 x(5) = 8.18215344456286079108182755123D+00 x(6) = 12.7341802917978137580126424582D+00 x(7) = 19.3957278622625403117125820576D+00 w(1) = 0.49647759753997234D+00 w(2) = 1.1776430608611976D+00 w(3) = 1.9182497816598063D+00 w(4) = 2.7718486362321113D+00 w(5) = 3.8412491224885148D+00 w(6) = 5.3806782079215330D+00 w(7) = 8.4054324868283103D+00 else if ( n .eq. 8 ) then x(1) = 0.170279632305100999788861856608D+00 x(2) = 0.903701776799379912186020223555D+00 x(3) = 2.25108662986613068930711836697D+00 x(4) = 4.26670017028765879364942182690D+00 x(5) = 7.04590540239346569727932548212D+00 x(6) = 10.7585160101809952240599567880D+00 x(7) = 15.7406786412780045780287611584D+00 x(8) = 22.8631317368892641057005342974D+00 w(1) = 0.43772341049291136D+00 w(2) = 1.0338693476655976D+00 w(3) = 1.6697097656587756D+00 w(4) = 2.3769247017585995D+00 w(5) = 3.2085409133479259D+00 w(6) = 4.2685755108251344D+00 w(7) = 5.8180833686719184D+00 w(8) = 8.9062262152922216D+00 else if ( n .eq. 9 ) then x(1) = 0.152322227731808247428107073127D+00 x(2) = 0.807220022742255847741419210952D+00 x(3) = 2.00513515561934712298303324701D+00 x(4) = 3.78347397333123299167540609364D+00 x(5) = 6.20495677787661260697353521006D+00 x(6) = 9.37298525168757620180971073215D+00 x(7) = 13.4662369110920935710978818397D+00 x(8) = 18.8335977889916966141498992996D+00 x(9) = 26.3740718909273767961410072937D+00 w(1) = 0.39143112431563987D+00 w(2) = 0.92180502852896307D+00 w(3) = 1.4801279099429154D+00 w(4) = 2.0867708075492613D+00 w(5) = 2.7729213897119713D+00 w(6) = 3.5916260680922663D+00 w(7) = 4.6487660021402037D+00 w(8) = 6.2122754197471348D+00 w(9) = 9.3632182377057980D+00 else if ( n .eq. 10 ) then x(1) = 0.137793470540492430830772505653D+00 x(2) = 0.729454549503170498160373121676D+00 x(3) = 1.80834290174031604823292007575D+00 x(4) = 3.40143369785489951448253222141D+00 x(5) = 5.55249614006380363241755848687D+00 x(6) = 8.33015274676449670023876719727D+00 x(7) = 11.8437858379000655649185389191D+00 x(8) = 16.2792578313781020995326539358D+00 x(9) = 21.9965858119807619512770901956D+00 x(10) = 29.9206970122738915599087933408D+00 w(1) = 0.35400973860699630D+00 w(2) = 0.83190230104358065D+00 w(3) = 1.3302885617493283D+00 w(4) = 1.8630639031111309D+00 w(5) = 2.4502555580830108D+00 w(6) = 3.1227641551351848D+00 w(7) = 3.9341526955615240D+00 w(8) = 4.9924148721930299D+00 w(9) = 6.5722024851307994D+00 w(10) = 9.7846958403746243D+00 else write ( *, '(a)' ) '' write ( *, '(a)' ) 'LAGUERRE_1_SET - Fatal error!' write ( *, '(a)' ) ' Illegal value of N.' stop 1 end if return end subroutine legendre_set ( n, x, w ) c*********************************************************************72 c cc LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature. c c Discussion: c c The integral: c c integral ( -1 <= x <= 1 ) f(x) dx c c The quadrature rule: c c sum ( 1 <= i <= n ) w(i) * f ( x(i) ) c c The quadrature rule is exact for polynomials through degree 2*N-1. c c The abscissas are the zeroes of the Legendre polynomial P(N)(X). c c Mathematica can compute the abscissas and weights of a Gauss-Legendre c rule of order N for the interval [A,B] with P digits of precision c by the commands: c c Needs["NumericalDifferentialEquationAnalysis`"] c GaussianQuadratureWeights [ n, a, b, p ] c c Licensing: c c This code is distributed under the GNU LGPL license. c c Modified: c c 20 April 2010 c c Author: c c John Burkardt c c Reference: c c Milton Abramowitz, Irene Stegun, c Handbook of Mathematical Functions, c National Bureau of Standards, 1964, c ISBN: 0-486-61272-4, c LC: QA47.A34. c c Vladimir Krylov, c Approximate Calculation of Integrals, c Dover, 2006, c ISBN: 0486445798, c LC: QA311.K713. c c Arthur Stroud, Don Secrest, c Gaussian Quadrature Formulas, c Prentice Hall, 1966, c LC: QA299.4G3S7. c c Stephen Wolfram, c The Mathematica Book, c Fourth Edition, c Cambridge University Press, 1999, c ISBN: 0-521-64314-7, c LC: QA76.95.W65. c c Daniel Zwillinger, editor, c CRC Standard Mathematical Tables and Formulae, c 30th Edition, c CRC Press, 1996, c ISBN: 0-8493-2479-3, c LC: QA47.M315. c c Parameters: c c Input, integer N, the order. c N must be between 1 and 10. c c Output, double precision X(N), the abscissas. c c Output, double precision W(N), the weights. c implicit none integer n double precision w(n) double precision x(n) if ( n .eq. 1 ) then x(1) = 0.000000000000000000000000000000D+00 w(1) = 2.000000000000000000000000000000D+00 else if ( n .eq. 2 ) then x(1) = -0.577350269189625764509148780502D+00 x(2) = 0.577350269189625764509148780502D+00 w(1) = 1.000000000000000000000000000000D+00 w(2) = 1.000000000000000000000000000000D+00 else if ( n .eq. 3 ) then x(1) = -0.774596669241483377035853079956D+00 x(2) = 0.000000000000000000000000000000D+00 x(3) = 0.774596669241483377035853079956D+00 w(1) = 0.555555555555555555555555555556D+00 w(2) = 0.888888888888888888888888888889D+00 w(3) = 0.555555555555555555555555555556D+00 else if ( n .eq. 4 ) then x(1) = -0.861136311594052575223946488893D+00 x(2) = -0.339981043584856264802665759103D+00 x(3) = 0.339981043584856264802665759103D+00 x(4) = 0.861136311594052575223946488893D+00 w(1) = 0.347854845137453857373063949222D+00 w(2) = 0.652145154862546142626936050778D+00 w(3) = 0.652145154862546142626936050778D+00 w(4) = 0.347854845137453857373063949222D+00 else if ( n .eq. 5 ) then x(1) = -0.906179845938663992797626878299D+00 x(2) = -0.538469310105683091036314420700D+00 x(3) = 0.000000000000000000000000000000D+00 x(4) = 0.538469310105683091036314420700D+00 x(5) = 0.906179845938663992797626878299D+00 w(1) = 0.236926885056189087514264040720D+00 w(2) = 0.478628670499366468041291514836D+00 w(3) = 0.568888888888888888888888888889D+00 w(4) = 0.478628670499366468041291514836D+00 w(5) = 0.236926885056189087514264040720D+00 else if ( n .eq. 6 ) then x(1) = -0.932469514203152027812301554494D+00 x(2) = -0.661209386466264513661399595020D+00 x(3) = -0.238619186083196908630501721681D+00 x(4) = 0.238619186083196908630501721681D+00 x(5) = 0.661209386466264513661399595020D+00 x(6) = 0.932469514203152027812301554494D+00 w(1) = 0.171324492379170345040296142173D+00 w(2) = 0.360761573048138607569833513838D+00 w(3) = 0.467913934572691047389870343990D+00 w(4) = 0.467913934572691047389870343990D+00 w(5) = 0.360761573048138607569833513838D+00 w(6) = 0.171324492379170345040296142173D+00 else if ( n .eq. 7 ) then x(1) = -0.949107912342758524526189684048D+00 x(2) = -0.741531185599394439863864773281D+00 x(3) = -0.405845151377397166906606412077D+00 x(4) = 0.000000000000000000000000000000D+00 x(5) = 0.405845151377397166906606412077D+00 x(6) = 0.741531185599394439863864773281D+00 x(7) = 0.949107912342758524526189684048D+00 w(1) = 0.129484966168869693270611432679D+00 w(2) = 0.279705391489276667901467771424D+00 w(3) = 0.381830050505118944950369775489D+00 w(4) = 0.417959183673469387755102040816D+00 w(5) = 0.381830050505118944950369775489D+00 w(6) = 0.279705391489276667901467771424D+00 w(7) = 0.129484966168869693270611432679D+00 else if ( n .eq. 8 ) then x(1) = -0.960289856497536231683560868569D+00 x(2) = -0.796666477413626739591553936476D+00 x(3) = -0.525532409916328985817739049189D+00 x(4) = -0.183434642495649804939476142360D+00 x(5) = 0.183434642495649804939476142360D+00 x(6) = 0.525532409916328985817739049189D+00 x(7) = 0.796666477413626739591553936476D+00 x(8) = 0.960289856497536231683560868569D+00 w(1) = 0.101228536290376259152531354310D+00 w(2) = 0.222381034453374470544355994426D+00 w(3) = 0.313706645877887287337962201987D+00 w(4) = 0.362683783378361982965150449277D+00 w(5) = 0.362683783378361982965150449277D+00 w(6) = 0.313706645877887287337962201987D+00 w(7) = 0.222381034453374470544355994426D+00 w(8) = 0.101228536290376259152531354310D+00 else if ( n .eq. 9 ) then x(1) = -0.968160239507626089835576203D+00 x(2) = -0.836031107326635794299429788D+00 x(3) = -0.613371432700590397308702039D+00 x(4) = -0.324253423403808929038538015D+00 x(5) = 0.000000000000000000000000000D+00 x(6) = 0.324253423403808929038538015D+00 x(7) = 0.613371432700590397308702039D+00 x(8) = 0.836031107326635794299429788D+00 x(9) = 0.968160239507626089835576203D+00 w(1) = 0.081274388361574411971892158111D+00 w(2) = 0.18064816069485740405847203124D+00 w(3) = 0.26061069640293546231874286942D+00 w(4) = 0.31234707704000284006863040658D+00 w(5) = 0.33023935500125976316452506929D+00 w(6) = 0.31234707704000284006863040658D+00 w(7) = 0.26061069640293546231874286942D+00 w(8) = 0.18064816069485740405847203124D+00 w(9) = 0.081274388361574411971892158111D+00 else if ( n .eq. 10 ) then x(1) = -0.973906528517171720077964012D+00 x(2) = -0.865063366688984510732096688D+00 x(3) = -0.679409568299024406234327365D+00 x(4) = -0.433395394129247190799265943D+00 x(5) = -0.148874338981631210884826001D+00 x(6) = 0.148874338981631210884826001D+00 x(7) = 0.433395394129247190799265943D+00 x(8) = 0.679409568299024406234327365D+00 x(9) = 0.865063366688984510732096688D+00 x(10) = 0.973906528517171720077964012D+00 w(1) = 0.066671344308688137593568809893D+00 w(2) = 0.14945134915058059314577633966D+00 w(3) = 0.21908636251598204399553493423D+00 w(4) = 0.26926671930999635509122692157D+00 w(5) = 0.29552422471475287017389299465D+00 w(6) = 0.29552422471475287017389299465D+00 w(7) = 0.26926671930999635509122692157D+00 w(8) = 0.21908636251598204399553493423D+00 w(9) = 0.14945134915058059314577633966D+00 w(10) = 0.066671344308688137593568809893D+00 else write ( *, '(a)' ) '' write ( *, '(a)' ) 'LEGENDRE_SET - Fatal error!' write ( *, '(a)' ) ' Illegal value of N.' stop 1 end if return end