2
$\begingroup$

I have five conditions con1 to con5 over the two-dimensional domain $t,u\in[o,2\pi]$. I want to compute this area $(con1 \cap con2 )\cup(con3 \cap con4 \cap con5)$ numerically. I try the two codes given below both with NIntegrate and Boole with the answers which are slightly different.

My question is:

Which of these results is more reliable? How can I modify the codes to get the most accurate result possible?

NIntegrate[ Boole[( con1 && con2 ) || (con3 && con4 && con5)], {t, 0, 2 Pi}, {u, 0, 2 Pi}, MaxRecursion -> 30, PrecisionGoal -> 30, Method -> "MonteCarlo"]
(*32.4243*)


NIntegrate[Boole[(con1 && con2) || (con3 && con4 && con5)], {t, 0, 2 Pi}, {u, 0, 2 Pi}]
(*34.3397*)

con1 = -4 (30 + 27 Cos[2 t] + 4 Cos[t - 5 u] + 3 Cos[2 t - 4 u] - 
      4 Cos[t] (-9 Cos[u] + 9 Cos[t] Cos[2 u] + 14 Cos[3 u]) - 
      18 Sin[2 t] Sin[2 u] + 4 Sin[t] (3 Sin[u] + 4 Sin[3 u])) <= 
   Cos[2 t - 8 u] - 6 Cos[2 t - 6 u] + 
    32 (-3 Cos[t] + Cos[t - 2 u]) Cos[u] + 
    24 (5 - 3 Cos[2 u]) Cos[2 u]^2 + 
    Cos[2 t] (27 - 66 Cos[2 u] + 92 Cos[4 u]) + 
    14 Sin[2 t] (3 Sin[2 u] - 5 Sin[4 u]) <= 
   4 (2 Abs[
        9 Cos[2 t] + 2 Cos[t - 5 u] + Cos[2 t - 4 u] - 
         2 (-7 + 4 Cos[t - 3 u] + 3 Cos[2 t - 2 u] - 4 Cos[t - u] + 
            Cos[2 u] + 9 Cos[t + 3 u])] + 
      Abs[6 + 9 Cos[2 t] + Cos[2 t - 4 u] - 4 Cos[t - 3 u] - 
        6 Cos[2 t - 2 u] + 20 Cos[t] Cos[u] - 10 Cos[2 u] - 
        4 Sin[t] Sin[u]] + 8 Cos[u]^2)    ;
con2 = 6 + 9 Cos[2 t] + Cos[t - 5 u] + Cos[2 t - 4 u] + 
    14 Cos[t] Cos[u] + 2 Sin[t] Sin[u] > 
   6 Cos[t - 3 u] + 6 Cos[2 t - 2 u] + 10 Cos[2 u] + 9 Cos[t + 3 u] ;
con3 = -32 Cos[u]^2 + (
    4 (9 Cos[2 t] + 2 Cos[t - 5 u] + Cos[2 t - 4 u] - 
       2 (-7 + 4 Cos[t - 3 u] + 3 Cos[2 t - 2 u] - 4 Cos[t - u] + 
          Cos[2 u] + 9 Cos[t + 3 u]))^2)/(
    9 Cos[2 t] + Cos[2 t - 4 u] - 
     2 (1 + 2 Cos[t - 3 u] + 3 Cos[2 t - 2 u] - 10 Cos[t] Cos[u] + 
        9 Cos[2 u] + 2 Sin[t] Sin[u])) <= 
   Cos[2 t - 8 u] - 6 Cos[2 t - 6 u] + 
    32 (-3 Cos[t] + Cos[t - 2 u]) Cos[u] + 
    24 (5 - 3 Cos[2 u]) Cos[2 u]^2 + 
    Cos[2 t] (27 - 66 Cos[2 u] + 92 Cos[4 u]) + 
    14 Sin[2 t] (3 Sin[2 u] - 5 Sin[4 u]) <= 
   4 (2 Abs[
        9 Cos[2 t] + 2 Cos[t - 5 u] + Cos[2 t - 4 u] - 
         2 (-7 + 4 Cos[t - 3 u] + 3 Cos[2 t - 2 u] - 4 Cos[t - u] + 
            Cos[2 u] + 9 Cos[t + 3 u])] + 
      Abs[6 + 9 Cos[2 t] + Cos[2 t - 4 u] - 4 Cos[t - 3 u] - 
        6 Cos[2 t - 2 u] + 20 Cos[t] Cos[u] - 10 Cos[2 u] - 
        4 Sin[t] Sin[u]] + 8 Cos[u]^2)    ;
con4 = (-3 Cos[t] + Cos[t - 2 u] - 2 Cos[u])^2 (9 Cos[2 t] + 
      Cos[2 t - 4 u] - 
      2 (1 + 2 Cos[t - 3 u] + 3 Cos[2 t - 2 u] - 10 Cos[t] Cos[u] + 
         9 Cos[2 u] + 2 Sin[t] Sin[u])) < 0   ;
con5 = Cos[
     u] (Sin[1/2 (t - 7 u)] + Sin[1/2 (3 t - 5 u)] - 
      12 Cos[(t + u)/2]^2 Sin[1/2 (t - 3 u)]) (Sin[1/2 (t - 3 u)] - 
      3 Sin[(t + u)/2])^2 (Sin[1/2 (t - 5 u)] + 
      3 Sin[1/2 (t + 3 u)]) < 0;
Improve this question
$\endgroup$
0

3 Answers 3

Reset to default
1
$\begingroup$

The Monte Carlo method

While the Monte Carlo method is not terribly accurate, it is fairly easy to understand. I don't know why the OP did not include the very useful information in the NIntegrate::maxp warning about the error estimate (the partitioning is suggested by the domain defined by the integrand; see below):

NIntegrate[
 Boole[(con1 && con2) || (con3 && con4 && con5)],
 {t, 0, 2 Pi}, {u, 0, 2 Pi},
 PrecisionGoal -> 30, 
 Method -> {"MonteCarlo",
   "MaxPoints" -> 50000,       (* default *)
   "Partitioning" -> {8, 8}}]  (* seems helpful *)

NIntegrate::maxp: The integral failed to converge after 50100 integrand evaluations. NIntegrate obtained 32.539178786473315` and 0.02103113210212079` for the integral and error estimates.

(*  32.5392  *)

So we know how accurate NIntegrate calculates the integral to be. If you believe it, it answers the OP's question, because the other integral is much farther away from the integral 32.539 than 0.021. We will discuss the OP's result of 32.4243 and why the other integral is less reliable below. First let's investigate the reliability of the Monte Carlo result. One nice(?) thing about Method -> "MonteCarlo" is that it does not BlockRandom[] the random number generator, so successive calls behave as a random variable. (You can control the reproducibility of the method through the "RandomSeed" suboption.) On a quick integral like this one, we can investigate the variability of the method.

(* Takes about 6.2 sec for me *)
Through[{Mean, StandardDeviation, MeanCI, MinMax, # &}@
   Table[Quiet@
     NIntegrate[Boole[(con1 && con2) || (con3 && con4 && con5)],
      {t, 0, 2 Pi}, {u, 0, 2 Pi},
      PrecisionGoal -> 30,
      Method -> {"MonteCarlo",
        "MaxPoints" -> 50000, "Partitioning" -> {8, 8}}],
    {100}]] // Labeled[
   Histogram[Last@#, 9],
   Grid@Transpose@{{"Mean", "S.D.", "C.I.", "Range"}, Most@#},
   Right] &

About the OP's 32.4243: Without the "Partitioning", the range in one trial of 100 runs was {32.3329, 32.6765}, which contains the OP's result, and the error estimate on a single run was about 0.067 corresponding with the wider range. The confidence interval was {32.4917, 32.5185}, which does not contain the OP's result.

Another consistency test is to raise the number of sample points. We know that the error in Monte Carlo should be proportional to 1/Sqrt[n] for n sample points.

NIntegrate[Boole[(con1 && con2) || (con3 && con4 && con5)],
 {t, 0, 2 Pi}, {u, 0, 2 Pi},
 PrecisionGoal -> 30,
 Method -> {"MonteCarlo", "MaxPoints" -> 1000000}]

NIntegrate::maxp: The integral failed to converge after 1000100 integrand evaluations. NIntegrate obtained 32.51101572370637` and 0.015049739619971662` for the integral and error estimates.

(*  32.511  *)

If we compare the error estimate of 0.67 for n = 50000 with the estimate 0.15 for n = 1000000, we calculate

0.067 Sqrt[50000/1000000]
(*  0.01498... *)

We don't expect strict equality because of the random sampling, but they're very close. The consistency supports the validity of the error estimates.

The (un)reliability of the Multidimensional Rule here

The theory supporting the default integration methods of NIntegrate[], such as the Multidimensional Rule used in the second integral in the OP, assume the integrand is analytic over the domain or at least within the interior of it. When there are singularities, a singularity handler may be invoked, but only when the singularity can be identified. If a singularity is not handled, one may expect the initial estimates for the integral and error to be inaccurate. With enough recursion in the global adaptive algorithm, the singularities might be confined to a small enough region that the error over that region contributes negligibly to the overall error. This is easiest when the singularities line up with coordinate boundaries.

In the OP's integral we have a bad singularity (jump discontinuity) along a complicated boundary. I don't think NIntegrate succeeds in partitioning the domain along this boundary.

Here is the region whose area the integral calculates:

(* Takes about 3 sec for me *)
bmesh = BoundaryDiscretizeRegion[
 ImplicitRegion[(con1 && con2) || (con3 && con4 && con5),
 {{t, 0, 2 Pi}, {u, 0, 2 Pi}}],
 MaxCellMeasure -> "Area" -> 0.001];

Show[bmesh, 
 GridLines -> {Subdivide[0, 2 Pi, 8], Subdivide[0, 2 Pi, 8]}, 
 Method -> {"GridLinesInFront" -> True}, Frame -> True]

Why "Partitioning" -> {8, 8} helps: We can see why "Partitioning" -> {8, 8} helps in the figure above. The squares at are all white or all blue have zero error, estimated and actual.

Region-based integration

By the way, we get another estimate of the integral from bmesh above:

Area@bmesh

(*  32.5073  *)

If we look at the region-based integration in NIntegrate[], there is no recursive refinement, based on error estimation of the integral. I don't believe there is a numerical error estimate at all, because there is only one integration rule applied and nothing to compare it to (some day, this assumption will be wrong, maybe today, but I think not yet). Therefore no error messages.

We can see below in these codes based on the other answers that PrecisionGoal is used in DiscretizeRegion. I do not know what error measure this controls. It is not related to NIntegrate and it does not seem to be related to the area of the region. Note also that the meshing routines work at machine precision. So while WorkingPrecision does not have its usual effect, it does change PrecisionGoal in NIntegration, which NIntegrate passes along to DiscretizeRegion.

(* These take 5-10 sec each *)
ir = ImplicitRegion[{(con1 && con2) || (con3 && con4 && con5), 
    0 <= u <= 2 Pi, 0 <= t <= 2 Pi}, {u, t}];
dr1 = Trace[
  NIntegrate[
   Boole[(con1 && con2) || (con3 && con4 && con5)],
   {u, t} \[Element] ir, WorkingPrecision -> 16],
  dr_DiscretizeRegion :> Return[Hold[dr], Trace],
  TraceInternal -> True]
(*
Hold[DiscretizeRegion[ImplicitRegion[<<1>>,{u,t}],
 PrecisionGoal -> 7., AccuracyGoal -> \[Infinity]]]
*)

dr2 = Trace[
  NIntegrate[
   Boole[(con1 && con2) || (con3 && con4 && con5)], {u, t} \[Element] 
    ir, WorkingPrecision -> 20],
  dr_DiscretizeRegion :> Return[Hold[dr], Trace],
  TraceInternal -> True];
(*
Hold[DiscretizeRegion[ImplicitRegion[<<1>>,{u,t}],
 PrecisionGoal -> 11., AccuracyGoal -> \[Infinity]]]
*)

dr3 = Trace[
   NIntegrate[1,
    {t, u} \[Element] 
     ImplicitRegion[(con1 && con2) || (con3 && con4 && con5),
      {{t, 0, 2 Pi}, {u, 0, 2 Pi}}]],
   dr_DiscretizeRegion :> Return[Hold[dr], Trace],
   TraceInternal -> True];
(*
Hold[DiscretizeRegion[ImplicitRegion[<<1>>,{u,t}],
 PrecisionGoal -> 7., AccuracyGoal -> \[Infinity]]]
*)

There is a limit to PrecisionGoal in DiscretizeRegion[]. If you set it too high (e.g. 14), you get this message:

DiscretizeRegion::drtol: Tolerance requested by the AccuracyGoal and PrecisionGoal options is too small to be achieved. Increasing to absolute tolerance 1.324082296537885`*^-7

One can also get this message in the following:

NIntegrate[1, {u, t} \[Element] ir, WorkingPrecision -> 32]

Note that PrecisionGoal -> 13 results in a computation that ran for several minutes before I killed it. You can reproduce the effect in NIntegrate[] as follows:

NIntegrate[1, {u, t} \[Element] ir, WorkingPrecision -> 13 + 9]

Update: Local adaptive integration

I thought about it earlier, but spent too much time on the above. But something about the region suggested it would be good to try the "LocalAdaptive" method:

PrintTemporary@Dynamic@Clock@Infinity;
AbsoluteTime[];
NIntegrate[
 Boole[(con1 && con2) || (con3 && con4 && con5)],
 {t, 0, 2 Pi}, {u, 0, 2 Pi},
 MaxRecursion -> 20,
 PrecisionGoal -> 6,
 Method -> "LocalAdaptive"]
AbsoluteTime[] - %%

(* PrecisionGoal -> integral (* time (sec) *)
   6 -> 32.50761472315272`   (*    5+s *)
   7 -> 32.50798110021814`   (*    6s  *)
   8 -> 32.508212967539805`  (*    9s  *)
   9 -> 32.5082805861034`    (*   16s  *)
  10 -> 32.50829781640405`   (*   38s  *)
  11 -> 32.50830725623989`   (*   90s  *)
  12 -> 32.508309526765856`  (*  339s  *)
  13 -> 32.50830998844479`   (* 1040s  *)
*)

We can see the rather slow convergence does not quite match the PrecisionGoal, but the convergence is steady.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you very much for the comprehensive answer. Sorry that I did not include the error since I am a beginner in Mathematica and did not think that error is that important. For the same reason (being a beginner in Mathematica ), now, with your answer, I cannot recognize which answer and which method is more accurate. Did you mean that the two already given answers are not accurate? $\endgroup$
    – Martha97
    Commented Jun 8, 2023 at 21:05
  • 1
    $\begingroup$ @Martha97 It's no problem about the error message. I was trying to emphasize that some of the messages have rather useful information. Some don't. People get used to them just meaning you did something wrong and ignore the details of the message. -- Sorry if the answer was confusing. "The other integral is much farther away from the integral 32.539 than 0.021" means the second integral (which internally uses the Multidimensional Rule) seems to be farther from the true value than the first integral (which uses Monte Carlo).... $\endgroup$
    – Michael E2
    Commented Jun 8, 2023 at 21:12
  • 1
    $\begingroup$ The value you got for Monte Carlo (32.4243) seems to be farther than normal from the average result of Monte Carlo which probably between 32.50 and 32.52. $\endgroup$
    – Michael E2
    Commented Jun 8, 2023 at 21:14
  • 1
    $\begingroup$ @Martha97 dr3 in the last section of my answer on Region-based integration shows the method used by cvgmt's solution. My Area@bmesh is an improvement of that approach. As I said, "I don't believe there is a numerical error estimate at all." That means Mma does not test how accurate the integral is and never gives an error about its accuracy, even if it is way off. OTOH, "MonteCarlo" tests the answer's accuracy, and issues an error when the accuracy does not meet PrecisionGoal. Further the error tells you that the answer is (probably) accurate to almost 3 digits. $\endgroup$
    – Michael E2
    Commented Jun 9, 2023 at 17:45
  • 1
    $\begingroup$ ¶ Summary: Region method -- you don't know how accurate it is, may be best, may be worst. Monte Carlo method -- you know how accurate it is, not great accuracy but pretty good and reliable. In the case at hand, the region methods are probably worse than Monte Carlo. (By "probably," I mean, yes, I could be wrong, but the evidence would lead me to pick Monte Carlo for my answer. In the time it took to write this, I ran MC for 10^6 points 100 times and got a 95% confidence interval for the integral to be between {32.5092, 32.5127}. $\endgroup$
    – Michael E2
    Commented Jun 9, 2023 at 17:45
2
$\begingroup$

To check if some result is accurate enough, one method is to increase "WorkingPrecision" and check if the result changes. The default "WorkingPrecision" is "MachinPrecisio"n, but without error control. We may therefore do the calculation with "WorkingPrecision" set to 16 and 20 and check if the results differ.

Further we max use "ImplicitRegion" to define the region:

ir = ImplicitRegion[{(con1 && con2) || (con3 && con4 && con5), 
    0 <= u <= 2 Pi, 0 <= t <= 2 Pi}, {u, t}];

With this:

NIntegrate[
 Boole[(con1 && con2) || (con3 && con4 && con5)], {u, t} \[Element] 
  ir, WorkingPrecision -> 16]

32.4927

NIntegrate[
 Boole[(con1 && con2) || (con3 && con4 && con5)], {u, t} \[Element] 
  ir, WorkingPrecision -> 20]

32.4927

As an increase of "WorkingPrecision" from 16 to 20 does not change the first 6 digits of the result, we may assume that these 6 digits are accurate.

Addendum

The integrand of "Boole[..]" is actually not necessary because Boole[..] evaluates to 1 in all the integration region . With an integrand of 1 we get:

NIntegrate[1, {u, t} \[Element] ir, WorkingPrecision -> 16]

32.4929

NIntegrate[1, {u, t} \[Element] ir, WorkingPrecision -> 20]

32.4929
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ It seems weird that the answers are in machine precision, when normally they should be precisions 16, 24 respectively. I think it's using mesh discretization for ir, which is available only in machine precision. NIntegrate really ought to warn about that. $\endgroup$
    – Michael E2
    Commented Jun 7, 2023 at 20:22
1
$\begingroup$

The two approachs without warning message.

NIntegrate[1, {t, u} ∈ 
  ImplicitRegion[(con1 && con2) || (con3 && con4 && con5), {{t, 0, 
     2 Pi}, {u, 0, 2 Pi}}]]

32.493

ImplicitRegion[(con1 && con2) || (con3 && con4 && con5), {{t, 0, 
     2 Pi}, {u, 0, 2 Pi}}] // DiscretizeRegion // Area

32.493

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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