I wanted to use Gauss-Legendre Quadrature to calculate an integral as follows: enter image description here

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.

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

  • thank you very much! I thought it was easy to solve for the roots of polynomial of large degree before. – Quere Sep 14 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 at 12:33
  • Yes!I noticed it.The result was really amazing. – Quere Sep 14 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. is computer-less Sep 24 at 19:45

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

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. is computer-less Sep 24 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 at 17:34
  • Argh, I switched up "comrade" and "colleague" again... I keep switching the general and specific cases. :D – J. M. is computer-less Sep 25 at 17:40

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.