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.
WorkingPrecision->15
is not going to improve your answers. You are probably confusingWorkingPrecision
withPrecisionGoal
. Voting to close. $\endgroup$integrand = 4^a r^3/((1+r^2-2r Cos[s])^4)ChebyshevU[b, -((-1+r)/Sqrt[1+r^2-2r Cos[s]])] ((1+r-2Sqrt[r]Cos[s/2])/Sqrt[1+r^2-2rCos[s]])^a Sin[s]^2; wp = MachinePrecision; wp = 24; pg = wp/2.6; toNInt[_[s1_,___,s2_] && _[r1_,___,r2_]] := NIntegrate[integrand, {s,s1,s2}, {r,r1,r2}, PrecisionGoal -> pg, WorkingPrecision -> wp]; Block[{b = 10, a = 20, echo}, subintegrals = Replace[Reduce[(1+r^2-2r 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}], HoldPattern@Or[i__] :> toNInt /@ {i}]; integral = SetPrecision[subintegrals, pg] // Total]
$\endgroup$