Double Numerical Integral - precision is not improving answer

Question

Consider this two-dimensional integral

 NIntegrate[
 r^3 4^a/ ((1 + r^2 - 2 r Cos[s])^4)
      ChebyshevU[b, -((-1 + r)/Sqrt[1 + r^2 - 2 r Cos[s]])] ((
     1 + r - 2 Sqrt[r] Cos[s/2])/Sqrt[1 + r^2 - 2 r Cos[s]])^
    a If[1 + r^2 - 2 r Cos[s] < 1 && 
      1 - r Cos[s] < 4/
       5 && (r^2 > 1 || 1/2 > r Cos[s] || r Cos[s] > 4/5), 1, 0] Sin[
     s]^2  , {r, 0, ∞}, {s, 0, π}]

which I would like to integrate for various values of $a$ and $b$.

If I compute this integral with WorkingPrecision->15 for $b=0,a=4$, then I get 13.9927829 with error 0.0177517. Similarly, for $b=10,a=20$ I get -0.212035 with error 0.003404. In both cases, if I improve the WorkingPrecision, PrecisionGoal, and/or AccuracyGoal (which are the usual tricks I know that help reduce error) though, then the error does not decrease! The problem does not seem to be oscillations, as the same problem occurs for both $b=0$ and $b=10$, and I do not think the integrand oscillates when $b=0$ (even though it might when $b=10$ due do the Chebyshev polynomial).

How can I improve the integral such that I can get arbitrary accurate results for any reasonable (i.e. not huge) value of $a,b\geq0$ (even if I have to wait longer to get more accuracy, I don't mind that)?

Answer 1

Code, with updated Reduce[] and integration Method:

integrand = 4^a r^3/((1 + r^2 - 2 r Cos[s])^4) *
   ChebyshevU[b, -((-1 + r)/Sqrt[1 + r^2 - 2 r Cos[s]])] *
   ((1 + r - 2 Sqrt[r] Cos[s/2])/Sqrt[1 + r^2 - 2 r Cos[s]])^a *
   (*Boole[1+r^2-2 r Cos[s]<1 && 1-r Cos[s]<4/5 && (r^2>1 || 1/2>r Cos[s] || r Cos[s]>4/5)]*)
   Sin[s]^2;
wp = MachinePrecision; (* set working precision *)
wp = 24;               (* set working precision *)
pg = wp/2.9;           (* set precision goal *)
toNInt[_[s1_, ___, s2_] && _[r1_, ___, r2_]] := (* integrate over a region *)
  NIntegrate[integrand, {s, s1, s2}, {r, r1, r2},
   PrecisionGoal -> pg, WorkingPrecision -> wp,
   Method -> {"MultidimensionalRule", "Generators" -> 9}];
PrintTemporary@Dynamic@{foo = Clock[Infinity]};  (* monitors time *)
Block[{b = 10, a = 20, echo},
  echo = (Print[foo, ": integrating ", #]; #) &; (* prints time each integral starts *)
  subintegrals = Replace[
    Reduce[#, r] & /@ 
     BooleanConvert@
      Reduce[(1 + r^2 - 2 r Cos[s] < 1 && 
          1 - r Cos[s] <  4/5 &&
          (r^2 > 1 || 1/2 > r Cos[s] || r Cos[s] > 4/5)) && 
        0 < s < Pi && r > 0, {r}, Reals],
    HoldPattern@Or[i__] :> toNInt@*echo /@ {i}
    ];
  integral = SetPrecision[subintegrals, pg] // Total;
  integralprec = Precision[integral];
  {integral, integralprec}
  ] // AbsoluteTiming
(*
  0.335269: integrating 0<s<ArcCos[4/5]&&Sec[s]/5<r<Sec[s]/2

  0.335269: integrating 0<s<ArcCos[4/5]&&(4 Sec[s])/5<r<2 Cos[s]

  0.335269: integrating ArcCos[4/5]<=s<π/3&&1<r<2 Cos[s]

  0.335269: integrating π/3<=s<ArcCos[1/Sqrt[10]]&&Sec[s]/5<r<2 Cos[s]

  0.335269: integrating ArcCos[4/5]<=s<π/3&&Sec[s]/5<r<Sec[s]/2

Out[]=
  {0.427674, {-0.02578, 3.90475}}
*)

Here are the results for increasing working precision. One can infer the accuracy of a previous result by comparing how many digits agree with the subsequent results. MachinePrecision results are obtained quickly, but increasing WorkingPrecision rapidly slows down the computation.

{time,
 {integral,                      pg}}

{0.427674`,                      (* wp=MachinePrecision *)
 {-0.0257816408128350764,        3.9047482670549747`}}
           ^ 6 digits of precision > pg=3.9
{4.366194`,                      (* wp=24 *)
 {-0.0257816195648316161,        6.679027290327976`}}
              ^ 9 digits of precision > pg=6.7
{94.906609`,                     (* wp=36 *)
 {-0.025781619559796699077118,   10.816958324724922`}}
                    ^ 15 digits of precision > pg=10.8
{369.352305`,                    (* wp=42 *)
 {-0.025781619559796682025612,   12.885923841966303`}}

The domain of integration defined by If[cond, 1, 0] can be broken down by Reduce into component subregions, if we solve for r in terms of s (indicated by the order of the variables {s, r} in the second argument of Reduce[]). It splits the domain into five subregions:

regs = Reduce[(1 + r^2 - 2 r Cos[s] < 1 && 
    1 - r Cos[s] < 4/5 &&
    (r^2 > 1 || 1/2 > r Cos[s] || r Cos[s] > 4/5)) && 
  0 < s < Pi && r > 0, {r}]
(*
(0 < s < ArcCos[4/5] && Sec[s]/5 < r < Sec[s]/2) ||
 (0 < s < ArcCos[4/5] && (4 Sec[s])/5 < r < 2 Cos[s]) ||
 (ArcCos[4/5] <= s < π/3 && 1 < r < 2 Cos[s]) ||
 (π/3 <= s < ArcCos[1/Sqrt[10]] && Sec[s]/5 < r < 2 Cos[s]) ||
 (ArcCos[4/5] <= s < π/3 && Sec[s]/5 < r < Sec[s]/2)
*)

Show[
 RegionPlot[regs,
  {s, 0, Pi/2}, {r, 0, 2.5}],
 ContourPlot[Evaluate@DeleteDuplicates@Flatten@Cases[regs,
      HoldPattern@Inequality[x1_, _, x_, _, x2_] | 
        Less[x1_, x_, x2_] :> {x == x1, x == x2}, Infinity],
  {s, 0, Pi/2}, {r, 0, 2.5}]
 ]

The function toNInt converts a subregion into a call to NIntegrate. The code assumes the output of Reduce is in form shown above.