Accuracy and WorkingPrecision

Question

Observe that

(*Gauss Laguerre*)
myglxw[n_Integer] := Block[{
    x,
    xis = Cases[NSolve[LaguerreL[n, x] == 0], _?NumericQ, Infinity]
    },

   {xis, xis/(n + 1)^2/LaguerreL[n + 1, xis]^2}

   ];

myglxw2[n_Integer] := Block[{
    x,
    xis = 
     Cases[NSolve[LaguerreL[n, x] == 0, WorkingPrecision -> 30], _?
       NumericQ, Infinity]
    },

   {xis, xis/(n + 1)^2/LaguerreL[n + 1, xis]^2}

   ];



myglxw2[10] - myglxw[10] // Chop
myglxw2[15] - myglxw[15] // Chop
myglxw2[20] - myglxw[20] // Chop

There differences look "small" (negligible).

myglxw2[50] - myglxw[50] // Chop
myglxw2[100] - myglxw[100] // Chop

Some really big differences.

Question: When to know I "have" to increase WorkingPrecision?

I only notice this because something went seriously wrong at some stage during a long calculation.

As a rule of thumb: Whenever you expect intermediate results near the order of machine precision, you should be cautious. You have that, since there are LaguerreL-intermediates in the order of 10^37, which get squared and divided by.
In any case you might try to stick to Solve here. — Jinxed, Jan 29 '15 at 20:10

score 3 · Accepted Answer · edited Jun 11 '16 at 06:47

If you can tolerate some slowdown, use arbitrary-precision numbers instead of machine-precision numbers. Unlike machine-precision numbers, arbitrary-precision numbers keep track of how error propagates through a calculation. Let's see what happens when we try your NSolve after setting the precision explicitly:

NSolve[SetPrecision[LaguerreL[50, x], $MachinePrecision] == 0, {x}]

(Note that you must use the number $MachinePrecision and not the symbol MachinePrecision.)

The tooltips for the error boxes read No significant digits are available to display due to a total loss of precision. Using SetPrecision[_, 30] instead does not result in any loss of precision.

However, the timing difference is significant:

Timing[Do[NSolve[LaguerreL[50, x] == 0, {x}], {100}];]
(* {0.562500,Null} *)
Timing[Do[
   NSolve[SetPrecision[LaguerreL[50, x], $MachinePrecision] == 
     0, {x}], {100}];]
(* {2.687500,Null} *)

score 2 · Answer 2 · edited Apr 13 '17 at 12:56

As already noted, due to the large range of variation between the nodes and weights of the Gauss-Laguerre rule, one would usually want to use arbitrary precision evaluation for high orders.

In any case, let me present two alternative approaches to generating the nodes and weights for Gauss-Laguerre quadrature. I'll be linking to the papers explaining these methods; read them if you want to understand the interplay between linear algebra and orthogonal polynomials that underlies the following routines.

Here's the Golub-Welsch algorithm:

gaussLaguerreGW[n_Integer?Positive, prec_: MachinePrecision] := 
     With[{od = Range[n - 1]},
          MapAt[Map[First, #]^2 &, 
                Eigensystem[N[SparseArray[{Band[{2, 1}] -> od, 
                                           Band[{1, 1}] -> Range[1, 2 n - 1, 2],
                                           Band[{1, 2}] -> od}], prec]], 2]]

Here's Laurie's algorithm:

gaussLaguerreL[n_Integer?Positive, prec_: MachinePrecision] := 
     With[{sq = Sqrt[Range[n]]},
          {Normal[Diagonal[#2]^2], First[#1]^2} & @@ 
          Take[SingularValueDecomposition[N[SparseArray[{Band[{1, 1}] -> sq,
                                                         Band[{2, 1}] -> Most[sq]}],
                                            prec], Tolerance -> 0], 2]]

(I had previously presented these two approaches in a previous answer, but specialized to the Gauss-Legendre case.)

The two methods are intimately related, in the sense that the eigenvalues of a symmetric positive definite matrix are related to the singular values of that matrix's Cholesky triangle. This works especially well in the Gauss-Laguerre case, since its Jacobi matrix (matrix of recurrence coefficients for the Laguerre polynomials) is symmetric positive definite, so one can easily use either of the approaches of Golub-Welsch or Laurie.

Accuracy and WorkingPrecision

2 Answers2

Linked