# Why does precision decrease after increasing PrecisionGoal and WorkingPrecision on this integral?

I am trying to calculate a messy double-integral to many digits of precision. As shown below, the precision seems to increase reliably as WorkingPrecision goes from 16 to 19. Then starting at 20, the result changes in the 5th decimal place. I know that PrecisionGoal automatically adjusts to be 10 less than WorkingPrecision if not set explicitly, which makes the behavior shown below even more surprising.

example[pg_, wp_, mei_: 2000, mr_: Automatic] :=
NIntegrate[
x^(-3/2) 1/(E^(x 6000) - 1) Exp[-x (Sqrt[1 + y^2] - 300)^2],
{x, 1/600^2, 2/600^2, 10/600^2, 100/600^2,
1000/600^2, \[Infinity]}, {y, 0,
Max[0, Abs[Sqrt[300^2 - 1]] - 5/x], Abs[Sqrt[300^2 - 1]],
Abs[Sqrt[300^2 - 1]] + 5/x, \[Infinity]},
PrecisionGoal -> pg, WorkingPrecision -> wp,
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> mei},
MaxRecursion -> mr]


The messy region of integration helps NIntegrate concentrate on sharp peaks in the integrand, which I found to be necessary when testing this and similar integrals.

example[Automatic, 16] // InputForm // Timing

{0.652093, 5.953074338937872380734835516.*^6}

example[Automatic, 17] // InputForm // Timing

{1.45922, 5.9530773637614705765128073465798783694943772717.*^6}

example[Automatic, 18] // InputForm // Timing

{2.85708, 5.953077381545432910401927528344222872977559818.*^6}

example[Automatic, 19] // InputForm // Timing

{6.62376, 5.9530773819531482121505108553712972664073315319.*^6}

example[Automatic, 20] // InputForm // Timing

{16.2174, 5.9532898242902800360842281566447028066369685920.*^6}

The previous line shows that the something breaks at WorkingPrecision->20 because the 5th digit suddenly changes.

example[10, 30] // InputForm // Timing

{15.5158, 5.9532898242902800360842281566447028066544613230.*^6}

The above line shows that it is really increasing PrecisionGoal to 10 that breaks something, and that increasing WorkingPrecision at this value of PrecisionGoal doesn't help.

I would like to be able to reliably compute this integral with up to 20 digits of precision. Any ideas what's causing this problem and how to fix it? I tried changing with the settings MaxErrorIncreases and MaxRecursion to no avail.

• @MichaelE2 in that line, the previously consistent 5.9530... becomes 5.9532..., a sudden change in the 5th digit after increasing from 9-digit precision to 10-digit precision. May 28 '20 at 23:24
• Yes, sorry I added the last line after writing that part. The last line was intended to show that increasing WorkingPrecision while holding PrecisionGoal fixed doesn't solve the problem. I'll update the post. May 28 '20 at 23:34
• Add MinRecursion -> 3 and you'll get a bigger shock. May 28 '20 at 23:42

If we take a look at the integrand (plotting with a change of variables to tangent, so that we can see the integrand at infinity), we notice that the support is mainly concentrated in a small section of a corner. A more judicious change of variables (on the left) scales the spike to cover most of the domain. This suggests the second change of variables might help the integration.

Plot3D[x^(-3/2) 1/(E^(x 6000) - 1) Exp[-x (Sqrt[1 + y^2] - 300)^2] Dt@
x*Dt@y /.
{x -> Tan[s], y -> Tan[Pi/2 - t]} /. _Dt -> 1 //
Evaluate,
{s, 1/600^2 // ArcTan[#] &, Pi/2}, {t, 0, Pi/2},
AxesLabel -> {x, y}, PlotRange -> All, MaxRecursion -> 3,
Ticks -> {Table[{ArcTan[x], x}, {x, {0, 0.5, 1, 2, 10, Infinity}}],
Table[{Pi/2 - ArcTan[y], y}, {y, {0, 0.5, 1, 2, 10, Infinity}}],
Automatic}]

Plot3D[x^(-3/2) 1/(E^(x 6000) - 1) Exp[-x (Sqrt[1 + y^2] - 300)^2] Dt@
x*Dt@y /.
{x -> Tan[s]/600^2,
y -> Tan[Pi/2 - t] 200} /. _Dt -> 1 // Evaluate,
{s, 1/600^2 // ArcTan[600^2 #] &, Pi/2}, {t, 0, Pi/2},
AxesLabel -> {x, y}, PlotRange -> All, MaxRecursion -> 3,
Ticks -> {Table[{ArcTan[600^2 x],
x}, {x, {0, 3.*^-6, 5.*^-6, 7.*^-6, 1.*^-5, 2.*^-5, Infinity}}],
Table[{Pi/2 - ArcTan[y/200],
y}, {y, {0, 50, 100, 200, 500, 1000, Infinity}}], Automatic}]


Here are two helper functions, which are slight modifications of the OP's example[] function. The first does the integral like the OP's, but with the ability to specify the integration rule via the suboption Method as well as other options.

(* Like the OP's original example[] but with Method and options *)
example2[pg_, wp_, mei_ : 2000, mr_ : Automatic, meth_ : Automatic,
opts : OptionsPattern[NIntegrate]] :=
NIntegrate[
x^(-3/2) 1/(E^(x 6000) - 1) Exp[-x (Sqrt[1 + y^2] - 300)^2],
{x, 1/600^2, 2/600^2, 10/600^2, 100/600^2,
1000/600^2, ∞},
{y, 0, Max[0, Abs[Sqrt[300^2 - 1]] - 5/x], Abs[Sqrt[300^2 - 1]],
Abs[Sqrt[300^2 - 1]] + 5/x, ∞},
PrecisionGoal -> pg, WorkingPrecision -> wp,
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> mei,
Method -> meth}, MaxRecursion -> mr, opts]

(* With tangent substitutions for x and y *)
example3[pg_, wp_, mei_ : 2000, mr_ : Automatic, meth_ : Automatic,
opts : OptionsPattern[NIntegrate]] :=
NIntegrate[
x^(-3/2) 1/(E^(x 6000) - 1) Exp[-x (Sqrt[1 + y^2] - 300)^2] Dt@x * Dt@y /.
{x -> Tan[s]/600^2, y -> 200 Tan[Pi/2 - t]} /.
_Dt ->  1 // Evaluate,
{t, Pi/2, 0},
{s, 1/600^2 // ArcTan[600^2 #] &, Pi/2},
PrecisionGoal -> pg, WorkingPrecision -> wp,
Method -> {"GlobalAdaptive", "MaxErrorIncreases" -> mei,
Method -> meth}, MaxRecursion -> mr, opts]


Neither result in the OP seems accurate. The true value appears to be close to $$1.5 \times 10^7$$. The OP tried (in effect) to address this problem by a manual subdivision of the interval, but more than that is needed. That can be obtained by increasing the order of the integration rule or the minimum number of subdivisions.

example2[Automatic, MachinePrecision, 2000,
20, {"GaussKronrod", "Points" -> 9}, MinRecursion -> 0] //
InputForm // AbsoluteTiming
(*  {0.089678, 8.268709756340054*^6}  *)

example2[Automatic, MachinePrecision, 2000,
20, {"GaussKronrod", "Points" -> 9}, MinRecursion -> 1] //
InputForm // AbsoluteTiming
(*  {0.161061, 1.5367709045839794*^7}  *)

example2[Automatic, MachinePrecision, 2000, 20, Automatic,
MinRecursion -> 5] // InputForm // AbsoluteTiming
(*  {1.36757, 1.5367709116318425*^7}  *)


With the tangent substitutions, no special subdivision is needed for NIntegrate to hone in on the correct value. The cartesian product of Gauss-Kronrod rules seems to consistently outperform the multidimensional rule (which is the usual Automatic rule chosen for a multiple integral).

example3[Automatic, MachinePrecision, 2000, 20, "GaussKronrod"] //
InputForm // AbsoluteTiming
(*  {0.026507, 1.5367708664590633*^7}  *)

(res20 = example3[20, 40, 2000, 20, "GaussKronrod"]) //
InputForm // AbsoluteTiming
(*
{9.09103,
1.536770897940587986735671340147671576809806580607507...40.*^7}
*)

(res25 = example3[25, 50, 2000, 20, "GaussKronrod"]) //
InputForm // AbsoluteTiming
(*
{13.7239,
1.536770897940587986735671340149774119056152243104002...50.*^7}
*)

(res30 = example3[30, 60, 2000, 20, "GaussKronrod"]) //
InputForm // AbsoluteTiming
(*
{24.4783,
1.536770897940587986735671340149774119156342012210349...60.*^7}
*)


Check the relative error:

({res20, res25} - res30)/res30
(*  {-1.368156014*10^-30, -6.51949937630*10^-38}  *)


The integral seems to be converging and have at least 20, 25, 30 digits of accuracy, respectively.

• Thank you for this unexpectedly thorough answer! Upon experimenting with this answer, I realized I made the dumb mistake of not subdividing my y region correctly: the two places where I had a 5/x should instead be Sqrt[5/x]. This partly resolves my problem, but your solution is still more accurate and orders of magnitude faster, so very much worthwhile. May 30 '20 at 4:59
• Interesting for a novice like me to see how much changing variables can affect computation. May 30 '20 at 5:00
• @WillG You're welcome. I often solve such problems as this by inserting division points like Sqrt[5/x]`. One thing to keep in mind is that for (nonoscillatory) integrals over an infinite domain, Mathematica will choose a change of variables that maps the infinite domain to finite one, but it does not examine the integrand to figure out what a good one might be. May 30 '20 at 15:49