0
$\begingroup$

The following call to NIntegrate throws the error:

Throw::sysexc: Uncaught SystemException returned to top level. Can be caught with Catch[[Ellipsis], _SystemException].

and returns SystemException["MemoryAllocationFailure"] as output.

NIntegrate[
 PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Sec[x]^2 /. {k -> 
    1/100 Tan[x], T -> 1/1000, a -> 27/100}, {x, 0, \[Pi]/2}, 
 PrecisionGoal -> 12]

Meanwhile, adding additionally WorkingPrecision->20 causes the evaluation to last too long and results in a kernel crash.

Any ideas what the root problem is and how to fix this? I really would like to evaluate this integral to many digits of precision (the issue goes away at lower precision), and I'm looking for a reliable way to calculate many integrals of this type with slightly different parameter values.

$\endgroup$
4
  • $\begingroup$ Try using WorkingPrecision->200 instead. $\endgroup$
    – Carl Woll
    Oct 30 '20 at 15:09
  • $\begingroup$ NIntegrate seems not to be able to calculate the integral with exact numbers. But it gives an answer if you change some constants to machine numbers. $\endgroup$ Oct 30 '20 at 15:10
  • $\begingroup$ Simplify the integrand: NIntegrate[PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Sec[x]^2 /. {k -> 1/100 Tan[x], T -> 1/1000, a -> 27/100} // Simplify, {x, 0, \[Pi]/2}, WorkingPrecision -> 20] $\endgroup$
    – Bob Hanlon
    Oct 30 '20 at 16:43
  • $\begingroup$ Thanks for the suggestions! They work. @CarlWoll any idea why increasing WorkingPrecision avoids issues involved with memory overload? I would have thought that would make things worse? $\endgroup$
    – WillG
    Oct 30 '20 at 18:15
2
$\begingroup$

I don't know what might be needed in other integrals, but for this one, MinRecursion -> 3 seems to work at various precision goals:

Table[
 NIntegrate[
  PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Sec[x]^2 /. {k -> 
     1/100 Tan[x], T -> 1/1000, a -> 27/100}, {x, 0, π/2}, 
  PrecisionGoal -> pg, MinRecursion -> 3, 
  Method -> {"GlobalAdaptive", 
    Method -> {"GaussKronrodRule", "Points" -> 21}, 
    "SymbolicProcessing" -> None}, WorkingPrecision -> 2 pg],
 {pg, {12, 16, 30}}]
(*
{-1.66625181179756853846696*10^-8, \
-1.6662518117975685384669577064196*10^-8, \
-1.66625181179756853846695769835962377225227776502356624286066*10^-8}
*)

It works at machine precision, too:

NIntegrate[
 PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Sec[x]^2 /. {k -> 
    1/100 Tan[x], T -> 1/1000, a -> 27/100}, {x, 0, π/2}, 
 PrecisionGoal -> 12, MinRecursion -> 3]

Local adaptive integration works, too:

NIntegrate[
 PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Sec[x]^2 /. {k -> 
    1/100 Tan[x], T -> 1/1000, a -> 27/100}, {x, 0, π/2}, 
 PrecisionGoal -> 12, Method -> "LocalAdaptive"]

Update: Yet another fix (suggests bad numerics near Pi/2):

NIntegrate[
 PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Sec[x]^2 /. 
  {k -> 1/100 Tan[x], T -> 1/1000, a -> 27/100},
 {x, 0, π/2 - 0.01, π/2}, PrecisionGoal -> 12]

Analysis and discussion

The first thing to do when something goes wrong is to plot the integrand (note that Dt[100 k] == Sec[x]^2 Dt[x] -- makes me wonder about the 100):

plot = Plot[PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Dt[100 k] /.
     {k -> 1/100 Tan[x], T -> 1/1000, a -> 27/100} /. Dt[x] -> 1 // 
  Evaluate, {x, 0, \[Pi]/2}]

Shortly after x == 1.5 the curve flattens dramatically, which caught my attention. In fact, there's some sort of underflow:

Cases[plot, Line[p_] :> Drop[p[[All, 2]], {10, -13}], Infinity]
(*
{{-3.02837*10^-9, -3.02854*10^-9, -3.02907*10^-9, -3.02925*10^-9, 
  -3.02946*10^-9, -3.02969*10^-9, -3.02994*10^-9, -3.03022*10^-9, 
  -3.03051*10^-9,
  ...,
  -6.103*10^-125, -2.84532*10^-138, -6.38854*10^-154, -1.20053*10^-194,
   0., 0., 0., 0., 0., 0., 0., 0.}}
*)

The OP's issue is certainly connected to this. Another tool to use is NIntegrateSamplingPoints. We'll reduce the AccuracyGoal in the OP's integral to less than 16, so that the integration will terminate without error.

Needs["Integration`NIntegrateUtilities`"]
NIntegrateSamplingPoints[
 NIntegrate[
  PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Dt[100 k] /.
    {k -> 1/100 Tan[x], T -> 1/1000, a -> 27/100} /. Dt[x] -> 1,
  {x, 0, \[Pi]/2}, PrecisionGoal -> 12, AccuracyGoal -> 15]
 ]

One can see the sampling clustering around x == 1.5 for the last hundred sampling points. Something seems to happen to the sampling when the sampling gets above 100. Below it, we see the typical Gauss-Kronrod sampling patterns; also at the top. Something special seems to be happening around 1.5, but I'm not sure what that is. My guess is that the extra sampling gets out of hand in the OP's integral (in which the default AccuracyGoal is Infinity). Note that we cannot use AccuracyGoal -> 15 to achieve the OP's aims, since an accuracy of 15 places past the decimal point in a result around 10^-8 is only around 7 digits of relative precision.

NIntegrateSamplingPoints also reveals what happens with WorkingPrecision -> 200 (@CarlWoll's suggestion). The normal default method uses Gauss-Kronrod with 5 Gauss points and 6 Kronrod points (it's always n and n+1). But for high precision integrations, this is increased and depends on the precision. The increase can be seen in the sampling plot:

Needs["Integration`NIntegrateUtilities`"]
NIntegrateSamplingPoints[
 NIntegrate[
  PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Dt[100 k] /.
    {k -> 1/100 Tan[x], T -> 1/1000, a -> 27/100} /. Dt[x] -> 1,
  {x, 0, \[Pi]/2}, PrecisionGoal -> 12, WorkingPrecision -> 200]
 ]

There are five GK samplings, and we can deduce that the number of Gauss points is 66:

Count[%, _Point, Infinity]/5.
(% - 1)/2
(*
  133.
  66.
*)

This suggests that we could use the GK method with "Points" -> 66, which we can even at machine precision:

NIntegrate[
 PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Dt[100 k] /.
   {k -> 1/100 Tan[x], T -> 1/1000, a -> 27/100} /. Dt[x] -> 1,
 {x, 0, \[Pi]/2}, PrecisionGoal -> 12, 
 Method -> {"GaussKronrodRule", "Points" -> 66}]

(*  -1.66625*10^-8  *)

A better substitution?

Finally, the issue is introduce by the tangent transformation k -> 1/100 Tan[x]. We could look for another one that doesn't have the flattening out around 1.5. A logarithmic one seems likely to be what we want.

Plot[PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Dt[100 k] /.
    {k -> 1/100 Log[1 + x/(1 - x)^2],   (* new transformation *)
     T -> 1/1000, a -> 27/100} /. Dt[x] -> 1 // Evaluate,
 {x, 0, 1}]      (* different interval for new transformation *)

This produces an integral without issues:

NIntegrate[
 PolyLog[-3, -E^(-(Sqrt[a (k^2 + 1)] - 1/2)/T)] Dt[100 k] /.
     {k -> 1/100 Log[1 + x/(1 - x)^2],
      T -> 1/1000, a -> 27/100} /. Dt[x] -> 1 // Evaluate,
 {x, 0, 1}, PrecisionGoal -> 12]

(*  -1.66625*10^-8  *)
$\endgroup$
2
  • $\begingroup$ Thanks for these strategies. I'm going to try implementing them to see if they work on the numerous generic integrals of similar form that I am dealing with. $\endgroup$
    – WillG
    Oct 30 '20 at 18:17
  • $\begingroup$ @WillG I added some more alternatives. $\endgroup$
    – Michael E2
    Oct 30 '20 at 22:27

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.