subroutine ncc_set ( n, x, w ) c*********************************************************************72 c cc ncc_set() sets abscissas and weights for Newton-Cotes closed quadrature. c c Discussion: c c The Newton-Cotes closed rules use equally spaced abscissas, and c hence may be used with tabulated function data. c c The rules are called "closed" because they include the endpoints. c As a favor, we include an order 1 rule, the midpoint rule, even c though this does not satisfy the requirement that the endpoints c be included! c c The higher order rules involve negative weights. These can produce c loss of accuracy due to the subtraction of large, nearly equal quantities. c c N = 1 is the midpoint rule (and is not really an NCC rule!) c N = 2 is the trapezoidal rule. c N = 3 is Simpson's rule. c N = 4 is Simpson's 3/8 rule. c N = 5 is Bode's rule. c c The Kopal reference for N = 12 lists c W(6) = 15494566.0D+00 / 43545600.0D+00 c but this results in a set of coeffients that don't add up to 2. c The correct value is c W(6) = 15493566.0D+00 / 43545600.0. 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 In Mathematica, the Newton-Cotes closed weights of order N in c the interval [-1,+1] can be computed by: c c Needs["NumericalDifferentialEquationAnalysis`"] c NewtonCotesWeights [ n, -1, 1, QuadratureType -> Closed ] c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 21 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 Johnson, c Quarterly Journal of Mathematics, c Volume 46, Number 52, 1915. c c Zdenek Kopal, c Numerical Analysis, c John Wiley, 1955, c LC: QA297.K6. 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 1 <= N <= 21. 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 c c 2 c x(1) = 0.00000000000000000000D+00 w(1) = 2.00000000000000000000D+00 else if ( n .eq. 2 ) then c c 1 c 1 c 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 c c 1 / 3 c 4 / 3 c 1 / 3 c 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 c c 1 / 4 c 3 / 4 c 3 / 4 c 1 / 4 c x(1) = -1.00000000000000000000D+00 x(2) = -0.33333333333333333333D+00 x(3) = 0.33333333333333333333D+00 x(4) = 1.00000000000000000000D+00 w(1) = 0.25000000000000000000D+00 w(2) = 0.75000000000000000000D+00 w(3) = 0.75000000000000000000D+00 w(4) = 0.25000000000000000000D+00 else if ( n .eq. 5 ) then c c 7 / 45 c 32 / 45 c 12 / 45 c 32 / 45 c 7 / 45 c x(1) = -1.00000000000000000000D+00 x(2) = -0.50000000000000000000D+00 x(3) = 0.00000000000000000000D+00 x(4) = 0.50000000000000000000D+00 x(5) = 1.00000000000000000000D+00 w(1) = 0.15555555555555555556D+00 w(2) = 0.71111111111111111111D+00 w(3) = 0.26666666666666666667D+00 w(4) = 0.71111111111111111111D+00 w(5) = 0.15555555555555555556D+00 else if ( n .eq. 6 ) then c c 19 / 144 c 75 / 144 c 50 / 144 c 50 / 144 c 75 / 144 c 19 / 144 c x(1) = -1.00000000000000000000D+00 x(2) = -0.60000000000000000000D+00 x(3) = -0.20000000000000000000D+00 x(4) = 0.20000000000000000000D+00 x(5) = 0.60000000000000000000D+00 x(6) = 1.00000000000000000000D+00 w(1) = 0.13194444444444444444D+00 w(2) = 0.52083333333333333333D+00 w(3) = 0.34722222222222222222D+00 w(4) = 0.34722222222222222222D+00 w(5) = 0.52083333333333333333D+00 w(6) = 0.13194444444444444444D+00 else if ( n .eq. 7 ) then c c 41 / 420 c 216 / 420 c 27 / 420 c 272 / 420 c 27 / 420 c 216 / 420 c 41 / 420 c x(1) = -1.00000000000000000000D+00 x(2) = -0.66666666666666666667D+00 x(3) = -0.33333333333333333333D+00 x(4) = 0.00000000000000000000D+00 x(5) = 0.33333333333333333333D+00 x(6) = 0.66666666666666666667D+00 x(7) = 1.00000000000000000000D+00 w(1) = 0.097619047619047619048D+00 w(2) = 0.51428571428571428571D+00 w(3) = 0.064285714285714285714D+00 w(4) = 0.64761904761904761905D+00 w(5) = 0.064285714285714285714D+00 w(6) = 0.51428571428571428571D+00 w(7) = 0.097619047619047619048D+00 else if ( n .eq. 8 ) then c c 751 / 8640 c 3577 / 8640 c 1323 / 8640 c 2989 / 8640 c 2989 / 8640 c 1323 / 8640 c 3577 / 8640 c 751 / 8640 c x(1) = -1.00000000000000000000D+00 x(2) = -0.71428571428571428571D+00 x(3) = -0.42857142857142857143D+00 x(4) = -0.14285714285714285714D+00 x(5) = 0.14285714285714285714D+00 x(6) = 0.42857142857142857143D+00 x(7) = 0.71428571428571428571D+00 x(8) = 1.00000000000000000000D+00 w(1) = 0.086921296296296296296D+00 w(2) = 0.41400462962962962963D+00 w(3) = 0.15312500000000000000D+00 w(4) = 0.34594907407407407407D+00 w(5) = 0.34594907407407407407D+00 w(6) = 0.15312500000000000000D+00 w(7) = 0.41400462962962962963D+00 w(8) = 0.086921296296296296296D+00 else if ( n .eq. 9 ) then c c 989 / 14175 c 5888 / 14175 c -928 / 14175 c 10496 / 14175 c -4540 / 14175 c 10496 / 14175 c -928 / 14175 c 5888 / 14175 c 989 / 14175 c x(1) = -1.00000000000000000000D+00 x(2) = -0.75000000000000000000D+00 x(3) = -0.50000000000000000000D+00 x(4) = -0.25000000000000000000D+00 x(5) = 0.00000000000000000000D+00 x(6) = 0.25000000000000000000D+00 x(7) = 0.50000000000000000000D+00 x(8) = 0.75000000000000000000D+00 x(9) = 1.00000000000000000000D+00 w(1) = 0.069770723104056437390D+00 w(2) = 0.41537918871252204586D+00 w(3) = -0.065467372134038800705D+00 w(4) = 0.74045855379188712522D+00 w(5) = -0.32028218694885361552D+00 w(6) = 0.74045855379188712522D+00 w(7) = -0.065467372134038800705D+00 w(8) = 0.41537918871252204586D+00 w(9) = 0.069770723104056437390D+00 else if ( n .eq. 10 ) then c c 2857 / 44800 c 15741 / 44800 c 1080 / 44800 c 19344 / 44800 c 5778 / 44800 c 5778 / 44800 c 19344 / 44800 c 1080 / 44800 c 15741 / 44800 c 2857 / 44800 c x(1) = -1.00000000000000000000D+00 x(2) = -0.77777777777777777778D+00 x(3) = -0.55555555555555555556D+00 x(4) = -0.33333333333333333333D+00 x(5) = -0.11111111111111111111D+00 x(6) = 0.11111111111111111111D+00 x(7) = 0.33333333333333333333D+00 x(8) = 0.55555555555555555556D+00 x(9) = 0.77777777777777777778D+00 x(10) = 1.00000000000000000000D+00 w(1) = 0.063772321428571428571D+00 w(2) = 0.35136160714285714286D+00 w(3) = 0.024107142857142857143D+00 w(4) = 0.43178571428571428571D+00 w(5) = 0.12897321428571428571D+00 w(6) = 0.12897321428571428571D+00 w(7) = 0.43178571428571428571D+00 w(8) = 0.024107142857142857143D+00 w(9) = 0.35136160714285714286D+00 w(10) = 0.063772321428571428571D+00 else if ( n .eq. 11 ) then c c 16067 / 299376 c 106300 / 299376 c - 48525 / 299376 c 272400 / 299376 c - 260550 / 299376 c 427368 / 299376 c - 260550 / 299376 c 272400 / 299376 c - 48525 / 299376 c 106300 / 299376 c 16067 / 299376 c x(1) = -1.00000000000000000000D+00 x(2) = -0.80000000000000000000D+00 x(3) = -0.60000000000000000000D+00 x(4) = -0.40000000000000000000D+00 x(5) = -0.20000000000000000000D+00 x(6) = 0.00000000000000000000D+00 x(7) = 0.20000000000000000000D+00 x(8) = 0.40000000000000000000D+00 x(9) = 0.60000000000000000000D+00 x(10) = 0.80000000000000000000D+00 x(11) = 1.00000000000000000000D+00 w(1) = 0.053668296723852279408D+00 w(2) = 0.35507188284966062744D+00 w(3) = -0.16208714125380792047D+00 w(4) = 0.90989257655924322591D+00 w(5) = -0.87031024531024531025D+00 w(6) = 1.4275292608625941959D+00 w(7) = -0.87031024531024531025D+00 w(8) = 0.90989257655924322591D+00 w(9) = -0.16208714125380792047D+00 w(10) = 0.35507188284966062744D+00 w(11) = 0.053668296723852279408D+00 else if ( n .eq. 12 ) then c c 2171465 / 43545600 c 13486539 / 43545600 c - 3237113 / 43545600 c 25226685 / 43545600 c - 9595542 / 43545600 c 15493566 / 43545600 c 15493566 / 43545600 c - 9595542 / 43545600 c 25226685 / 43545600 c - 3237113 / 43545600 c 13486539 / 43545600 c 2171465 / 43545600 c x(1) = -1.00000000000000000000D+00 x(2) = -0.81818181818181818182D+00 x(3) = -0.63636363636363636364D+00 x(4) = -0.45454545454545454545D+00 x(5) = -0.27272727272727272727D+00 x(6) = -0.090909090909090909091D+00 x(7) = 0.090909090909090909091D+00 x(8) = 0.27272727272727272727D+00 x(9) = 0.45454545454545454545D+00 x(10) = 0.63636363636363636364D+00 x(11) = 0.81818181818181818182D+00 x(12) = 1.00000000000000000000D+00 w(1) = 0.049866461823927101705D+00 w(2) = 0.30971071704144620811D+00 w(3) = -0.074338463587595532040D+00 w(4) = 0.57931650958994708995D+00 w(5) = -0.22035617835097001764D+00 w(6) = 0.35580095348324514991D+00 w(7) = 0.35580095348324514991D+00 w(8) = -0.22035617835097001764D+00 w(9) = 0.57931650958994708995D+00 w(10) = -0.074338463587595532040D+00 w(11) = 0.30971071704144620811D+00 w(12) = 0.049866461823927101705D+00 else if ( n .eq. 13 ) then c c 1364651 / 31531500 c 9903168 / 31531500 c - 7587864 / 31531500 c 35725120 / 31531500 c - 51491295 / 31531500 c 87516288 / 31531500 c - 87797136 / 31531500 c 87516288 / 31531500 c - 51491295 / 31531500 c 35725120 / 31531500 c - 7587864 / 31531500 c 9903168 / 31531500 c 1364651 / 31531500 c x(1) = -1.00000000000000000000D+00 x(2) = -0.83333333333333333333D+00 x(3) = -0.66666666666666666667D+00 x(4) = -0.50000000000000000000D+00 x(5) = -0.33333333333333333333D+00 x(6) = -0.16666666666666666667D+00 x(7) = 0.00000000000000000000D+00 x(8) = 0.16666666666666666667D+00 x(9) = 0.33333333333333333333D+00 x(10) = 0.50000000000000000000D+00 x(11) = 0.66666666666666666667D+00 x(12) = 0.83333333333333333333D+00 x(13) = 1.00000000000000000000D+00 w(1) = 0.043278974993260707546D+00 w(2) = 0.31407221350078492936D+00 w(3) = -0.24064392750107035821D+00 w(4) = 1.1329977958549387121D+00 w(5) = -1.6330112744398458684D+00 w(6) = 2.7755193378050520908D+00 w(7) = -2.7844262404262404262D+00 w(8) = 2.7755193378050520908D+00 w(9) = -1.6330112744398458684D+00 w(10) = 1.1329977958549387121D+00 w(11) = -0.24064392750107035821D+00 w(12) = 0.31407221350078492936D+00 w(13) = 0.043278974993260707546D+00 else if ( n .eq. 14 ) then c c 6137698213 / 150885504000 c 42194238652 / 150885504000 c - 23361540993 / 150885504000 c 116778274403 / 150885504000 c - 113219777650 / 150885504000 c 154424590209 / 150885504000 c - 32067978834 / 150885504000 c - 32067978834 / 150885504000 c 154424590209 / 150885504000 c - 113219777650 / 150885504000 c 116778274403 / 150885504000 c - 23361540993 / 150885504000 c 42194238652 / 150885504000 c 6137698213 / 150885504000 c x(1) = -1.00000000000000000000D+00 x(2) = -0.84615384615384615385D+00 x(3) = -0.69230769230769230769D+00 x(4) = -0.53846153846153846154D+00 x(5) = -0.38461538461538461538D+00 x(6) = -0.23076923076923076923D+00 x(7) = -0.076923076923076923077D+00 x(8) = 0.076923076923076923077D+00 x(9) = 0.23076923076923076923D+00 x(10) = 0.38461538461538461538D+00 x(11) = 0.53846153846153846154D+00 x(12) = 0.69230769230769230769D+00 x(13) = 0.84615384615384615385D+00 x(14) = 1.00000000000000000000D+00 w(1) = 0.040669438210247155353D+00 w(2) = 0.27975217053157074652D+00 w(3) = -0.15542374057682837445D+00 w(4) = 0.77579230848776566369D+00 w(5) = -0.75384763266423526013D+00 w(6) = 1.0273523591123107492D+00 w(7) = -0.21429490310083068020D+00 w(8) = -0.21429490310083068020D+00 w(9) = 1.0273523591123107492D+00 w(10) = -0.75384763266423526013D+00 w(11) = 0.77579230848776566369D+00 w(12) = -0.15542374057682837445D+00 w(13) = 0.27975217053157074652D+00 w(14) = 0.040669438210247155353D+00 else if ( n .eq. 15 ) then c c 90241897 / 2501928000 c 710986864 / 2501928000 c - 770720657 / 2501928000 c 3501442784 / 2501928000 c - 6625093363 / 2501928000 c 12630121616 / 2501928000 c - 16802270373 / 2501928000 c 19534438464 / 2501928000 c - 16802270373 / 2501928000 c 12630121616 / 2501928000 c - 6625093363 / 2501928000 c 3501442784 / 2501928000 c - 770720657 / 2501928000 c 710986864 / 2501928000 c 90241897 / 2501928000 c x(1) = -1.00000000000000000000D+00 x(2) = -0.85714285714285714286D+00 x(3) = -0.71428571428571428571D+00 x(4) = -0.57142857142857142857D+00 x(5) = -0.42857142857142857143D+00 x(6) = -0.28571428571428571429D+00 x(7) = -0.14285714285714285714D+00 x(8) = 0.00000000000000000000D+00 x(9) = 0.14285714285714285714D+00 x(10) = 0.28571428571428571429D+00 x(11) = 0.42857142857142857143D+00 x(12) = 0.57142857142857142857D+00 x(13) = 0.71428571428571428571D+00 x(14) = 0.85714285714285714286D+00 x(15) = 1.00000000000000000000D+00 w(1) = 0.036068942431596752584D+00 w(2) = 0.28417558938546592868D+00 w(3) = -0.30805069410470645039D+00 w(4) = 1.3994978208805369299D+00 w(5) = -2.6479952112930507992D+00 w(6) = 5.0481555088715582543D+00 w(7) = -6.7157289790113864188D+00 w(8) = 7.8077540456799716059D+00 w(9) = -6.7157289790113864188D+00 w(10) = 5.0481555088715582543D+00 w(11) = -2.6479952112930507992D+00 w(12) = 1.3994978208805369299D+00 w(13) = -0.30805069410470645039D+00 w(14) = 0.28417558938546592868D+00 w(15) = 0.036068942431596752584D+00 else if ( n .eq. 16 ) then c c 105930069 / 3099672576 c 796661595 / 3099672576 c - 698808195 / 3099672576 c 3143332755 / 3099672576 c - 4688522055 / 3099672576 c 7385654007 / 3099672576 c - 6000998415 / 3099672576 c 3056422815 / 3099672576 c 3056422815 / 3099672576 c - 6000998415 / 3099672576 c 7385654007 / 3099672576 c - 4688522055 / 3099672576 c 3143332755 / 3099672576 c - 698808195 / 3099672576 c 796661595 / 3099672576 c 105930069 / 3099672576 c x(1) = -1.00000000000000000000D+00 x(2) = -0.86666666666666666667D+00 x(3) = -0.73333333333333333333D+00 x(4) = -0.60000000000000000000D+00 x(5) = -0.46666666666666666667D+00 x(6) = -0.33333333333333333333D+00 x(7) = -0.20000000000000000000D+00 x(8) = -0.066666666666666666667D+00 x(9) = 0.066666666666666666667D+00 x(10) = 0.20000000000000000000D+00 x(11) = 0.33333333333333333333D+00 x(12) = 0.46666666666666666667D+00 x(13) = 0.60000000000000000000D+00 x(14) = 0.73333333333333333333D+00 x(15) = 0.86666666666666666667D+00 x(16) = 1.00000000000000000000D+00 w(1) = 0.034174599543251887002D+00 w(2) = 0.25701475735481036820D+00 w(3) = -0.22544581011901045383D+00 w(4) = 1.0140854164204471124D+00 w(5) = -1.5125862296882804695D+00 w(6) = 2.3827206990135980091D+00 w(7) = -1.9360104229924960952D+00 w(8) = 0.98604699046767964179D+00 w(9) = 0.98604699046767964179D+00 w(10) = -1.9360104229924960952D+00 w(11) = 2.3827206990135980091D+00 w(12) = -1.5125862296882804695D+00 w(13) = 1.0140854164204471124D+00 w(14) = -0.22544581011901045383D+00 w(15) = 0.25701475735481036820D+00 w(16) = 0.034174599543251887002D+00 else if ( n .eq. 17 ) then c c 15043611773 / 488462349375 c 127626606592 / 488462349375 c - 179731134720 / 488462349375 c 832211855360 / 488462349375 c - 1929498607520 / 488462349375 c 4177588893696 / 488462349375 c - 6806534407936 / 488462349375 c 9368875018240 / 488462349375 c - 10234238972220 / 488462349375 c 9368875018240 / 488462349375 c - 6806534407936 / 488462349375 c 4177588893696 / 488462349375 c - 1929498607520 / 488462349375 c 832211855360 / 488462349375 c - 179731134720 / 488462349375 c 127626606592 / 488462349375 c 15043611773 / 488462349375 c x(1) = -1.00000000000000000000D+00 x(2) = -0.87500000000000000000D+00 x(3) = -0.75000000000000000000D+00 x(4) = -0.62500000000000000000D+00 x(5) = -0.50000000000000000000D+00 x(6) = -0.37500000000000000000D+00 x(7) = -0.25000000000000000000D+00 x(8) = -0.12500000000000000000D+00 x(9) = 0.00000000000000000000D+00 x(10) = 0.12500000000000000000D+00 x(11) = 0.25000000000000000000D+00 x(12) = 0.37500000000000000000D+00 x(13) = 0.50000000000000000000D+00 x(14) = 0.62500000000000000000D+00 x(15) = 0.75000000000000000000D+00 x(16) = 0.87500000000000000000D+00 x(17) = 1.00000000000000000000D+00 w(1) = 0.030797894233299012495D+00 w(2) = 0.26128238288028031086D+00 w(3) = -0.36795289329867605622D+00 w(4) = 1.7037379778090086905D+00 w(5) = -3.9501480717783930427D+00 w(6) = 8.5525299934402953388D+00 w(7) = -13.934614237197880038D+00 w(8) = 19.180342211078732848D+00 w(9) = -20.951950514333334128D+00 w(10) = 19.180342211078732848D+00 w(11) = -13.934614237197880038D+00 w(12) = 8.5525299934402953388D+00 w(13) = -3.9501480717783930427D+00 w(14) = 1.7037379778090086905D+00 w(15) = -0.36795289329867605622D+00 w(16) = 0.26128238288028031086D+00 w(17) = 0.030797894233299012495D+00 else if ( n .eq. 18 ) then c c 55294720874657 / 1883051089920000 c 450185515446285 / 1883051089920000 c - 542023437008852 / 1883051089920000 c 2428636525764260 / 1883051089920000 c - 4768916800123440 / 1883051089920000 c 8855416648684984 / 1883051089920000 c - 10905371859796660 / 1883051089920000 c 10069615750132836 / 1883051089920000 c - 3759785974054070 / 1883051089920000 c - 3759785974054070 / 1883051089920000 c 10069615750132836 / 1883051089920000 c - 10905371859796660 / 1883051089920000 c 8855416648684984 / 1883051089920000 c - 4768916800123440 / 1883051089920000 c 2428636525764260 / 1883051089920000 c - 542023437008852 / 1883051089920000 c 450185515446285 / 1883051089920000 c 55294720874657 / 1883051089920000 c x(1) = -1.00000000000000000000D+00 x(2) = -0.88235294117647058824D+00 x(3) = -0.76470588235294117647D+00 x(4) = -0.64705882352941176471D+00 x(5) = -0.52941176470588235294D+00 x(6) = -0.41176470588235294118D+00 x(7) = -0.29411764705882352941D+00 x(8) = -0.17647058823529411765D+00 x(9) = -0.058823529411764705882D+00 x(10) = 0.058823529411764705882D+00 x(11) = 0.17647058823529411765D+00 x(12) = 0.29411764705882352941D+00 x(13) = 0.41176470588235294118D+00 x(14) = 0.52941176470588235294D+00 x(15) = 0.64705882352941176471D+00 x(16) = 0.76470588235294117647D+00 x(17) = 0.88235294117647058824D+00 x(18) = 1.00000000000000000000D+00 w(1) = 0.029364429446790078519D+00 w(2) = 0.23907238516051669677D+00 w(3) = -0.28784319231183443641D+00 w(4) = 1.2897348026109258587D+00 w(5) = -2.5325477495812627261D+00 w(6) = 4.7026959045817496499D+00 w(7) = -5.7913308450170443690D+00 w(8) = 5.3475000248456540826D+00 w(9) = -1.9966457597354948350D+00 w(10) = -1.9966457597354948350D+00 w(11) = 5.3475000248456540826D+00 w(12) = -5.7913308450170443690D+00 w(13) = 4.7026959045817496499D+00 w(14) = -2.5325477495812627261D+00 w(15) = 1.2897348026109258587D+00 w(16) = -0.28784319231183443641D+00 w(17) = 0.23907238516051669677D+00 w(18) = 0.029364429446790078519D+00 else if ( n .eq. 19 ) then c c 203732352169 / 7604556960000 c 1848730221900 / 7604556960000 c - 3212744374395 / 7604556960000 c 15529830312096 / 7604556960000 c - 42368630685840 / 7604556960000 c 103680563465808 / 7604556960000 c - 198648429867720 / 7604556960000 c 319035784479840 / 7604556960000 c - 419127951114198 / 7604556960000 c 461327344340680 / 7604556960000 c - 419127951114198 / 7604556960000 c 319035784479840 / 7604556960000 c - 198648429867720 / 7604556960000 c 103680563465808 / 7604556960000 c - 42368630685840 / 7604556960000 c 15529830312096 / 7604556960000 c - 3212744374395 / 7604556960000 c 1848730221900 / 7604556960000 c 203732352169 / 7604556960000 c x(1) = -1.00000000000000000000D+00 x(2) = -0.88888888888888888889D+00 x(3) = -0.77777777777777777778D+00 x(4) = -0.66666666666666666667D+00 x(5) = -0.55555555555555555556D+00 x(6) = -0.44444444444444444444D+00 x(7) = -0.33333333333333333333D+00 x(8) = -0.22222222222222222222D+00 x(9) = -0.11111111111111111111D+00 x(10) = 0.00000000000000000000D+00 x(11) = 0.11111111111111111111D+00 x(12) = 0.22222222222222222222D+00 x(13) = 0.33333333333333333333D+00 x(14) = 0.44444444444444444444D+00 x(15) = 0.55555555555555555556D+00 x(16) = 0.66666666666666666667D+00 x(17) = 0.77777777777777777778D+00 x(18) = 0.88888888888888888889D+00 x(19) = 1.00000000000000000000D+00 w(1) = 0.026790824664820447344D+00 w(2) = 0.24310820888374278151D+00 w(3) = -0.42247620621346493274D+00 w(4) = 2.0421742376029227612D+00 w(5) = -5.5714791681749728126D+00 w(6) = 13.634004454324976218D+00 w(7) = -26.122288374274995239D+00 w(8) = 41.953237533490708445D+00 w(9) = -55.115367445968607749D+00 w(10) = 60.664591871329740161D+00 w(11) = -55.115367445968607749D+00 w(12) = 41.953237533490708445D+00 w(13) = -26.122288374274995239D+00 w(14) = 13.634004454324976218D+00 w(15) = -5.5714791681749728126D+00 w(16) = 2.0421742376029227612D+00 w(17) = -0.42247620621346493274D+00 w(18) = 0.24310820888374278151D+00 w(19) = 0.026790824664820447344D+00 else if ( n .eq. 20 ) then c c 69028763155644023 / 2688996956405760000 c 603652082270808125 / 2688996956405760000 c - 926840515700222955 / 2688996956405760000 c 4301581538450500095 / 2688996956405760000 c - 10343692234243192788 / 2688996956405760000 c 22336420328479961316 / 2688996956405760000 c - 35331888421114781580 / 2688996956405760000 c 43920768370565135580 / 2688996956405760000 c - 37088370261379851390 / 2688996956405760000 c 15148337305921759574 / 2688996956405760000 c 15148337305921759574 / 2688996956405760000 c - 37088370261379851390 / 2688996956405760000 c 43920768370565135580 / 2688996956405760000 c - 35331888421114781580 / 2688996956405760000 c 22336420328479961316 / 2688996956405760000 c - 10343692234243192788 / 2688996956405760000 c 4301581538450500095 / 2688996956405760000 c - 926840515700222955 / 2688996956405760000 c 603652082270808125 / 2688996956405760000 c 69028763155644023 / 2688996956405760000 c x(1) = -1.00000000000000000000D+00 x(2) = -0.89473684210526315789D+00 x(3) = -0.78947368421052631579D+00 x(4) = -0.68421052631578947368D+00 x(5) = -0.57894736842105263158D+00 x(6) = -0.47368421052631578947D+00 x(7) = -0.36842105263157894737D+00 x(8) = -0.26315789473684210526D+00 x(9) = -0.15789473684210526316D+00 x(10) = -0.052631578947368421053D+00 x(11) = 0.052631578947368421053D+00 x(12) = 0.15789473684210526316D+00 x(13) = 0.26315789473684210526D+00 x(14) = 0.36842105263157894737D+00 x(15) = 0.47368421052631578947D+00 x(16) = 0.57894736842105263158D+00 x(17) = 0.68421052631578947368D+00 x(18) = 0.78947368421052631579D+00 x(19) = 0.89473684210526315789D+00 x(20) = 1.00000000000000000000D+00 w(1) = 0.025670822345560078100D+00 w(2) = 0.22448968595251886556D+00 w(3) = -0.34467890099030890987D+00 w(4) = 1.5996974366978074270D+00 w(5) = -3.8466730910952978835D+00 w(6) = 8.3065993344729824120D+00 w(7) = -13.139430424771119113D+00 w(8) = 16.333513604742678295D+00 w(9) = -13.792641220001198577D+00 w(10) = 5.6334527526463774045D+00 w(11) = 5.6334527526463774045D+00 w(12) = -13.792641220001198577D+00 w(13) = 16.333513604742678295D+00 w(14) = -13.139430424771119113D+00 w(15) = 8.3065993344729824120D+00 w(16) = -3.8466730910952978835D+00 w(17) = 1.5996974366978074270D+00 w(18) = -0.34467890099030890987D+00 w(19) = 0.22448968595251886556D+00 w(20) = 0.025670822345560078100D+00 else if ( n .eq. 21 ) then x(1) = -1.00000000000000000000D+00 x(2) = -0.90000000000000000000D+00 x(3) = -0.80000000000000000000D+00 x(4) = -0.70000000000000000000D+00 x(5) = -0.60000000000000000000D+00 x(6) = -0.50000000000000000000D+00 x(7) = -0.40000000000000000000D+00 x(8) = -0.30000000000000000000D+00 x(9) = -0.20000000000000000000D+00 x(10) = -0.10000000000000000000D+00 x(11) = 0.00000000000000000000D+00 x(12) = 0.10000000000000000000D+00 x(13) = 0.20000000000000000000D+00 x(14) = 0.30000000000000000000D+00 x(15) = 0.40000000000000000000D+00 x(16) = 0.50000000000000000000D+00 x(17) = 0.60000000000000000000D+00 x(18) = 0.70000000000000000000D+00 x(19) = 0.80000000000000000000D+00 x(20) = 0.90000000000000000000D+00 x(21) = 1.00000000000000000000D+00 w(1) = 0.023650546498063206389D+00 w(2) = 0.22827543528921394997D+00 w(3) = -0.47295674102285392846D+00 w(4) = 2.4123737869637513288D+00 w(5) = -7.5420634534306609355D+00 w(6) = 20.673596439879602287D+00 w(7) = -45.417631687959024596D+00 w(8) = 83.656114844387109207D+00 w(9) = -128.15055898030800930D+00 w(10) = 165.59456694494570344D+00 w(11) = -180.01073427048578932D+00 w(12) = 165.59456694494570344D+00 w(13) = -128.15055898030800930D+00 w(14) = 83.656114844387109207D+00 w(15) = -45.417631687959024596D+00 w(16) = 20.673596439879602287D+00 w(17) = -7.5420634534306609355D+00 w(18) = 2.4123737869637513288D+00 w(19) = -0.47295674102285392846D+00 w(20) = 0.22827543528921394997D+00 w(21) = 0.023650546498063206389D+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NCC_SET - Fatal error!' write ( *, '(a,i8)' ) ' Illegal value of N = ', n write ( *, '(a)' ) ' Legal values are 1 through 21.' stop end if return end subroutine r8mat_transpose ( m, n, a, at ) c*********************************************************************72 c cc r8mat_transpose() makes a transposed copy of a matrix. c c Discussion: c c An R8MAT is an array of R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 13 June 2011 c c Author: c c John Burkardt c c Parameters: c c Input, integer M, N, the number of rows and columns of the matrix A. c c Input, double precision A(N,N), the matrix to be transposed. c c Output, double precision AT(N,M), the matrix to be transposed. c implicit none integer m integer n double precision a(m,n) double precision at(n,m) integer i integer j do j = 1, m do i = 1, n at(i,j) = a(j,i) end do end do return end function r8vec_max ( n, a ) c*********************************************************************72 c cc r8vec_max() returns the maximum value in an R8VEC. c c Discussion: c c An R8VEC is a vector of R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 May 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries in the array. c c Input, double precision A(N), the array. c c Output, double precision R8VEC_MAX, the value of the largest entry. c implicit none integer n double precision a(n) integer i double precision r8_huge parameter ( r8_huge = 1.79769313486231571D+308 ) double precision r8vec_max double precision value value = - r8_huge do i = 1, n value = max ( value, a(i) ) end do r8vec_max = value return end function r8vec_sum ( n, v1 ) c*********************************************************************72 c cc r8vec_sum() sums the entries of an R8VEC. c c Discussion: c c An R8VEC is a vector of R8's. c c In FORTRAN90, the system routine SUM should be called c directly. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 22 July 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the dimension of the vectors. c c Input, double precision V1(N), the vector. c c Output, double precision R8VEC_SUM, the sum of the entries. c implicit none integer n integer i double precision r8vec_sum double precision v1(n) double precision value value = 0.0D+00 do i = 1, n value = value + v1(i) end do r8vec_sum = value return end subroutine r8vec_uniform_ab ( n, a, b, seed, r ) c*********************************************************************72 c cc r8vec_uniform_ab() returns a scaled pseudorandom R8VEC. c c Discussion: c c Each dimension ranges from A to B. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 29 January 2005 c c Author: c c John Burkardt c c Reference: c c Paul Bratley, Bennett Fox, Linus Schrage, c A Guide to Simulation, c Second Edition, c Springer, 1987, c ISBN: 0387964673, c LC: QA76.9.C65.B73. c c Bennett Fox, c Algorithm 647: c Implementation and Relative Efficiency of Quasirandom c Sequence Generators, c ACM Transactions on Mathematical Software, c Volume 12, Number 4, December 1986, pages 362-376. c c Pierre L'Ecuyer, c Random Number Generation, c in Handbook of Simulation, c edited by Jerry Banks, c Wiley, 1998, c ISBN: 0471134031, c LC: T57.62.H37. c c Peter Lewis, Allen Goodman, James Miller, c A Pseudo-Random Number Generator for the System/360, c IBM Systems Journal, c Volume 8, Number 2, 1969, pages 136-143. c c Parameters: c c Input, integer N, the number of entries in the vector. c c Input, double precision A, B, the lower and upper limits. c c Input/output, integer SEED, the "seed" value, which should NOT be 0. c On output, SEED has been updated. c c Output, double precision R(N), the vector of pseudorandom values. c implicit none integer n double precision a double precision b integer i integer i4_huge parameter ( i4_huge = 2147483647 ) integer k integer seed double precision r(n) if ( seed .eq. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8VEC_UNIFORM_AB - Fatal error!' write ( *, '(a)' ) ' Input value of SEED = 0.' stop 1 end if do i = 1, n k = seed / 127773 seed = 16807 * ( seed - k * 127773 ) - k * 2836 if ( seed .lt. 0 ) then seed = seed + i4_huge end if r(i) = a + ( b - a ) * dble ( seed ) * 4.656612875D-10 end do return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end subroutine weights_ls ( d, a, b, n, x, w ) c*********************************************************************72 c cc WEIGHTS_LS computes weights for a least squares quadrature rule. c c Discussion: c c We suppose that we wish to estimate integrals of functions f(x) c over the interval [a,b], given values of the function at N points X, c and that we wish to do this in such a way that our estimate is c based on the least squares polynomial approximant p(x) of degree D. c c The resulting quadrature rule will only be exact if N = D + 1 c and f(x) is a polynomial of degree D or less. c c This function is based in part on ideas of Mac Hyman, of Tulane c University. c c Method: c c The approximant p(x) to f(x) has coefficients C(1:D+1). c The values of p(x) at the abscissas are found by V * C = P. c The values of f(x) at the abscissas are in the vector F. c We want to minimize ||P-F|| = ||V * C - F||. c We set up the D+1 x D+1 system V' * V * C = V' * F. c c Now, for quadrature, we have Q * C = W * F, c so formally Q * inv ( V' * V ) * V' = W. c c Hence, to solve for the weights W, we solve two Vandermonde systems: c c 1) Solve Q' = V' * S for S, an N vector. c 2) Solve S = V * R for R, a D+1 vector. c 3) Compute W = R' * V' for W, an N vector. c c so that, by construction (and writing inv() for matrices which c actually don't have inverses): c c W = R' * V' c W' = V * R c W' = V * inv ( V ) * S c W' = V * inv ( V ) * inv ( V' ) * Q' c W = Q * inv ( V ) * inv ( V' ) * V' c W = Q * inv ( V' * V ) * V' c c Thus, knowing the values Q and the abscissas X, we can construct c V and solve for W, and the value of W depends (as it should) only c on N and X, but not on f(x). c c Having W, we can estimate the integral of any f(x) c if we have its values F(1:n) at the points X(1:n), by the formula: c c I(f) \approx Q(f) = sum ( 1 <= i <= n ) W(i) * F(i) c c A simple test is to use N points associated with a known quadrature c rule, and specify D = N-1. Assuming the rule is exact for polynomials c of degree D, this procedure should determine the same weights as c those used for the given quadrature rule. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 15 April 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer D, the desired degree of the polynomial c approximant. 0 <= D. c c Input, double precision A, B, the left and right endpoints of c the interval. c c Input, integer N, the number of data points or abscissas. c c Input, double precision X(N), the data points. c c Output, double precision W(N), the quadrature weights. c implicit none integer d integer n double precision a double precision ai double precision b double precision bi integer i integer j double precision q(d+1) double precision r(d+1) double precision s(n) double precision v(n,d+1) double precision vt(d+1,n) double precision w(n) double precision x(n) c c Integrate the basis functions x^0, x^1, ..., x^d over [a,b]. c bi = b ai = a do i = 1, d + 1 q(i) = ( bi - ai ) / dble ( i ) ai = ai * a bi = bi * b end do c c Form the N by D+1 Vandermonde matrix V. c do i = 1, n v(i,1) = 1.0D+00 do j = 2, d + 1 v(i,j) = v(i,j-1) * x(i) end do end do ! ! Solve the D+1 by N system V' * S = Q for vector S of length N. ! call r8mat_transpose ( n, d + 1, v, vt ) call qr_solve ( d + 1, n, vt, q, s ) ! ! Solve the N by D+1 system V * R = S for vector R, of length D+1. ! call qr_solve ( n, d + 1, v, s, r ) ! ! Compute W = R' * V'; ! call r8mat_mv ( n, d + 1, v, r, w ) return end