# Double Numerical Integral - precision is not improving answer

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)?

• With respect to what is or is not oscillating, do you mean the integrand or the integral? (It cannot be the integral, I think, since it is a number, not a function.) My plots of the integrand look nothing like the plots above. Sep 8, 2019 at 14:08
• WorkingPrecision->15 is not going to improve your answers. You are probably confusing WorkingPrecision with PrecisionGoal. Voting to close. Sep 8, 2019 at 14:40
• RE Michael E2: when i talk about oscillations, i am of course talking about the integrand (how could an integral be oscillating?!), and these are plots of the integrand. how could you be plotting the integral for $b=0, a=4$? Sep 8, 2019 at 18:21
• [Site tip: Use @user to notify the users of your reply comment.] I can improve the integral that's posted, but, what I said before, "My plots of the integrand look nothing like the plots above." I don't want to write up an answer only to have you tell me that you made a mistake with the integrand. Could you post the codes you used for the plots, so I can compare them to the integral? Sep 8, 2019 at 18:52
• 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] Sep 8, 2019 at 20:57

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.

• Have to go to my real job now. Will respond to queries when I have the time. Sep 10, 2019 at 12:03
• Wow this is super useful! If I understand correctly, the main advantage of your method is that NIntegrate is sampling points ONLY within the region of the integral, whereas in my naive method i was telling Mathematica to integrate everywhere in r,s, but then afterward putting an If statement in the integrand so in practice only a small percentage of sampling points where actually being used? Sep 10, 2019 at 12:42
• @esches Yes, when sampling points cross a boundary, it adds a large amount to the error estimation. NIntegrate has to work much harder to refine the sampling around the boundaries to get the error within goals. Sep 10, 2019 at 15:22
• I think there might be a small mistake with your answer. When you do reduce to get 4 regions, the last two regions are overlapping, so your integral is double counting. You can fix this for the fourth region by adding the condition $r<1$. do you agree? Sep 12, 2019 at 11:00
• @cphys Thanks. You could also just edit such corrections. Oct 14, 2020 at 15:45