Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50?

Question

I wanted to use Gauss-Legendre Quadrature to calculate an integral as follows：

When n=10 and some other number(except odd numbers),the numerical result is the same as theoretical result.

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 -> nodes[[i]]}, {i, 1, Length@nodes}];
Sum[coef[[i]] f[nodes[[i]]], {i, 1, Length@nodes}]
Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -> True]

-2.11450175075
-2.11450175075

But,when n>=50,the code will give wrong numerical result.

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 -> nodes[[i]]}, {i, 1, Length@nodes}];
Sum[coef[[i]] f[nodes[[i]]], {i, 1, Length@nodes}]
Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -> True]

-43.8873164858
-2.11450175075

That is very puzzling.

Related: https://mathematica.stackexchange.com/questions/72807/accuracy-and-workingprecision (Gauss-Laguerre instead of Gauss-Legendre). — Michael E2, Sep 17 '18 at 03:08

Henrik Schumacher · Answer 1 · 2018-09-16T13:51:07.477

6

NSolve 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):

Max[Abs[L[x] /. sol]]

164.177

That seems to be a precision problem (the polynomial hase huge coefficients but the powers of x tend to be very small so that this easily lead to catastrophic cancellation. You will get better results with higher WorkingPrecision:

sol = NSolve[L[x] == 0, x, WorkingPrecision -> 100];
Max[Abs[L[x] /. sol]]

0.*10^-92

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

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 -> True];
a - b

2.6225212*10^-190

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)...

edited Sep 16 '18 at 13:51

answered Sep 14 '18 at 05:56

Henrik Schumacher

106,770
7
179
309

thank you very much! I thought it was easy to solve for the roots of polynomial of large degree before. – Quere Sep 14 '18 at 12:25
@Quere In the meantime, I found out that using higher WorkingPrecision helps to make your method more precise. – Henrik Schumacher Sep 14 '18 at 12:33
Yes!I noticed it.The result was really amazing. – Quere Sep 14 '18 at 12:37
1

Not really worth writing as a separate answer, so: the trick of using even-order Gauss-Legendre quadrature to evaluate principal value integrals with a singularity at the center was published by Piessens, but I seem to recall this being known as folklore. In case the pole is not exactly in the middle of the integration interval, one can either split the integral so that there is a regular part and a part symmetric about the pole, or one can use a substitution that centers the pole, allowing the use of Gaussian quadrature. – J. M.'s missing motivation Sep 24 '18 at 19:45

Stratus · Answer 2 · 2018-10-03T23:20:17.293

If you add WorkingPrecision -> 60 in the NSolve function you get your solution.

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 -> 60];
nodes = (x /. sol);
coef = Table[A[x] /. {x -> nodes[[i]]}, {i, 1, Length@nodes}];
Sum[coef[[i]] f[nodes[[i]]], {i, 1, Length@nodes}]
Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -> True]

-2.11450175075145702914368470979175591804810787513962

-2.1145

Or you can use the built in functions i mma

<< NumericalDifferentialEquationAnalysis`
n = 50;
{pts, w} = Transpose[GaussianQuadratureWeights[n, -1, 1]];
Integrate[E^(-x)/x, {x, -1, 1.}, PrincipalValue -> True]
Sum[w[[i]] f[pts[[i]]], {i, 1, Length@pts}]

-2.1145

-2.1145

score 4 · Answer 3 · answered Sep 15 '18 at 03:23

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."

(*
 * 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}] -> -b/a, Band[{2, 1}] -> beta, 
       Band[{1, 2}] -> 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}
*)

The integration weights can be computed from these nodes as in the OP.

As a terminological note: the tridiagonal matrix with the orthogonal polynomial's recurrence coefficients is called a Jacobi matrix. (The comrade matrix, as you might recall, is a Jacobi matrix for the Chebyshev polynomials, with a correction term representing the Chebyshev series coefficients.) I talked about this here. — J. M.'s missing motivation, Sep 24 '18 at 11:27
@J.M. Barnett (1975) seems to introduce the name "comrade matrix" of a polynomial with respect to an arbitrary family of polynomials that satisfies a three-term recurrence and presents it as a generalization of the Chebyshev colleague matrix of Good (1961). But, you know, I'm not that confident about current terminology in this area. -- Yes, I found the link to that Q&A shortly after posting an answer and put it under the question above. — Michael E2, Sep 25 '18 at 17:34
Argh, I switched up "comrade" and "colleague" again... I keep switching the general and specific cases. :D — J. M.'s missing motivation, Sep 25 '18 at 17:40

Why Gauss-Legendre Quadrature should keep the number of integral points less than about 50?

3 Answers3

Linked