#! /usr/bin/env python # def kronrod_set ( n ): #*****************************************************************************80 # ## KRONROD_SET sets abscissas and weights for Gauss-Kronrod quadrature. # # Discussion: # # The integral: # # Integral ( -1 <= X <= 1 ) F(X) dX # # The quadrature rule: # # Sum ( 1 <= I <= N ) W(I) * F ( X(I) ) # # A Kronrod rule is used in conjunction with a lower order # Gauss rule, and provides an efficient error estimation. # # The error may be estimated as the difference in the two integral # approximations. # # The efficiency comes about because the Kronrod uses the abscissas # of the Gauss rule, thus saving on the number of function evaluations # necessary. If the Kronrod rule were replaced by a Gauss rule of # the same order, a higher precision integral estimate would be # made, but the function would have to be evaluated at many more # points. # # The Gauss Kronrod pair of rules involves an ( N + 1 ) / 2 # point Gauss-Legendre rule and an N point Kronrod rule. # Thus, the 15 point Kronrod rule should be paired with the # Gauss-Legendre 7 point rule. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 18 June 2015 # # Author: # # John Burkardt # # Reference: # # R Piessens, E de Doncker-Kapenger, C Ueberhuber, D Kahaner, # QUADPACK, A Subroutine Package for Automatic Integration, # Springer Verlag, 1983. # # Parameters: # # Input, integer N, the order. # N may be 15, 21, 31 or 41, corresponding to Gauss-Legendre rules of # order 7, 10, 15 or 20. # # Output, real X(N), the abscissas. # # Output, real W(N), the weights. # import numpy as np from sys import exit if ( n == 15 ): x = np.array ( [ \ - 0.9914553711208126, \ - 0.9491079123427585, \ - 0.8648644233597691, \ - 0.7415311855993944, \ - 0.5860872354676911, \ - 0.4058451513773972, \ - 0.2077849550789850, \ 0.0, \ 0.2077849550789850, \ 0.4058451513773972, \ 0.5860872354676911, \ 0.7415311855993944, \ 0.8648644233597691, \ 0.9491079123427585, \ 0.9914553711208126 ] ) w = np.array ( [ \ 0.2293532201052922E-01, \ 0.6309209262997855E-01, \ 0.1047900103222502, \ 0.1406532597155259, \ 0.1690047266392679, \ 0.1903505780647854, \ 0.2044329400752989, \ 0.2094821410847278, \ 0.2044329400752989, \ 0.1903505780647854, \ 0.1690047266392679, \ 0.1406532597155259, \ 0.1047900103222502, \ 0.6309209262997855E-01, \ 0.2293532201052922E-01 ] ) elif ( n == 21 ): x = np.array ( [ \ - 0.9956571630258081, \ - 0.9739065285171717, \ - 0.9301574913557082, \ - 0.8650633666889845, \ - 0.7808177265864169, \ - 0.6794095682990244, \ - 0.5627571346686047, \ - 0.4333953941292472, \ - 0.2943928627014602, \ - 0.1488743389816312, \ 0.0, \ 0.1488743389816312, \ 0.2943928627014602, \ 0.4333953941292472, \ 0.5627571346686047, \ 0.6794095682990244, \ 0.7808177265864169, \ 0.8650633666889845, \ 0.9301574913557082, \ 0.9739065285171717, \ 0.9956571630258081 ] ) w = np.array ( [ \ 0.1169463886737187E-01, \ 0.3255816230796473E-01, \ 0.5475589657435200E-01, \ 0.7503967481091995E-01, \ 0.9312545458369761E-01, \ 0.1093871588022976, \ 0.1234919762620659, \ 0.1347092173114733, \ 0.1427759385770601, \ 0.1477391049013385, \ 0.1494455540029169, \ 0.1477391049013385, \ 0.1427759385770601, \ 0.1347092173114733, \ 0.1234919762620659, \ 0.1093871588022976, \ 0.9312545458369761E-01, \ 0.7503967481091995E-01, \ 0.5475589657435200E-01, \ 0.3255816230796473E-01, \ 0.1169463886737187E-01 ] ) elif ( n == 31 ): x = np.array ( [ \ - 0.9980022986933971, \ - 0.9879925180204854, \ - 0.9677390756791391, \ - 0.9372733924007059, \ - 0.8972645323440819, \ - 0.8482065834104272, \ - 0.7904185014424659, \ - 0.7244177313601700, \ - 0.6509967412974170, \ - 0.5709721726085388, \ - 0.4850818636402397, \ - 0.3941513470775634, \ - 0.2991800071531688, \ - 0.2011940939974345, \ - 0.1011420669187175, \ 0.0, \ 0.1011420669187175, \ 0.2011940939974345, \ 0.2991800071531688, \ 0.3941513470775634, \ 0.4850818636402397, \ 0.5709721726085388, \ 0.6509967412974170, \ 0.7244177313601700, \ 0.7904185014424659, \ 0.8482065834104272, \ 0.8972645323440819, \ 0.9372733924007059, \ 0.9677390756791391, \ 0.9879925180204854, \ 0.9980022986933971 ] ) w = np.array ( [ \ 0.5377479872923349E-02, \ 0.1500794732931612E-01, \ 0.2546084732671532E-01, \ 0.3534636079137585E-01, \ 0.4458975132476488E-01, \ 0.5348152469092809E-01, \ 0.6200956780067064E-01, \ 0.6985412131872826E-01, \ 0.7684968075772038E-01, \ 0.8308050282313302E-01, \ 0.8856444305621177E-01, \ 0.9312659817082532E-01, \ 0.9664272698362368E-01, \ 0.9917359872179196E-01, \ 0.1007698455238756, \ 0.1013300070147915, \ 0.1007698455238756, \ 0.9917359872179196E-01, \ 0.9664272698362368E-01, \ 0.9312659817082532E-01, \ 0.8856444305621177E-01, \ 0.8308050282313302E-01, \ 0.7684968075772038E-01, \ 0.6985412131872826E-01, \ 0.6200956780067064E-01, \ 0.5348152469092809E-01, \ 0.4458975132476488E-01, \ 0.3534636079137585E-01, \ 0.2546084732671532E-01, \ 0.1500794732931612E-01, \ 0.5377479872923349E-02 ] ) elif ( n == 41 ): x = np.array ( [ \ - 0.9988590315882777, \ - 0.9931285991850949, \ - 0.9815078774502503, \ - 0.9639719272779138, \ - 0.9408226338317548, \ - 0.9122344282513259, \ - 0.8782768112522820, \ - 0.8391169718222188, \ - 0.7950414288375512, \ - 0.7463319064601508, \ - 0.6932376563347514, \ - 0.6360536807265150, \ - 0.5751404468197103, \ - 0.5108670019508271, \ - 0.4435931752387251, \ - 0.3737060887154196, \ - 0.3016278681149130, \ - 0.2277858511416451, \ - 0.1526054652409227, \ - 0.7652652113349733E-01, \ 0.0, \ 0.7652652113349733E-01, \ 0.1526054652409227, \ 0.2277858511416451, \ 0.3016278681149130, \ 0.3737060887154196, \ 0.4435931752387251, \ 0.5108670019508271, \ 0.5751404468197103, \ 0.6360536807265150, \ 0.6932376563347514, \ 0.7463319064601508, \ 0.7950414288375512, \ 0.8391169718222188, \ 0.8782768112522820, \ 0.9122344282513259, \ 0.9408226338317548, \ 0.9639719272779138, \ 0.9815078774502503, \ 0.9931285991850949, \ 0.9988590315882777 ] ) w = np.array ( [ \ 0.3073583718520532E-02, \ 0.8600269855642942E-02, \ 0.1462616925697125E-01, \ 0.2038837346126652E-01, \ 0.2588213360495116E-01, \ 0.3128730677703280E-01, \ 0.3660016975820080E-01, \ 0.4166887332797369E-01, \ 0.4643482186749767E-01, \ 0.5094457392372869E-01, \ 0.5519510534828599E-01, \ 0.5911140088063957E-01, \ 0.6265323755478117E-01, \ 0.6583459713361842E-01, \ 0.6864867292852162E-01, \ 0.7105442355344407E-01, \ 0.7303069033278667E-01, \ 0.7458287540049919E-01, \ 0.7570449768455667E-01, \ 0.7637786767208074E-01, \ 0.7660071191799966E-01, \ 0.7637786767208074E-01, \ 0.7570449768455667E-01, \ 0.7458287540049919E-01, \ 0.7303069033278667E-01, \ 0.7105442355344407E-01, \ 0.6864867292852162E-01, \ 0.6583459713361842E-01, \ 0.6265323755478117E-01, \ 0.5911140088063957E-01, \ 0.5519510534828599E-01, \ 0.5094457392372869E-01, \ 0.4643482186749767E-01, \ 0.4166887332797369E-01, \ 0.3660016975820080E-01, \ 0.3128730677703280E-01, \ 0.2588213360495116E-01, \ 0.2038837346126652E-01, \ 0.1462616925697125E-01, \ 0.8600269855642942E-02, \ 0.3073583718520532E-02 ] ) else: print ( '' ) print ( 'KRONROD_SET - Fatal error!' ) print ( ' Illegal value of N = %d' % ( n ) ) print ( ' Legal values are 15, 21, 31 or 41.' ) exit ( 'KRONROD_SET - Fatal error!' ) return x, w def kronrod_set_test ( ): #*****************************************************************************80 # ## KRONROD_SET_TEST tests KRONROD_SET. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 18 June 2015 # # Author: # # John Burkardt # import numpy as np import platform from legendre_set import legendre_set nl_test = np.array ( [ 7, 10, 15, 20 ] ) print ( '' ) print ( 'KRONROD_SET_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' KRONROD_SET sets up a Kronrod quadrature rule,' ) print ( ' This is used following a lower order Legendre rule.' ) for test in range ( 0, 4 ): print ( '' ) print ( ' Legendre/Kronrod quadrature pair #%d' % ( test ) ) print ( ' W X' ) print ( '' ) nl = nl_test[test] xl, wl = legendre_set ( nl ) for i in range ( 0, nl ): print ( ' %8d %24.16g %24.16g' % ( i, xl[i], wl[i] ) ) print ( '' ) nk = 2 * nl + 1 xk, wk = kronrod_set ( nk ) for i in range ( 0, nk ): print ( ' %8d %24.16g %24.16g' % ( i, xk[i], wk[i] ) ) # # Terminate. # print ( '' ) print ( 'KRONROD_SET_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) kronrod_set_test ( ) timestamp ( )