NIntegrate, how exact is the answer?

Question

Since the precision is very important in my case, I am curious how precise is the MMA result for this integral:

t = 0.2; w = 1; a = 4/w^2;

XX[lf_, mf_, l_, mg_, p_] := Abs[YY[lf, mf, l, mg, p]]^2;
YY[lf_, mf_, l_, mg_, p_] := WignerD[{lf, mf, l}, 
t] Exp[-((a r^2)/4)] Sum[((Abs[mg - l] + p)!/((p - j)! (Abs[mg - l] + j)! j!)) (Sqrt[2]/w)^(2 j + mg - l) (2/a)^(2 j + mg - l + 1) (Gamma[.5 (Abs[mg - mf] + mg - l + 2 + 2 j)]/(r Gamma[Abs[mg - mf] + 1])) a^(0.5 (2 j + mg - l + 1)) ((a r^2)/2)^(1/2 Abs[mg - mf] + .5) Hypergeometric1F1[(Abs[mg - mf] - 2 j - mg + l)/2, Abs[mg - mf] + 1, (a r^2)/4], {j, 0, p}];

MM[lf_, l_, mg_, p_] := Sum[XX[lf, ii, l, mg, p], {ii, -lf, lf, 1}];

NIntegrate[r M[2, 1, 3, 0], {r, 0, \[Infinity]}, MaxPoints -> 500000,PrecisionGoal -> 10, AccuracyGoal -> 5, MaxRecursion -> 50,WorkingPrecision -> 10]

MMA returns 4.159915211

But the result is working precision dependent (and not only):

If I put working precision 17 it says that the precision of the integrand is less than working precision. What does it mean? Can I still trust this result?

I wonder if there are some tools that could give me the error bar on how accurate this integration is since it seems like depending on the evaluation parameters I get different results. How to figure out the optimal set of parameters for this integral?

Answer 1

The precision of all numbers must be at least WorkingPrecision. Use exact numbers when possible to support subsequent use of any WorkingPrecision

t = 1/5; (* use exact number rather than 0.2 *)

w = 1; a = 4/w^2;

XX[lf_, mf_, l_, mg_, p_] := Abs[YY[lf, mf, l, mg, p]]^2;

(* Replace all 0.5 by 1/2 *)

YY[lf_, mf_, l_, mg_, p_] := 
  WignerD[{lf, mf, l}, 
    t] Exp[-((a r^2)/
       4)] Sum[((Abs[mg - l] + 
          p)!/((p - j)! (Abs[mg - l] + j)! j!)) (Sqrt[2]/w)^(2 j + 
        mg - l) (2/a)^(2 j + mg - l + 
        1) (Gamma[(Abs[mg - mf] + mg - l + 2 + 2 j)/
         2]/(r Gamma[Abs[mg - mf] + 1])) a^((2 j + mg - l + 1)/
        2) ((a r^2)/2)^(1/2 Abs[mg - mf] + 
        1/2) Hypergeometric1F1[(Abs[mg - mf] - 2 j - mg + l)/2, 
      Abs[mg - mf] + 1, (a r^2)/4], {j, 0, p}];

MM[lf_, l_, mg_, p_] := Sum[XX[lf, ii, l, mg, p], {ii, -lf, lf, 1}];

m = MM[2, 1, 3, 0] // Simplify[#, r >= 0] &; 
(* simplify to reduce the number of operations used. 
   Each operation tends to reduce precision *)

Grid[data =
  {#, NIntegrate[r m, {r, 0, ∞},
      MaxPoints -> 500,
      MaxRecursion -> 5,
      WorkingPrecision -> #]} & /@ Range[10, 20],
 Alignment -> "."]