# include # include # include # include # include # include # include "quadrature_least_squares.h" # include "qr_solve.h" # include "r8lib.h" /******************************************************************************/ void ncc_set ( int order, double x[], double w[] ) /******************************************************************************/ /* Purpose: NCC_SET sets abscissas and weights for closed Newton-Cotes quadrature. Discussion: The closed Newton-Cotes 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. ORDER = 1 is the midpoint rule (and is not really an NCC rule.) ORDER = 2 is the trapezoidal rule. ORDER = 3 is Simpson's rule. ORDER = 4 is Simpson's 3/8 rule. ORDER = 5 is Bode's rule. The Kopal reference for ORDER = 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 <= ORDER ) W(I) * F ( X(I) ) In Mathematica, the Newton-Cotes closed weights can be computed by: Needs["NumericalDifferentialEquationAnalysis`"] NewtonCotesWeights [ order, -1, 1, QuadratureType -> Closed ] Licensing: This code is distributed under the MIT license. Modified: 20 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, int ORDER, the order. ORDER must be between 1 and 21. Output, double X[ORDER], the abscissas. Output, double W[ORDER], the weights. */ { if ( order == 1 ) { /* 2 */ x[0] = 0.00000000000000000000; w[0] = 2.00000000000000000000; } else if ( order == 2 ) { /* 1 1 */ x[0] = -1.00000000000000000000; x[1] = 1.00000000000000000000; w[0] = 1.00000000000000000000; w[1] = 1.00000000000000000000; } else if ( order == 3 ) { /* 1 / 3 4 / 3 1 / 3 */ x[0] = -1.00000000000000000000; x[1] = 0.00000000000000000000; x[2] = 1.00000000000000000000; w[0] = 0.33333333333333333333; w[1] = 1.33333333333333333333; w[2] = 0.33333333333333333333; } else if ( order == 4 ) { /* 1 / 4 3 / 4 3 / 4 1 / 4 */ x[0] = -1.00000000000000000000; x[1] = -0.33333333333333333333; x[2] = 0.33333333333333333333; x[3] = 1.00000000000000000000; w[0] = 0.25000000000000000000; w[1] = 0.75000000000000000000; w[2] = 0.75000000000000000000; w[3] = 0.25000000000000000000; } else if ( order == 5 ) { /* 7 / 45 32 / 45 12 / 45 32 / 45 7 / 45 */ x[0] = -1.00000000000000000000; x[1] = -0.50000000000000000000; x[2] = 0.00000000000000000000; x[3] = 0.50000000000000000000; x[4] = 1.00000000000000000000; w[0] = 0.15555555555555555556; w[1] = 0.71111111111111111111; w[2] = 0.26666666666666666667; w[3] = 0.71111111111111111111; w[4] = 0.15555555555555555556; } else if ( order == 6 ) { /* 19 / 144 75 / 144 50 / 144 50 / 144 75 / 144 19 / 144 */ x[0] = -1.00000000000000000000; x[1] = -0.60000000000000000000; x[2] = -0.20000000000000000000; x[3] = 0.20000000000000000000; x[4] = 0.60000000000000000000; x[5] = 1.00000000000000000000; w[0] = 0.13194444444444444444; w[1] = 0.52083333333333333333; w[2] = 0.34722222222222222222; w[3] = 0.34722222222222222222; w[4] = 0.52083333333333333333; w[5] = 0.13194444444444444444; } else if ( order == 7 ) { /* 41 / 420 216 / 420 27 / 420 272 / 420 27 / 420 216 / 420 41 / 420 */ x[0] = -1.00000000000000000000; x[1] = -0.66666666666666666667; x[2] = -0.33333333333333333333; x[3] = 0.00000000000000000000; x[4] = 0.33333333333333333333; x[5] = 0.66666666666666666667; x[6] = 1.00000000000000000000; w[0] = 0.097619047619047619048; w[1] = 0.51428571428571428571; w[2] = 0.064285714285714285714; w[3] = 0.64761904761904761905; w[4] = 0.064285714285714285714; w[5] = 0.51428571428571428571; w[6] = 0.097619047619047619048; } else if ( order == 8 ) { /* 751 / 8640 3577 / 8640 1323 / 8640 2989 / 8640 2989 / 8640 1323 / 8640 3577 / 8640 751 / 8640 */ x[0] = -1.00000000000000000000; x[1] = -0.71428571428571428571; x[2] = -0.42857142857142857143; x[3] = -0.14285714285714285714; x[4] = 0.14285714285714285714; x[5] = 0.42857142857142857143; x[6] = 0.71428571428571428571; x[7] = 1.00000000000000000000; w[0] = 0.086921296296296296296; w[1] = 0.41400462962962962963; w[2] = 0.15312500000000000000; w[3] = 0.34594907407407407407; w[4] = 0.34594907407407407407; w[5] = 0.15312500000000000000; w[6] = 0.41400462962962962963; w[7] = 0.086921296296296296296; } else if ( order == 9 ) { /* 989 / 14175 5888 / 14175 -928 / 14175 10496 / 14175 -4540 / 14175 10496 / 14175 -928 / 14175 5888 / 14175 989 / 14175 */ x[0] = -1.00000000000000000000; x[1] = -0.75000000000000000000; x[2] = -0.50000000000000000000; x[3] = -0.25000000000000000000; x[4] = 0.00000000000000000000; x[5] = 0.25000000000000000000; x[6] = 0.50000000000000000000; x[7] = 0.75000000000000000000; x[8] = 1.00000000000000000000; w[0] = 0.069770723104056437390; w[1] = 0.41537918871252204586; w[2] = -0.065467372134038800705; w[3] = 0.74045855379188712522; w[4] = -0.32028218694885361552; w[5] = 0.74045855379188712522; w[6] = -0.065467372134038800705; w[7] = 0.41537918871252204586; w[8] = 0.069770723104056437390; } else if ( order == 10 ) { /* 2857 / 44800 15741 / 44800 1080 / 44800 19344 / 44800 5778 / 44800 5778 / 44800 19344 / 44800 1080 / 44800 15741 / 44800 2857 / 44800 */ x[0] = -1.00000000000000000000; x[1] = -0.77777777777777777778; x[2] = -0.55555555555555555556; x[3] = -0.33333333333333333333; x[4] = -0.11111111111111111111; x[5] = 0.11111111111111111111; x[6] = 0.33333333333333333333; x[7] = 0.55555555555555555556; x[8] = 0.77777777777777777778; x[9] = 1.00000000000000000000; w[0] = 0.063772321428571428571; w[1] = 0.35136160714285714286; w[2] = 0.024107142857142857143; w[3] = 0.43178571428571428571; w[4] = 0.12897321428571428571; w[5] = 0.12897321428571428571; w[6] = 0.43178571428571428571; w[7] = 0.024107142857142857143; w[8] = 0.35136160714285714286; w[9] = 0.063772321428571428571; } else if ( order == 11 ) { /* 16067 / 299376 106300 / 299376 - 48525 / 299376 272400 / 299376 - 260550 / 299376 427368 / 299376 - 260550 / 299376 272400 / 299376 - 48525 / 299376 106300 / 299376 16067 / 299376 */ x[0] = -1.00000000000000000000; x[1] = -0.80000000000000000000; x[2] = -0.60000000000000000000; x[3] = -0.40000000000000000000; x[4] = -0.20000000000000000000; x[5] = 0.00000000000000000000; x[6] = 0.20000000000000000000; x[7] = 0.40000000000000000000; x[8] = 0.60000000000000000000; x[9] = 0.80000000000000000000; x[10] = 1.00000000000000000000; w[0] = 0.053668296723852279408; w[1] = 0.35507188284966062744; w[2] = -0.16208714125380792047; w[3] = 0.90989257655924322591; w[4] = -0.87031024531024531025; w[5] = 1.4275292608625941959; w[6] = -0.87031024531024531025; w[7] = 0.90989257655924322591; w[8] = -0.16208714125380792047; w[9] = 0.35507188284966062744; w[10] = 0.053668296723852279408; } else if ( order == 12 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.81818181818181818182; x[2] = -0.63636363636363636364; x[3] = -0.45454545454545454545; x[4] = -0.27272727272727272727; x[5] = -0.090909090909090909091; x[6] = 0.090909090909090909091; x[7] = 0.27272727272727272727; x[8] = 0.45454545454545454545; x[9] = 0.63636363636363636364; x[10] = 0.81818181818181818182; x[11] = 1.00000000000000000000; w[0] = 0.049866461823927101705; w[1] = 0.30971071704144620811; w[2] = -0.074338463587595532040; w[3] = 0.57931650958994708995; w[4] = -0.22035617835097001764; w[5] = 0.35580095348324514991; w[6] = 0.35580095348324514991; w[7] = -0.22035617835097001764; w[8] = 0.57931650958994708995; w[9] = -0.074338463587595532040; w[10] = 0.30971071704144620811; w[11] = 0.049866461823927101705; } else if ( order == 13 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.83333333333333333333; x[2] = -0.66666666666666666667; x[3] = -0.50000000000000000000; x[4] = -0.33333333333333333333; x[5] = -0.16666666666666666667; x[6] = 0.00000000000000000000; x[7] = 0.16666666666666666667; x[8] = 0.33333333333333333333; x[9] = 0.50000000000000000000; x[10] = 0.66666666666666666667; x[11] = 0.83333333333333333333; x[12] = 1.00000000000000000000; w[0] = 0.043278974993260707546; w[1] = 0.31407221350078492936; w[2] = -0.24064392750107035821; w[3] = 1.1329977958549387121; w[4] = -1.6330112744398458684; w[5] = 2.7755193378050520908; w[6] = -2.7844262404262404262; w[7] = 2.7755193378050520908; w[8] = -1.6330112744398458684; w[9] = 1.1329977958549387121; w[10] = -0.24064392750107035821; w[11] = 0.31407221350078492936; w[12] = 0.043278974993260707546; } else if ( order == 14 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.84615384615384615385; x[2] = -0.69230769230769230769; x[3] = -0.53846153846153846154; x[4] = -0.38461538461538461538; x[5] = -0.23076923076923076923; x[6] = -0.076923076923076923077; x[7] = 0.076923076923076923077; x[8] = 0.23076923076923076923; x[9] = 0.38461538461538461538; x[10] = 0.53846153846153846154; x[11] = 0.69230769230769230769; x[12] = 0.84615384615384615385; x[13] = 1.00000000000000000000; w[0] = 0.040669438210247155353; w[1] = 0.27975217053157074652; w[2] = -0.15542374057682837445; w[3] = 0.77579230848776566369; w[4] = -0.75384763266423526013; w[5] = 1.0273523591123107492; w[6] = -0.21429490310083068020; w[7] = -0.21429490310083068020; w[8] = 1.0273523591123107492; w[9] = -0.75384763266423526013; w[10] = 0.77579230848776566369; w[11] = -0.15542374057682837445; w[12] = 0.27975217053157074652; w[13] = 0.040669438210247155353; } else if ( order == 15 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.85714285714285714286; x[2] = -0.71428571428571428571; x[3] = -0.57142857142857142857; x[4] = -0.42857142857142857143; x[5] = -0.28571428571428571429; x[6] = -0.14285714285714285714; x[7] = 0.00000000000000000000; x[8] = 0.14285714285714285714; x[9] = 0.28571428571428571429; x[10] = 0.42857142857142857143; x[11] = 0.57142857142857142857; x[12] = 0.71428571428571428571; x[13] = 0.85714285714285714286; x[14] = 1.00000000000000000000; w[0] = 0.036068942431596752584; w[1] = 0.28417558938546592868; w[2] = -0.30805069410470645039; w[3] = 1.3994978208805369299; w[4] = -2.6479952112930507992; w[5] = 5.0481555088715582543; w[6] = -6.7157289790113864188; w[7] = 7.8077540456799716059; w[8] = -6.7157289790113864188; w[9] = 5.0481555088715582543; w[10] = -2.6479952112930507992; w[11] = 1.3994978208805369299; w[12] = -0.30805069410470645039; w[13] = 0.28417558938546592868; w[14] = 0.036068942431596752584; } else if ( order == 16 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.86666666666666666667; x[2] = -0.73333333333333333333; x[3] = -0.60000000000000000000; x[4] = -0.46666666666666666667; x[5] = -0.33333333333333333333; x[6] = -0.20000000000000000000; x[7] = -0.066666666666666666667; x[8] = 0.066666666666666666667; x[9] = 0.20000000000000000000; x[10] = 0.33333333333333333333; x[11] = 0.46666666666666666667; x[12] = 0.60000000000000000000; x[13] = 0.73333333333333333333; x[14] = 0.86666666666666666667; x[15] = 1.00000000000000000000; w[0] = 0.034174599543251887002; w[1] = 0.25701475735481036820; w[2] = -0.22544581011901045383; w[3] = 1.0140854164204471124; w[4] = -1.5125862296882804695; w[5] = 2.3827206990135980091; w[6] = -1.9360104229924960952; w[7] = 0.98604699046767964179; w[8] = 0.98604699046767964179; w[9] = -1.9360104229924960952; w[10] = 2.3827206990135980091; w[11] = -1.5125862296882804695; w[12] = 1.0140854164204471124; w[13] = -0.22544581011901045383; w[14] = 0.25701475735481036820; w[15] = 0.034174599543251887002; } else if ( order == 17 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.87500000000000000000; x[2] = -0.75000000000000000000; x[3] = -0.62500000000000000000; x[4] = -0.50000000000000000000; x[5] = -0.37500000000000000000; x[6] = -0.25000000000000000000; x[7] = -0.12500000000000000000; x[8] = 0.00000000000000000000; x[9] = 0.12500000000000000000; x[10] = 0.25000000000000000000; x[11] = 0.37500000000000000000; x[12] = 0.50000000000000000000; x[13] = 0.62500000000000000000; x[14] = 0.75000000000000000000; x[15] = 0.87500000000000000000; x[16] = 1.00000000000000000000; w[0] = 0.030797894233299012495; w[1] = 0.26128238288028031086; w[2] = -0.36795289329867605622; w[3] = 1.7037379778090086905; w[4] = -3.9501480717783930427; w[5] = 8.5525299934402953388; w[6] = -13.934614237197880038; w[7] = 19.180342211078732848; w[8] = -20.951950514333334128; w[9] = 19.180342211078732848; w[10] = -13.934614237197880038; w[11] = 8.5525299934402953388; w[12] = -3.9501480717783930427; w[13] = 1.7037379778090086905; w[14] = -0.36795289329867605622; w[15] = 0.26128238288028031086; w[16] = 0.030797894233299012495; } else if ( order == 18 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.88235294117647058824; x[2] = -0.76470588235294117647; x[3] = -0.64705882352941176471; x[4] = -0.52941176470588235294; x[5] = -0.41176470588235294118; x[6] = -0.29411764705882352941; x[7] = -0.17647058823529411765; x[8] = -0.058823529411764705882; x[9] = 0.058823529411764705882; x[10] = 0.17647058823529411765; x[11] = 0.29411764705882352941; x[12] = 0.41176470588235294118; x[13] = 0.52941176470588235294; x[14] = 0.64705882352941176471; x[15] = 0.76470588235294117647; x[16] = 0.88235294117647058824; x[17] = 1.00000000000000000000; w[0] = 0.029364429446790078519; w[1] = 0.23907238516051669677; w[2] = -0.28784319231183443641; w[3] = 1.2897348026109258587; w[4] = -2.5325477495812627261; w[5] = 4.7026959045817496499; w[6] = -5.7913308450170443690; w[7] = 5.3475000248456540826; w[8] = -1.9966457597354948350; w[9] = -1.9966457597354948350; w[10] = 5.3475000248456540826; w[11] = -5.7913308450170443690; w[12] = 4.7026959045817496499; w[13] = -2.5325477495812627261; w[14] = 1.2897348026109258587; w[15] = -0.28784319231183443641; w[16] = 0.23907238516051669677; w[17] = 0.029364429446790078519; } else if ( order == 19 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.88888888888888888889; x[2] = -0.77777777777777777778; x[3] = -0.66666666666666666667; x[4] = -0.55555555555555555556; x[5] = -0.44444444444444444444; x[6] = -0.33333333333333333333; x[7] = -0.22222222222222222222; x[8] = -0.11111111111111111111; x[9] = 0.00000000000000000000; x[10] = 0.11111111111111111111; x[11] = 0.22222222222222222222; x[12] = 0.33333333333333333333; x[13] = 0.44444444444444444444; x[14] = 0.55555555555555555556; x[15] = 0.66666666666666666667; x[16] = 0.77777777777777777778; x[17] = 0.88888888888888888889; x[18] = 1.00000000000000000000; w[0] = 0.026790824664820447344; w[1] = 0.24310820888374278151; w[2] = -0.42247620621346493274; w[3] = 2.0421742376029227612; w[4] = -5.5714791681749728126; w[5] = 13.634004454324976218; w[6] = -26.122288374274995239; w[7] = 41.953237533490708445; w[8] = -55.115367445968607749; w[9] = 60.664591871329740161; w[10] = -55.115367445968607749; w[11] = 41.953237533490708445; w[12] = -26.122288374274995239; w[13] = 13.634004454324976218; w[14] = -5.5714791681749728126; w[15] = 2.0421742376029227612; w[16] = -0.42247620621346493274; w[17] = 0.24310820888374278151; w[18] = 0.026790824664820447344; } else if ( order == 20 ) { /* 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[0] = -1.00000000000000000000; x[1] = -0.89473684210526315789; x[2] = -0.78947368421052631579; x[3] = -0.68421052631578947368; x[4] = -0.57894736842105263158; x[5] = -0.47368421052631578947; x[6] = -0.36842105263157894737; x[7] = -0.26315789473684210526; x[8] = -0.15789473684210526316; x[9] = -0.052631578947368421053; x[10] = 0.052631578947368421053; x[11] = 0.15789473684210526316; x[12] = 0.26315789473684210526; x[13] = 0.36842105263157894737; x[14] = 0.47368421052631578947; x[15] = 0.57894736842105263158; x[16] = 0.68421052631578947368; x[17] = 0.78947368421052631579; x[18] = 0.89473684210526315789; x[19] = 1.00000000000000000000; w[0] = 0.025670822345560078100; w[1] = 0.22448968595251886556; w[2] = -0.34467890099030890987; w[3] = 1.5996974366978074270; w[4] = -3.8466730910952978835; w[5] = 8.3065993344729824120; w[6] = -13.139430424771119113; w[7] = 16.333513604742678295; w[8] = -13.792641220001198577; w[9] = 5.6334527526463774045; w[10] = 5.6334527526463774045; w[11] = -13.792641220001198577; w[12] = 16.333513604742678295; w[13] = -13.139430424771119113; w[14] = 8.3065993344729824120; w[15] = -3.8466730910952978835; w[16] = 1.5996974366978074270; w[17] = -0.34467890099030890987; w[18] = 0.22448968595251886556; w[19] = 0.025670822345560078100; } else if ( order == 21 ) { x[0] = -1.00000000000000000000; x[1] = -0.90000000000000000000; x[2] = -0.80000000000000000000; x[3] = -0.70000000000000000000; x[4] = -0.60000000000000000000; x[5] = -0.50000000000000000000; x[6] = -0.40000000000000000000; x[7] = -0.30000000000000000000; x[8] = -0.20000000000000000000; x[9] = -0.10000000000000000000; x[10] = 0.00000000000000000000; x[11] = 0.10000000000000000000; x[12] = 0.20000000000000000000; x[13] = 0.30000000000000000000; x[14] = 0.40000000000000000000; x[15] = 0.50000000000000000000; x[16] = 0.60000000000000000000; x[17] = 0.70000000000000000000; x[18] = 0.80000000000000000000; x[19] = 0.90000000000000000000; x[20] = 1.00000000000000000000; w[0] = 0.023650546498063206389; w[1] = 0.22827543528921394997; w[2] = -0.47295674102285392846; w[3] = 2.4123737869637513288; w[4] = -7.5420634534306609355; w[5] = 20.673596439879602287; w[6] = -45.417631687959024596; w[7] = 83.656114844387109207; w[8] = -128.15055898030800930; w[9] = 165.59456694494570344; w[10] = -180.01073427048578932; w[11] = 165.59456694494570344; w[12] = -128.15055898030800930; w[13] = 83.656114844387109207; w[14] = -45.417631687959024596; w[15] = 20.673596439879602287; w[16] = -7.5420634534306609355; w[17] = 2.4123737869637513288; w[18] = -0.47295674102285392846; w[19] = 0.22827543528921394997; w[20] = 0.023650546498063206389; } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "NCC_SET - Fatal error!\n" ); fprintf ( stderr, " Illegal value of ORDER = %d\n", order ); fprintf ( stderr, " Legal values are 1 through 21.\n" ); exit ( 1 ); } return; } /******************************************************************************/ double *weights_ls ( int d, double a, double b, int n, double x[] ) /******************************************************************************/ /* Purpose: 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, int D, the desired degree of the polynomial approximant. 0 <= D. Input, double A, B, the left and right endpoints of the interval. Input, int N, the number of data points or abscissas. Input, double X[N], the data points. Output, double WEIGHTS_LS[N], the quadrature weights. */ { double ai; double bi; int i; int j; double *q; double *r; double *s; double *v; double *vt; double *w; /* Integrate the basis functions x^0, x^1, ..., x^d over [a,b]. */ q = ( double * ) malloc ( ( d + 1 ) * sizeof ( double ) ); bi = b; ai = a; for ( i = 0; i < d + 1; i++ ) { q[i] = ( bi - ai ) / ( double ) ( i + 1 ); ai = ai * a; bi = bi * b; } /* Form the N by D+1 Vandermonde matrix V. */ v = ( double * ) malloc ( n * ( d + 1 ) * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { v[i+0*n] = 1.0; for ( j = 1; j < d + 1; j++ ) { v[i+j*n] = v[i+(j-1)*n] * x[i]; } } /* Solve the D+1 by N system V' * S = Q for vector S of length N. */ vt = r8mat_transpose_new ( n, d + 1, v ); s = qr_solve ( d + 1, n, vt, q ); /* Solve the N by D+1 system V * R = S for vector R, of length D+1. */ r = qr_solve ( n, d + 1, v, s ); /* Compute W = R' * V'; */ w = r8mat_mv_new ( n, d + 1, v, r ); /* Free memory. */ free ( q ); free ( r ); free ( s ); free ( v ); free ( vt ); return w; }