Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50? - Mathematica Stack Exchange most recent 30 from mathematica.stackexchange.com 2019-08-21T23:47:41Z https://mathematica.stackexchange.com/feeds/question/181868 http://www.creativecommons.org/licenses/by-sa/3.0/rdf https://mathematica.stackexchange.com/q/181868 5 Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50? Quere https://mathematica.stackexchange.com/users/51355 2018-09-14T05:29:51Z 2018-10-03T23:20:17Z <p>I wanted to use Gauss-Legendre Quadrature to calculate an integral as follows： <a href="https://i.stack.imgur.com/uvm23.png" rel="noreferrer"><img src="https://i.stack.imgur.com/uvm23.png" alt="enter image description here"></a></p> <p>When n=10 and some other number(except odd numbers),the numerical result is the same as theoretical result.</p> <pre><code>Clear["Global`*"]; n = 10; L[x_] := D[(x^2 - 1)^n, {x, n}]/(n!*2^n); A[x_] := 2/((1 - x^2) (D[L[x], x])^2); f[x_] := E^(-x)/x; sol = NSolve[L[x] == 0., x]; nodes = (x /. sol); coef = Table[A[x] /. {x -&gt; nodes[[i]]}, {i, 1, Length@nodes}]; Sum[coef[[i]] f[nodes[[i]]], {i, 1, Length@nodes}] Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -&gt; True] -2.11450175075 -2.11450175075 </code></pre> <p>But,when n>=50,the code will give wrong numerical result.</p> <pre><code>Clear["Global`*"]; n = 50; L[x_] := D[(x^2 - 1)^n, {x, n}]/(n!*2^n); A[x_] := 2/((1 - x^2) (D[L[x], x])^2); f[x_] := E^(-x)/x; sol = NSolve[L[x] == 0., x]; nodes = (x /. sol); coef = Table[A[x] /. {x -&gt; nodes[[i]]}, {i, 1, Length@nodes}]; Sum[coef[[i]] f[nodes[[i]]], {i, 1, Length@nodes}] Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -&gt; True] -43.8873164858 -2.11450175075 </code></pre> <p>That is very puzzling.</p> https://mathematica.stackexchange.com/questions/181868/-/181869#181869 6 Answer by Henrik Schumacher for Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50? Henrik Schumacher https://mathematica.stackexchange.com/users/38178 2018-09-14T05:56:52Z 2018-09-16T13:51:07Z <p><code>NSolve</code> has severe issues to solve for the roots (which is not surprising since finding the roots of a polynomial of degree 50 is a nontrivial task):</p> <pre><code>Max[Abs[L[x] /. sol]] </code></pre> <blockquote> <p>164.177</p> </blockquote> <p>That seems to be a precision problem (the polynomial hase huge coefficients but the powers of <code>x</code> tend to be very small so that this easily lead to catastrophic cancellation. You will get better results with higher <code>WorkingPrecision</code>:</p> <pre><code>sol = NSolve[L[x] == 0, x, WorkingPrecision -&gt; 100]; Max[Abs[L[x] /. sol]] </code></pre> <blockquote> <p>0.*10^-92</p> </blockquote> <p>So the real art here is to determine the roots in a numerically stable way. Actually, Gauß-Legendre quadrature rule is already built into the system. For example, see</p> <pre><code>n = 50; prec = 200; {nodes, coef, bla} = NIntegrate`GaussRuleData[n, prec]; (* compute all quadrature points by convex combinations*) nodes = (-1) * (1 - nodes) + nodes * 1; coef = 2 coef; a = coef.f[nodes]; b = Integrate[f[x], {x, -1, 1}, PrincipalValue -&gt; True]; a - b </code></pre> <blockquote> <p>2.6225212*10^-190</p> </blockquote> <p>PS.: I was quite astonished that Gauß quadrature works so well with this singular integral. Then it came back to my mind that it is constructed such that it neglects contributions of all odd functions (because they would integrate to 0 anyways)...</p> https://mathematica.stackexchange.com/questions/181868/-/181905#181905 5 Answer by Diogo for Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50? Diogo https://mathematica.stackexchange.com/users/36260 2018-09-14T15:08:06Z 2018-10-03T23:20:17Z <p>If you add <code>WorkingPrecision -&gt; 60</code> in the NSolve function you get your solution. </p> <pre><code>n = 50; A[x_] := 2/((1 - x^2) (D[LegendreP[n, x], x])^2); f[x_] := E^(-x)/x; sol = NSolve[LegendreP[n, x] == 0, x, WorkingPrecision -&gt; 60]; nodes = (x /. sol); coef = Table[A[x] /. {x -&gt; nodes[[i]]}, {i, 1, Length@nodes}]; Sum[coef[[i]] f[nodes[[i]]], {i, 1, Length@nodes}] Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -&gt; True] -2.11450175075145702914368470979175591804810787513962 -2.1145 </code></pre> <p>Or you can use the built in functions i mma</p> <pre><code>&lt;&lt; NumericalDifferentialEquationAnalysis` n = 50; {pts, w} = Transpose[GaussianQuadratureWeights[n, -1, 1]]; Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -&gt; True] Sum[w[[i]] f[pts[[i]]], {i, 1, Length@pts}] -2.1145 -2.1145 </code></pre> https://mathematica.stackexchange.com/questions/181868/-/181927#181927 3 Answer by Michael E2 for Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50? Michael E2 https://mathematica.stackexchange.com/users/4999 2018-09-15T03:23:59Z 2018-09-15T03:23:59Z <p>A more stable way of computing the roots of a polynomial family that satisfies a three-term recurrence is the method of Golub and Welsch (1969), which computes the eigenvalues of a matrix based on the recurrence, sometimes called the "comrade matrix."</p> <pre><code>(* * n-point Gauss quadrature (Golub-Welsch 1969) *) (* p[n][x] == (a[n]x+b[n])p[n-1][x] - c[n]p[n-2][x] *) ClearAll[comradeMatrix]; comradeMatrix[{a_, b_, c_}, {n_, n0_Integer}, prec_: Infinity] := Block[{n = Range@n0}, With[{beta = Sqrt[Rest@c/(Most@a*Rest@a)]}, N[SparseArray[{Band[{1, 1}] -&gt; -b/a, Band[{2, 1}] -&gt; beta, Band[{1, 2}] -&gt; beta}, {n0, n0}], prec] ]]; comradeMatrix["Legendre", n0_Integer, prec_: Infinity] := Module[{n}, comradeMatrix[{2 n - 1, 0, n - 1}/n, {n, n0}, prec]]; Sort@Eigenvalues@N@comradeMatrix["Legendre", 50] (* {-0.998866, -0.994032, -0.985354, -0.972864, -0.956611, -0.936657, \ -0.913079, -0.885968, -0.85543, -0.821582, -0.784556, -0.744494, \ -0.701552, -0.655896, -0.607703, -0.557158, -0.504458, -0.449806, \ -0.393414, -0.3355, -0.276288, -0.216007, -0.154891, -0.0931747, \ -0.0310983, 0.0310983, 0.0931747, 0.154891, 0.216007, 0.276288, \ 0.3355, 0.393414, 0.449806, 0.504458, 0.557158, 0.607703, 0.655896, \ 0.701552, 0.744494, 0.784556, 0.821582, 0.85543, 0.885968, 0.913079, \ 0.936657, 0.956611, 0.972864, 0.985354, 0.994032, 0.998866} *) </code></pre> <p>The integration weights can be computed from these nodes as in the OP.</p>