2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Sep 8 at 14:08
  • $\begingroup$ WorkingPrecision->15 is not going to improve your answers. You are probably confusing WorkingPrecision with PrecisionGoal. Voting to close. $\endgroup$ – Anton Antonov Sep 8 at 14:40
  • $\begingroup$ 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$? $\endgroup$ – esches Sep 8 at 18:21
  • 1
    $\begingroup$ [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? $\endgroup$ – Michael E2 Sep 8 at 18:52
  • 1
    $\begingroup$ 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$ – Michael E2 Sep 8 at 20:57
4
$\begingroup$

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}]
 ]

enter image description here

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

$\endgroup$
  • $\begingroup$ Have to go to my real job now. Will respond to queries when I have the time. $\endgroup$ – Michael E2 Sep 10 at 12:03
  • $\begingroup$ 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? $\endgroup$ – esches Sep 10 at 12:42
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Sep 10 at 15:22
  • $\begingroup$ 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? $\endgroup$ – esches Sep 12 at 11:00
  • $\begingroup$ @esches Oops, you're right. I didn't realize Reduce would produce overlapping regions. It works if the domain Reals is specified. I'll update with the accurate results when it's done running. $\endgroup$ – Michael E2 Sep 12 at 12:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.