Message Boards Message Boards

GROUPS:

[✓] Improve the results of NIntegrate - methods?

Posted 10 months ago
1145 Views
|
3 Replies
|
2 Total Likes
|

The following result of NIntegrate is almost consistent (second decimal place) across different methods (e.g. GlobalAdaptive versus DuffyCoordinates), and different values of PrecisionGoal and AccuracyGoal, but the Mathematica warnings will not go away. There is no singularity as far as I can see (the integral is regularised).

Limit[Block[{p4 = p1, p3 = p1, \[Omega] = 0}, NIntegrate[(p3*p4*(p1*p3 + p1*p4 + p3*p4)*Log[((p1 + p3)^2 + \[Omega]^2)/((-p1 + p3)^2 + \[Omega]^2)]*Log[((p1 + p4)^2 + \[Omega]^2)/((-p1 + p4)^2 + \[Omega]^2)]*Sin[p1])/(p1*(1 + p3^2)^2*(1 + p4^2)^2), {p1, 0, 10}, {p3, 0, 10}, 
    {p4, 0, 10}, PrecisionGoal -> 30, WorkingPrecision -> 40]], \[Omega] -> 1.*^-10]

enter image description here

    11.6717784706425703112853459084272623029371526368002822614702`40.

Any suggestions on how to not get the warnings, or otherwise believe the result - since they are very similar with different methods?

3 Replies

3-digit precision:

PrintTemporary@Dynamic@Clock[Infinity]; (* timer *)
Block[{(*p4=p1,p3=p1,*) \[Omega] = 0},    (* <-- hmm? *)
  (* symmetric w.r.t. p3 <-> p4, so divide region in half and multiply by 2 *)
  2 NIntegrate[(p3*p4*(p1*p3 + p1*p4 + p3*p4)*
       Log[((p1 + p3)^2 + \[Omega]^2)/((-p1 + p3)^2 + \[Omega]^2)]*
       Log[((p1 + p4)^2 + \[Omega]^2)/((-p1 + p4)^2 + \[Omega]^2)]*
       Sin[p1])/(p1*(1 + p3^2)^2*(1 + p4^2)^2),
    {p3, 0, 10}, {p4, 0, p3}, {p1, 0, p4, p3, 10},  (* near singularities at p1 = p3, p4; switching order helps *)
    PrecisionGoal -> 3,  (* 3-digit accuracy sought *)
    Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0, "Points" -> 21}, (* somewhat inefficient but reliable *)
    MinRecursion -> 1, MaxRecursion -> 2]  (* probably unnecessary -- just to make it quit soon *)
    ] // AbsoluteTiming

(*  {34.883, 11.6841}  *)

Ditto for 6-digit precision:

PrintTemporary@Dynamic@Clock[Infinity];
Block[{(*p4=p1,p3=p1,*)\[Omega] = 0}, 
  2 NIntegrate[(p3*p4*(p1*p3 + p1*p4 + p3*p4)*
       Log[((p1 + p3)^2 + \[Omega]^2)/((-p1 + p3)^2 + \[Omega]^2)]*
       Log[((p1 + p4)^2 + \[Omega]^2)/((-p1 + p4)^2 + \[Omega]^2)]*
       Sin[p1])/(p1*(1 + p3^2)^2*(1 + p4^2)^2),
    {p3, 0, 10}, {p4, 0, p3}, {p1, 0, p4, p3, 10},
    PrecisionGoal -> 6, 
    Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0, "Points" -> 21},
    MinRecursion -> 1, MaxRecursion -> 2]
    ] // AbsoluteTiming

(*  {333.226, 11.6855}  *)

The result seems stable. Now if you really want 30 digits of precision (PrecisionGoal -> 30), you're going to have to wait a while and increase "MaxErrorIncreases".

Posted 10 months ago

Thank you, Michael!

Are the results coming to be different for the two integrals, for the same limits?

Hi Arny,

I think you're asking why the lower-order digits are different in the two results I obtained.

The setting PrecisionGoal -> 3 tells NIntegrate to refine the integration until the error estimate indicates that the first three digits are correct (roughly speaking). Similarly, the setting PrecisionGoal -> 6 tells it to go until the first six digits are accurate. So one would expect the two answers to agree to three digits (or more). They agree to four. That suggests that the integration approach (a product of Gauss-Kronrod rules) is converging on the true answer.

In your original question, the error message indicates that the first two digits are supposed to be accurate (the message reports both the value and the error estimate). Those two digits, 11. agree with the first two digits of my answer; so all three are consistent within the error estimate of each result. In fact the third digit of your result 11.6 seems correct; this is consistent with the error estimate (and typical), because the error estimate is supposed to be an upper bound.

Here is another method with PrecisionGoal -> 6. Raising MinRecursion (from the default 0) sometimes gets rid of a slow convergence warning, but that warning can usually be ignored provided there are no other errors. The "MultidimensionalRule" is probably what your integral used, but the default number of generators is 5.

PrintTemporary@Dynamic@Clock[Infinity];
Block[{(*p4=p1,p3=p1,*)\[Omega] = 0}, 
  2 NIntegrate[(p3*p4*(p1*p3 + p1*p4 + p3*p4)*
       Log[((p1 + p3)^2 + \[Omega]^2)/((-p1 + p3)^2 + \[Omega]^2)]*
       Log[((p1 + p4)^2 + \[Omega]^2)/((-p1 + p4)^2 + \[Omega]^2)]*
       Sin[p1])/(p1*(1 + p3^2)^2*(1 + p4^2)^2),
    {p3, 0, 10}, {p4, 0, p3}, {p1, 0, p4, p3, 10},
    PrecisionGoal -> 6, 
    Method -> {"MultidimensionalRule", "SymbolicProcessing" -> 0,  "Generators" -> 9},
    MinRecursion -> 2]] // AbsoluteTiming

(*  {18.832607`, 11.68553522194798`}  *)

Here is the full result of my Gauss-Kronrod six-digit result:

{333.225781`, 11.685534968838665`}

The setting PrecisionGoal -> 6 means that normally they should agree to at least six digits. In fact, they agree to eight (if rounded). It looks to me that the integral is well-behaved.

In your original integral, the slow convergence probably has something to do with the "MaxErrorIncreases" problem. The main issue is trying to compute a thirty-digit answer (PrecisionGoal -> 30). I suspect raising MinRecursion would help. You would have to raise it to a pretty large amount, I would think. One could also try raising "MaxErrorIncreases", or a combination of both. You can see that your integral still wasn't very close (error approx. 0.188) when it hit the internal limit.

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract