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

Question

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.9530743389378723807348355`16.*^6}

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

{1.45922, 5.95307736376147057651280734657987836949437727`17.*^6}

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

{2.85708, 5.9530773815454329104019275283442228729775598`18.*^6}

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

{6.62376, 5.95307738195314821215051085537129726640733153`19.*^6}

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

{16.2174, 5.95328982429028003608422815664470280663696859`20.*^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.95328982429028003608422815664470280665446132`30.*^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.

Answer 1

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.