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

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.

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.