0
$\begingroup$

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

enter image description here

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?

$\endgroup$
  • 4
    $\begingroup$ (1) Please include the message name and the tag error when asking about error messages. (2) Try t = 2/10 instead of the MachinePrecision number 0.2. Read some of the tutorials on precision for more explanation. (3) See (14500), (75426) for getting the error estimate. $\endgroup$ – Michael E2 Dec 7 '16 at 1:35
4
$\begingroup$

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

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ how do I choose the parameters like MaxPoints and MaxRecursion. For example I had MaxPoints ->500000, you took 500. Why do you think 500 is enough? $\endgroup$ – MsTais Dec 7 '16 at 3:19
  • 2
    $\begingroup$ Start with a "reasonable" number (500 in this case) from a performance standpoint and get a result. Double the value (1000) and compare the results. Since there is no change, then 500 was good enough. $\endgroup$ – Bob Hanlon Dec 7 '16 at 3:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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