subroutine ncc_set ( n, x, w ) !*****************************************************************************80 ! !! NCC_SET sets abscissas and weights for Newton-Cotes closed quadrature. ! ! Discussion: ! ! The Newton-Cotes closed rules use equally spaced abscissas, and ! hence may be used with tabulated function data. ! ! The rules are called "closed" because they include the endpoints. ! As a favor, we include an order 1 rule, the midpoint rule, even ! though this does not satisfy the requirement that the endpoints ! be included! ! ! The higher order rules involve negative weights. These can produce ! loss of accuracy due to the subtraction of large, nearly equal quantities. ! ! N = 1 is the midpoint rule (and is not really an NCC rule!) ! N = 2 is the trapezoidal rule. ! N = 3 is Simpson's rule. ! N = 4 is Simpson's 3/8 rule. ! N = 5 is Bode's rule. ! ! The Kopal reference for N = 12 lists ! W(6) = 15494566.0D+00 / 43545600.0D+00 ! but this results in a set of coeffients that don't add up to 2. ! The correct value is ! W(6) = 15493566.0D+00 / 43545600.0. ! ! The integral: ! ! Integral ( -1 <= X <= 1 ) F(X) dx ! ! The quadrature rule: ! ! Sum ( 1 <= I <= N ) W(I) * F ( X(I) ) ! ! In Mathematica, the Newton-Cotes closed weights of order N in ! the interval [-1,+1] can be computed by: ! ! Needs["NumericalDifferentialEquationAnalysis`"] ! NewtonCotesWeights [ n, -1, 1, QuadratureType -> Closed ] ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 April 2010 ! ! 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. ! ! Johnson, ! Quarterly Journal of Mathematics, ! Volume 46, Number 52, 1915. ! ! Zdenek Kopal, ! Numerical Analysis, ! John Wiley, 1955, ! LC: QA297.K6. ! ! Stephen Wolfram, ! The Mathematica Book, ! Fourth Edition, ! Cambridge University Press, 1999, ! ISBN: 0-521-64314-7, ! LC: QA76.95.W65. ! ! Daniel Zwillinger, editor, ! CRC Standard Mathematical Tables and Formulae, ! 30th Edition, ! CRC Press, 1996, ! ISBN: 0-8493-2479-3, ! LC: QA47.M315. ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the order. ! 1 <= N <= 21. ! ! Output, real ( kind = 8 ) X(N), the abscissas. ! ! Output, real ( kind = 8 ) W(N), the weights. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) w(n) real ( kind = 8 ) x(n) if ( n == 1 ) then ! ! 2 ! x(1) = 0.00000000000000000000D+00 w(1) = 2.00000000000000000000D+00 else if ( n == 2 ) then ! ! 1 ! 1 ! x(1) = -1.00000000000000000000D+00 x(2) = 1.00000000000000000000D+00 w(1) = 1.00000000000000000000D+00 w(2) = 1.00000000000000000000D+00 else if ( n == 3 ) then ! ! 1 / 3 ! 4 / 3 ! 1 / 3 ! 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 == 4 ) then ! ! 1 / 4 ! 3 / 4 ! 3 / 4 ! 1 / 4 ! 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 == 5 ) then ! ! 7 / 45 ! 32 / 45 ! 12 / 45 ! 32 / 45 ! 7 / 45 ! 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 == 6 ) then ! ! 19 / 144 ! 75 / 144 ! 50 / 144 ! 50 / 144 ! 75 / 144 ! 19 / 144 ! 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 == 7 ) then ! ! 41 / 420 ! 216 / 420 ! 27 / 420 ! 272 / 420 ! 27 / 420 ! 216 / 420 ! 41 / 420 ! 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 == 8 ) then ! ! 751 / 8640 ! 3577 / 8640 ! 1323 / 8640 ! 2989 / 8640 ! 2989 / 8640 ! 1323 / 8640 ! 3577 / 8640 ! 751 / 8640 ! 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 == 9 ) then ! ! 989 / 14175 ! 5888 / 14175 ! -928 / 14175 ! 10496 / 14175 ! -4540 / 14175 ! 10496 / 14175 ! -928 / 14175 ! 5888 / 14175 ! 989 / 14175 ! 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 == 10 ) then ! ! 2857 / 44800 ! 15741 / 44800 ! 1080 / 44800 ! 19344 / 44800 ! 5778 / 44800 ! 5778 / 44800 ! 19344 / 44800 ! 1080 / 44800 ! 15741 / 44800 ! 2857 / 44800 ! 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 == 11 ) then ! ! 16067 / 299376 ! 106300 / 299376 ! - 48525 / 299376 ! 272400 / 299376 ! - 260550 / 299376 ! 427368 / 299376 ! - 260550 / 299376 ! 272400 / 299376 ! - 48525 / 299376 ! 106300 / 299376 ! 16067 / 299376 ! 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 == 12 ) then ! ! 2171465 / 43545600 ! 13486539 / 43545600 ! - 3237113 / 43545600 ! 25226685 / 43545600 ! - 9595542 / 43545600 ! 15493566 / 43545600 ! 15493566 / 43545600 ! - 9595542 / 43545600 ! 25226685 / 43545600 ! - 3237113 / 43545600 ! 13486539 / 43545600 ! 2171465 / 43545600 ! 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 == 13 ) then ! ! 1364651 / 31531500 ! 9903168 / 31531500 ! - 7587864 / 31531500 ! 35725120 / 31531500 ! - 51491295 / 31531500 ! 87516288 / 31531500 ! - 87797136 / 31531500 ! 87516288 / 31531500 ! - 51491295 / 31531500 ! 35725120 / 31531500 ! - 7587864 / 31531500 ! 9903168 / 31531500 ! 1364651 / 31531500 ! 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 == 14 ) then ! ! 6137698213 / 150885504000 ! 42194238652 / 150885504000 ! - 23361540993 / 150885504000 ! 116778274403 / 150885504000 ! - 113219777650 / 150885504000 ! 154424590209 / 150885504000 ! - 32067978834 / 150885504000 ! - 32067978834 / 150885504000 ! 154424590209 / 150885504000 ! - 113219777650 / 150885504000 ! 116778274403 / 150885504000 ! - 23361540993 / 150885504000 ! 42194238652 / 150885504000 ! 6137698213 / 150885504000 ! 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 == 15 ) then ! ! 90241897 / 2501928000 ! 710986864 / 2501928000 ! - 770720657 / 2501928000 ! 3501442784 / 2501928000 ! - 6625093363 / 2501928000 ! 12630121616 / 2501928000 ! - 16802270373 / 2501928000 ! 19534438464 / 2501928000 ! - 16802270373 / 2501928000 ! 12630121616 / 2501928000 ! - 6625093363 / 2501928000 ! 3501442784 / 2501928000 ! - 770720657 / 2501928000 ! 710986864 / 2501928000 ! 90241897 / 2501928000 ! 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 == 16 ) then ! ! 105930069 / 3099672576 ! 796661595 / 3099672576 ! - 698808195 / 3099672576 ! 3143332755 / 3099672576 ! - 4688522055 / 3099672576 ! 7385654007 / 3099672576 ! - 6000998415 / 3099672576 ! 3056422815 / 3099672576 ! 3056422815 / 3099672576 ! - 6000998415 / 3099672576 ! 7385654007 / 3099672576 ! - 4688522055 / 3099672576 ! 3143332755 / 3099672576 ! - 698808195 / 3099672576 ! 796661595 / 3099672576 ! 105930069 / 3099672576 ! 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 == 17 ) then ! ! 15043611773 / 488462349375 ! 127626606592 / 488462349375 ! - 179731134720 / 488462349375 ! 832211855360 / 488462349375 ! - 1929498607520 / 488462349375 ! 4177588893696 / 488462349375 ! - 6806534407936 / 488462349375 ! 9368875018240 / 488462349375 ! - 10234238972220 / 488462349375 ! 9368875018240 / 488462349375 ! - 6806534407936 / 488462349375 ! 4177588893696 / 488462349375 ! - 1929498607520 / 488462349375 ! 832211855360 / 488462349375 ! - 179731134720 / 488462349375 ! 127626606592 / 488462349375 ! 15043611773 / 488462349375 ! 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 == 18 ) then ! ! 55294720874657 / 1883051089920000 ! 450185515446285 / 1883051089920000 ! - 542023437008852 / 1883051089920000 ! 2428636525764260 / 1883051089920000 ! - 4768916800123440 / 1883051089920000 ! 8855416648684984 / 1883051089920000 ! - 10905371859796660 / 1883051089920000 ! 10069615750132836 / 1883051089920000 ! - 3759785974054070 / 1883051089920000 ! - 3759785974054070 / 1883051089920000 ! 10069615750132836 / 1883051089920000 ! - 10905371859796660 / 1883051089920000 ! 8855416648684984 / 1883051089920000 ! - 4768916800123440 / 1883051089920000 ! 2428636525764260 / 1883051089920000 ! - 542023437008852 / 1883051089920000 ! 450185515446285 / 1883051089920000 ! 55294720874657 / 1883051089920000 ! 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 == 19 ) then ! ! 203732352169 / 7604556960000 ! 1848730221900 / 7604556960000 ! - 3212744374395 / 7604556960000 ! 15529830312096 / 7604556960000 ! - 42368630685840 / 7604556960000 ! 103680563465808 / 7604556960000 ! - 198648429867720 / 7604556960000 ! 319035784479840 / 7604556960000 ! - 419127951114198 / 7604556960000 ! 461327344340680 / 7604556960000 ! - 419127951114198 / 7604556960000 ! 319035784479840 / 7604556960000 ! - 198648429867720 / 7604556960000 ! 103680563465808 / 7604556960000 ! - 42368630685840 / 7604556960000 ! 15529830312096 / 7604556960000 ! - 3212744374395 / 7604556960000 ! 1848730221900 / 7604556960000 ! 203732352169 / 7604556960000 ! 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 == 20 ) then ! ! 69028763155644023 / 2688996956405760000 ! 603652082270808125 / 2688996956405760000 ! - 926840515700222955 / 2688996956405760000 ! 4301581538450500095 / 2688996956405760000 ! - 10343692234243192788 / 2688996956405760000 ! 22336420328479961316 / 2688996956405760000 ! - 35331888421114781580 / 2688996956405760000 ! 43920768370565135580 / 2688996956405760000 ! - 37088370261379851390 / 2688996956405760000 ! 15148337305921759574 / 2688996956405760000 ! 15148337305921759574 / 2688996956405760000 ! - 37088370261379851390 / 2688996956405760000 ! 43920768370565135580 / 2688996956405760000 ! - 35331888421114781580 / 2688996956405760000 ! 22336420328479961316 / 2688996956405760000 ! - 10343692234243192788 / 2688996956405760000 ! 4301581538450500095 / 2688996956405760000 ! - 926840515700222955 / 2688996956405760000 ! 603652082270808125 / 2688996956405760000 ! 69028763155644023 / 2688996956405760000 ! 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 == 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 weights_ls ( d, a, b, n, x, w ) !*****************************************************************************80 ! !! WEIGHTS_LS computes weights for a least squares quadrature rule. ! ! Discussion: ! ! We suppose that we wish to estimate integrals of functions f(x) ! over the interval [a,b], given values of the function at N points X, ! and that we wish to do this in such a way that our estimate is ! based on the least squares polynomial approximant p(x) of degree D. ! ! The resulting quadrature rule will only be exact if N = D + 1 ! and f(x) is a polynomial of degree D or less. ! ! This function is based in part on ideas of Mac Hyman, of Tulane ! University. ! ! Method: ! ! The approximant p(x) to f(x) has coefficients C(1:D+1). ! The values of p(x) at the abscissas are found by V * C = P. ! The values of f(x) at the abscissas are in the vector F. ! We want to minimize ||P-F|| = ||V * C - F||. ! We set up the D+1 x D+1 system V' * V * C = V' * F. ! ! Now, for quadrature, we have Q * C = W * F, ! so formally Q * inv ( V' * V ) * V' = W. ! ! Hence, to solve for the weights W, we solve two Vandermonde systems: ! ! 1) Solve Q' = V' * S for S, an N vector. ! 2) Solve S = V * R for R, a D+1 vector. ! 3) Compute W = R' * V' for W, an N vector. ! ! so that, by construction (and writing inv() for matrices which ! actually don't have inverses): ! ! W = R' * V' ! W' = V * R ! W' = V * inv ( V ) * S ! W' = V * inv ( V ) * inv ( V' ) * Q' ! W = Q * inv ( V ) * inv ( V' ) * V' ! W = Q * inv ( V' * V ) * V' ! ! Thus, knowing the values Q and the abscissas X, we can construct ! V and solve for W, and the value of W depends (as it should) only ! on N and X, but not on f(x). ! ! Having W, we can estimate the integral of any f(x) ! if we have its values F(1:n) at the points X(1:n), by the formula: ! ! I(f) \approx Q(f) = sum ( 1 <= i <= n ) W(i) * F(i) ! ! A simple test is to use N points associated with a known quadrature ! rule, and specify D = N-1. Assuming the rule is exact for polynomials ! of degree D, this procedure should determine the same weights as ! those used for the given quadrature rule. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 April 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) D, the desired degree of the polynomial ! approximant. 0 <= D. ! ! Input, real ( kind = 8 ) A, B, the left and right endpoints of ! the interval. ! ! Input, integer ( kind = 4 ) N, the number of data points or abscissas. ! ! Input, real ( kind = 8 ) X(N), the data points. ! ! Output, real ( kind = 8 ) W(N), the quadrature weights. ! implicit none integer ( kind = 4 ) d integer ( kind = 4 ) n real ( kind = 8 ) a real ( kind = 8 ) ai real ( kind = 8 ) b real ( kind = 8 ) bi integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) q(d+1) real ( kind = 8 ) r(d+1) real ( kind = 8 ) s(n) real ( kind = 8 ) v(n,d+1) real ( kind = 8 ) vt(d+1,n) real ( kind = 8 ) w(n) real ( kind = 8 ) x(n) ! ! Integrate the basis functions x^0, x^1, ..., x^d over [a,b]. ! bi = b ai = a do i = 1, d + 1 q(i) = ( bi - ai ) / real ( i, kind = 8 ) ai = ai * a bi = bi * b end do ! ! Form the N by D+1 Vandermonde matrix V. ! v(1:n,1) = 1.0D+00 do j = 2, d + 1 v(1:n,j) = v(1:n,j-1) * x(1:n) end do ! ! Solve the D+1 by N system V' * S = Q for vector S of length N. ! vt = transpose ( v ) 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