5
$\begingroup$

I have this two-variable function $$f(x,y)= (8 \cos (x+y)+7)\cos \left(\frac{x}{2}\right)+\cos \frac{x-2 y}{2}+2 \cos \left(\frac{3 x}{2}\right) $$ where $0<x,y<\pi$. I want to calculate numerically the area for which the function is negative $f(x,y)<0$. I use this code

NIntegrate[
 Boole[2 Cos[(3 x)/2] + Cos[1/2 (x - 2 y)] + 
    Cos[x/2] (7 + 8 Cos[x + y]) < 0], {x, 0, Pi}, {y, 0, Pi}]

and it gives the answer $3.49458$, but Mathematica gives the following warnings. Are there any other ways to calculate this value that is more reliable and more accurate than this method?

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 3.494581434480605 and 0.002397336775896384 for the integral and error estimates.

Improve this question
$\endgroup$
6
  • $\begingroup$ Your f[x, y] is a constant rather than a function of x and y. $\endgroup$
    – Bob Hanlon
    Commented Jun 17, 2021 at 19:03
  • $\begingroup$ @BobHanlon Thanks, it was a typo mistake. $\endgroup$
    – user80186
    Commented Jun 17, 2021 at 23:21
  • $\begingroup$ There is an exact solution to your integral $-2 \text{Li}_2\left(-\frac{1}{2}\right)+\text{Li}_2\left(\frac{1}{4} \left(i \sqrt{3}+1\right)\right)+\text{Li}_2\left(\frac{1}{4} \left(1-i \sqrt{3}\right)\right)+\frac{2 \pi ^2}{9}$ $\endgroup$
    – yarchik
    Commented Jun 18, 2021 at 13:12
  • $\begingroup$ @yarchik I think you should add that and its derivation as an answer. $\endgroup$
    – Kiro
    Commented Jun 18, 2021 at 14:03
  • 4
    $\begingroup$ @SaraChem @Kiro I was simply a bit more patient than @eyorble. I used u = Integrate[upperCurve, {x, 0, Pi}] and v = Integrate[lowerCurve, {x, 0, 2 Pi/3}], followed by u-v. The lowerCurve and the upperCurve are defined in eyorble's post. $\endgroup$
    – yarchik
    Commented Jun 18, 2021 at 19:46

4 Answers 4

Reset to default
4
$\begingroup$

Let us name the function of interest f[x,y]:

f[x_, y_] := (8 Cos[x + y] + 7) Cos[x/2] + Cos[(x - 2 y)/2] + 2 Cos[3 x/2]

Attempt to ContourPlot to find the zero lines:

ContourPlot[f[x, y] == 0, {x, 0, Pi}, {y, 0, Pi}]

Contour plot of zero lines for f[x,y]

Simple numerical evaluation determines that the negative region is the inside of these two curves. Notice that the ContourPlot is quite rapid and has very clean lines. Interesting coincidence. Perhaps there exists an analytical solution to these curve lines?

sol = Solve[{f[x, y] == 0}]

This returns a list of 4 possible curves, while also stating that some solutions may be missing. By manual inspection (such as by using Plot), we can find that the 2nd and 4th solutions are of interest to us, so we shall label them:

upperCurve = y /. sol[[2]];
lowerCurve = y /. sol[[4]];
Plot[{upperCurve, lowerCurve}, {x, 0, Pi}, PlotRange -> {0, Pi}]

Plots of the 2nd (blue) and 4th (yellow) curves

Checking the curves manually by plotting them against the original ContourPlot, we see that upperCurve matches the upper line for the whole domain, and that lowerCurve matches the lower line up until it reaches its minimum.

Find the minimum of the lowerCurve:

FindMinimum[{lowerCurve, 0 < x < Pi}, x, WorkingPrecision -> 25]

{1.872299341324760554288429*10^-8, {x -> 2.094395111754692173633430}}

The warning about a small imaginary part is of little concern here, but you can increase the WorkingPrecision and PrecisionGoal if you would like more digits.

Michael Seifert also pointed out that an exact form can be found for this solution by applying TrigFactor to f[x,y]:

TrigFactor[f[x,y]]

2 (Cos[x/2 - y/2] + 2 Cos[x/2 + y/2]) (2 Cos[x + y/2] + Cos[y/2]) == 0

The lower line happens to correspond to the second variable factor in this expression, and its minimum is found when y is set to 0 and solved.

Solve[{(2 Cos[x + y/2] + Cos[y/2]) == 0 /. y -> 0, 0 < x < Pi}, x]

{{ x -> 2 Pi/ 3 }}

Integrate the area below the 2nd curve over the whole domain minus the area under the 4th curve for 0 through 2.094...

NIntegrate[upperCurve, {x, 0, Pi}, WorkingPrecision -> 25] - 
 NIntegrate[lowerCurve, {x, 0, 2.09439511175469217363342977478168904781`25.},
   WorkingPrecision -> 25]

3.49805583366099845069196

Or with the exact form, we can see a slightly different answer:

NIntegrate[upperCurve, {x, 0, Pi}, WorkingPrecision -> 25] - 
 NIntegrate[lowerCurve, {x, 0, 2 Pi/3}, WorkingPrecision -> 25]

3.49805583366099836305434

While this method is not universally applicable, it does work for this function and is much faster than Area or direct application of NIntegrate and Boole for high precisions. As yarchik notes, you can swap NIntegrate for Integrate here to acquire an exact solution, though it takes a bit longer to evaluate.

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ Using TrigFactor on f[x,y] shows that the two curves are $\cos((x-y)/2) + 2 \cos((x+y)/2) = 0$ (upper) and $2 \cos (x+y/2)+ \cos(y/2) = 0$ (lower). I don't know if that really helps, though. $\endgroup$
    – Michael Seifert
    Commented Jun 17, 2021 at 20:10
  • 1
    $\begingroup$ If nothing else, it shows that the intersection of the lower curve with the $x$-axis is when $2 \cos x + 1 = 0$, or $x = 2 \pi/3$. $\endgroup$
    – Michael Seifert
    Commented Jun 17, 2021 at 20:12
  • $\begingroup$ There are a lot of inconsistencies in your answer. What is y that you are integrating? $\endgroup$
    – yarchik
    Commented Jun 18, 2021 at 11:51
  • $\begingroup$ @yarchik The y in the integration should be read as y /. sol[[2]] (which is the upper curve from the contour plot) or y /. sol[[4]] (which is the lower curve from the contour plot). These are just explicit forms (in x) of the curves where f[x,y]==0, selected to match the curves we are interested in. $\endgroup$
    – eyorble
    Commented Jun 18, 2021 at 12:02
  • $\begingroup$ @yarchik I have modified the answer to hopefully clarify what's being integrated here. $\endgroup$
    – eyorble
    Commented Jun 18, 2021 at 12:11
1
$\begingroup$ 
Clear["Global`*"]

RegionPlot[
 2 Cos[(3 x)/2] + Cos[1/2 (x - 2 y)] + Cos[x/2] (7 + 8 Cos[x + y]) < 0,
 {x, 0, Pi}, {y, 0, Pi},
 Frame -> True]

enter image description here

rgn = ImplicitRegion[{2 Cos[(3 x)/2] + Cos[1/2 (x - 2 y)] + 
       Cos[x/2] (7 + 8 Cos[x + y]) < 0 && 0 < x < Pi && 
     0 < y < Pi}, {x, y}];

Area[rgn, WorkingPrecision -> MachinePrecision]

(* 3.49805 *)

Area[rgn, WorkingPrecision -> 15]

(* 3.49805 *)
Improve this answer
$\endgroup$
1
$\begingroup$

A similar way is

reg=ImplicitRegion[
   2 Cos[(3 x)/2] + Cos[1/2 (x - 2 y)] + Cos[x/2] (7 + 8 Cos[x + y]) <
     0, {{x, 0, Pi}, {y, 0, Pi}}];
NIntegrate[1, {x, y} ∈ reg]

3.49485

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks. Does WorkingPrecission matter in this method? I see a small difference in the different answers here obtained by different methods others suggested; if I need the most accurate one, which one is more reliable? $\endgroup$
    – user80186
    Commented Jun 18, 2021 at 12:54
0
$\begingroup$

Method ->"LocalAdaptive" evaluates without errormessage

NIntegrate[Boole[2 Cos[(3 x)/2] + Cos[1/2 (x - 2 y)] +Cos[x/2] (7 + 8 Cos[x + y]) < 0], {x, 0, Pi}, {y,0, Pi}, Method -> "LocalAdaptive"]
(*3.49818*)
Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks. Does WorkingPrecission matter in this method? I see a small difference in the different answers here obtained by different methods others suggested; if I need the most accurate one, which one is more reliable? $\endgroup$
    – user80186
    Commented Jun 18, 2021 at 12:55
  • 1
    $\begingroup$ @SaraChem Sorry, I was offline some days. If you add the option , IntegrationMonitor :> ((errors = Through[#1@"Error"]) &) inside NIntegrate , it's possible to evaluate the integrationerror Total@errors afterwards. Might be helpful. $\endgroup$
    – Ulrich Neumann
    Commented Jun 20, 2021 at 14:27

Your Answer

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