#! /usr/bin/env python # def radau_set ( n ): #*****************************************************************************80 # ## RADAU_SET sets abscissas and weights for Radau quadrature. # # Discussion: # # The Radau rule is distinguished by the fact that the left endpoint # (-1) is always an abscissa. # # The integral: # # Integral ( -1 <= X <= 1 ) F(X) dX # # The quadrature rule: # # Sum ( 1 <= I <= N ) W(I) * F ( X(I) ) # # The quadrature rule will integrate exactly all polynomials up to # X^(2*N-2). # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 12 October 2005 # # Author: # # John Burkardt # # Reference: # # Milton Abramowitz, Irene Stegun, # Handbook of Mathematical Functions, # National Bureau of Standards, 1964, # ISBN: 0-486-61272-4, # LC: QA47.A34. # # Arthur Stroud, Don Secrest, # Gaussian Quadrature Formulas, # Prentice Hall, 1966, # LC: QA299.4G3S7. # # Daniel Zwillinger, editor, # CRC Standard Mathematical Tables and Formulae, # 30th Edition, # CRC Press, 1996, # ISBN: 0-8493-2479-3. # # Parameters: # # Input, integer N, the order. # N must be between 1 and 15. # # Output, real X(N), the abscissas. # # Output, real W(N), the weights. # import numpy as np from sys import exit if ( n == 1 ): x = np.array ( [ \ - 1.0 ] ) w = np.array ( [ \ 2.0 ] ) elif ( n == 2 ): x = np.array ( [ \ - 1.0, \ 1.0 / 3.0 ] ) w = np.array ( [ \ 0.5, \ 1.5 ] ) elif ( n == 3 ): x = np.array ( [ \ - 1.0, \ - 0.289897948556635619639456814941, \ 0.689897948556635619639456814941 ] ) w = np.array ( [ \ 0.222222222222222222222222222222, \ 1.02497165237684322767762689304, \ 0.752806125400934550100150884739 ] ) elif ( n == 4 ): x = np.array ( [ \ - 1.0, \ - 0.575318923521694112050483779752, \ 0.181066271118530578270147495862, \ 0.822824080974592105208907712461 ] ) w = np.array ( [ \ 0.125, \ 0.657688639960119487888578442146, \ 0.776386937686343761560464613780, \ 0.440924422353536750550956944074 ] ) elif ( n == 5 ): x = np.array ( [ \ - 1.0, \ - 0.720480271312438895695825837750, \ - 0.167180864737833640113395337326, \ 0.446313972723752344639908004629, \ 0.885791607770964635613757614892 ] ) w = np.array ( [ \ 0.08, \ 0.446207802167141488805120436457, \ 0.623653045951482508163709823153, \ 0.562712030298924120384345300681, \ 0.287427121582451882646824439708 ] ) elif ( n == 6 ): x = np.array ( [ \ - 1.0, \ - 0.802929828402347147753002204224, \ - 0.390928546707272189029229647442, \ 0.124050379505227711989974959990, \ 0.603973164252783654928415726409, \ 0.920380285897062515318386619813 ] ) w = np.array ( [ \ 0.0555555555555555555555555555556, \ 0.319640753220510966545779983796, \ 0.485387188468969916159827915587, \ 0.520926783189574982570229406570, \ 0.416901334311907738959406382743, \ 0.201588385253480840209200755749 ] ) elif ( n == 7 ): x = np.array ( [ \ - 1.0, \ - 0.853891342639482229703747931639, \ - 0.538467724060109001833766720231, \ - 0.117343037543100264162786683611, \ 0.326030619437691401805894055838, \ 0.703842800663031416300046295008, \ 0.941367145680430216055899446174 ] ) w = np.array ( [ \ 0.0408163265306122448979591836735, \ 0.239227489225312405787077480770, \ 0.380949873644231153805938347876, \ 0.447109829014566469499348953642, \ 0.424703779005955608398308039150, \ 0.318204231467301481744870434470, \ 0.148988471112020635866497560418 ] ) elif ( n == 8 ): x = np.array ( [ \ - 1.0, \ - 0.887474878926155707068695617935, \ - 0.639518616526215270024840114382, \ - 0.294750565773660725252184459658, \ 0.0943072526611107660028971153047, \ 0.468420354430821063046421216613, \ 0.770641893678191536180719525865, \ 0.955041227122575003782349000858 ] ) w = np.array ( [ \ 0.03125, \ 0.185358154802979278540728972699, \ 0.304130620646785128975743291400, \ 0.376517545389118556572129261442, \ 0.391572167452493593082499534004, \ 0.347014795634501280228675918422, \ 0.249647901329864963257869293513, \ 0.114508814744257199342353728520 ] ) elif ( n == 9 ): x = np.array ( [ \ - 1.0, \ - 0.910732089420060298533757956283, \ - 0.711267485915708857029562959544, \ - 0.426350485711138962102627520502, \ - 0.0903733696068532980645444599064, \ 0.256135670833455395138292079035, \ 0.571383041208738483284917464837, \ 0.817352784200412087992517083851, \ 0.964440169705273096373589797925 ] ) w = np.array ( [ \ 0.0246913580246913580246913580247, \ 0.147654019046315385819588499802, \ 0.247189378204593052361239794969, \ 0.316843775670437978338000849642, \ 0.348273002772966594071991031186, \ 0.337693966975929585803724239792, \ 0.286386696357231171146705637752, \ 0.200553298024551957421165090417, \ 0.907145049232829170128934984159E-01 ] ) elif ( n == 10 ): x = np.array ( [ \ - 1.0, \ - 0.927484374233581078117671398464, \ - 0.763842042420002599615429776011, \ - 0.525646030370079229365386614293, \ - 0.236234469390588049278459503207, \ 0.760591978379781302337137826389E-01, \ 0.380664840144724365880759065541, \ 0.647766687674009436273648507855, \ 0.851225220581607910728163628088, \ 0.971175180702246902734346518378 ] ) w = np.array ( [ \ 0.02, \ 0.120296670557481631517310522702, \ 0.204270131879000675555788672223, \ 0.268194837841178696058554475262, \ 0.305859287724422621016275475401, \ 0.313582457226938376695902847302, \ 0.290610164832918311146863077963, \ 0.239193431714379713376571966160, \ 0.164376012736921475701681668908, \ 0.736170054867584989310512940790E-01 ] ) elif ( n == 11 ): x = np.array ( [ \ - 1.0, \ - 0.939941935677027005913871284731, \ - 0.803421975580293540697597956820, \ - 0.601957842073797690275892603234, \ - 0.351888923353330214714301017870, \ - 0.734775314313212657461903554238E-01, \ 0.210720306228426314076095789845, \ 0.477680647983087519467896683890, \ 0.705777100713859519144801128840, \ 0.876535856245703748954741265611, \ 0.976164773135168806180508826082 ] ) w = np.array ( [ \ 0.165289256198347107438016528926E-01, \ 0.998460819079680638957534695802E-01, \ 0.171317619206659836486712649042, \ 0.228866123848976624401683231126, \ 0.267867086189684177806638163355, \ 0.285165563941007337460004408915, \ 0.279361333103383045188962195720, \ 0.250925377697128394649140267633, \ 0.202163108540024418349931754266, \ 0.137033682133202256310153880580, \ 0.609250978121311347072183268883E-01 ] ) elif ( n == 12 ): x = np.array ( [ \ - 1.0, \ - 0.949452759204959300493337627077, \ - 0.833916773105189706586269254036, \ - 0.661649799245637148061133087811, \ - 0.444406569781935851126642615609, \ - 0.196994559534278366455441427346, \ 0.637247738208319158337792384845E-01, \ 0.319983684170669623532789532206, \ 0.554318785912324288984337093085, \ 0.750761549711113852529400825472, \ 0.895929097745638894832914608454, \ 0.979963439076639188313950540264 ] ) w = np.array ( [ \ 0.138888888888888888888888888888E-01, \ 0.841721349386809762415796536813E-01, \ 0.145563668853995128522547654706, \ 0.196998534826089634656049637969, \ 0.235003115144985839348633985940, \ 0.256991338152707776127974253598, \ 0.261465660552133103438074715743, \ 0.248121560804009959403073107079, \ 0.217868879026192438848747482023, \ 0.172770639313308564306065766966, \ 0.115907480291738392750341908272, \ 0.512480992072692974680229451351E-01 ] ) elif ( n == 13 ): x = np.array ( [ \ - 1.0, \ - 0.956875873668299278183813833834, \ - 0.857884202528822035697620310269, \ - 0.709105087529871761580423832811, \ - 0.519197779050454107485205148087, \ - 0.299201300554509985532583446686, \ - 0.619016986256353412578604857936E-01, \ 0.178909837597084635021931298881, \ 0.409238231474839556754166331248, \ 0.615697890940291918017885487543, \ 0.786291018233046684731786459135, \ 0.911107073689184553949066402429, \ 0.982921890023145161262671078244 ] ) w = np.array ( [ \ 0.118343195266272189349112426036E-01, \ 0.719024162924955289397537405641E-01, \ 0.125103834331152358133769287976, \ 0.171003460470616642463758674512, \ 0.206960611455877074631132560829, \ 0.230888862886995434012203758668, \ 0.241398342287691148630866924129, \ 0.237878547660712031342685189180, \ 0.220534229288451464691077164199, \ 0.190373715559631732254759820746, \ 0.149150950090000205151491864242, \ 0.992678068818470859847363877478E-01, \ 0.437029032679020748288533846051E-01 ] ) elif ( n == 14 ): x = np.array ( [ \ - 1.0, \ - 0.962779269978024297120561244319, \ - 0.877048918201462024795266773531, \ - 0.747389642613378838735429134263, \ - 0.580314056546874971105726664999, \ - 0.384202003439203313794083903375, \ - 0.168887928042680911008441695622, \ 0.548312279917645496498107146428E-01, \ 0.275737205435522399182637403545, \ 0.482752918588474966820418534355, \ 0.665497977216884537008955042481, \ 0.814809550601994729434217249123, \ 0.923203722520643299246334950272, \ 0.985270697947821356698617003172 ] ) w = np.array ( [ \ 0.102040816326530612244897959184E-01, \ 0.621220169077714601661329164668E-01, \ 0.108607722744362826826720935229, \ 0.149620539353121355950520836946, \ 0.183127002125729654123867302103, \ 0.207449763335175672668082886489, \ 0.221369811499570948931671683021, \ 0.224189348002707794238414632220, \ 0.215767100604618851381187446115, \ 0.196525518452982430324613091930, \ 0.167429727891086278990102277038, \ 0.129939668737342347807425737146, \ 0.859405354429804030893077310866E-01, \ 0.377071632698969142774627282919E-01 ] ) elif ( n == 15 ): x = np.array ( [ \ - 1.0, \ - 0.967550468197200476562456018282, \ - 0.892605400120550767066811886849, \ - 0.778685617639031079381743321893, \ - 0.630779478886949283946148437224, \ - 0.455352905778529370872053455981, \ - 0.260073376740807915768961188263, \ - 0.534757226797460641074538896258E-01, \ 0.155410685384859484319182024964, \ 0.357456512022127651195319205174, \ 0.543831458701484016930711802760, \ 0.706390264637572540152679669478, \ 0.838029000636089631215097384520, \ 0.932997190935973719928072142859, \ 0.987166478414363086378359071811 ] ) w = np.array ( [ \ 0.00888888888888888888888888888889, \ 0.0542027800486444943382142368018, \ 0.0951295994604808992038477266346, \ 0.131875462504951632186262157944, \ 0.162854477303832629448732245828, \ 0.186715145839450908083795103799, \ 0.202415187030618429872703310435, \ 0.209268608147694581430889790306, \ 0.206975960249553755479027321787, \ 0.195637503045116116473556617575, \ 0.175748872642447685670310440476, \ 0.148179527003467253924682058743, \ 0.114135203489752753013075582569, \ 0.0751083927605064397329716653914, \ 0.0328643915845935322530428528231 ] ) else: print ( '' ) print ( 'RADAU_SET - Fatal error!' ) print ( ' Illegal value of N = %d' % ( n ) ) print ( ' Legal values are 1 to 15.' ) exit ( 'RADAU_SET - Fatal error!' ) return x, w def radau_set_test ( ): #*****************************************************************************80 # ## RADAU_SET_TEST tests RADAU_SET. # # Licensing: # # This code is distributed under the GNU LGPL license. # # Modified: # # 10 June 2015 # # Author: # # John Burkardt # import platform print ( '' ) print ( 'RADAU_SET_TEST' ) print ( ' Python version: %s' % ( platform.python_version ( ) ) ) print ( ' RADAU_SET sets a Radau rule.' ) print ( '' ) print ( ' I X W' ) for n in range ( 4, 12, 3 ): x, w = radau_set ( n ) print ( '' ) for i in range ( 0, n ): print ( ' %8d %24.16g %24.16g' % ( i, x[i], w[i] ) ) # # Terminate. # print ( '' ) print ( 'RADAU_SET_TEST:' ) print ( ' Normal end of execution.' ) return if ( __name__ == '__main__' ): from timestamp import timestamp timestamp ( ) radau_set_test ( ) timestamp ( )