Mathematica Asked on January 20, 2021
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.
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 *)
Answered by Michael E2 on January 20, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP