# How to calculate the total area of $\frac{\cos (2 y) \cos (2 x-y)+\cos (3 y)}{\cos (2 x-y)+\cos (y)}$ enclosed in the range $[0,1]$?

I have this two-variable function $$f(x,y)=\frac{\cos (2 y) \cos (2 x-y)+\cos (3 y)}{\cos (2 x-y)+\cos (y)}$$ where $$0. I would like to calculate numerically the total area of the surface where $$0, in other words the area enclosed in the range $$[0,1]$$?

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

• The question is unclearly formulated. What area do you mean? May 27 at 8:06
• Indeed is not clear if you mean the area of the x,y domain where that relation holds, of the area of the surface defined by f May 28 at 0:05

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


we define a function that is 1 in the allowed region and 0 otherwise. Such function is calla a predicate:

f1[x_, y_] = If[0 <= f[x, y] <= 1, 1, 0];

NIntegrate[f1[x, y], {x, 0, 2 Pi}, {y, 0, 2 Pi},
Method -> "MonteCarlo"]

(* 12.8533 *)


You may ask why I used: Method -> "MonteCarlo"? Sharp edges pose a problem for ordinary numerical methods. Not so for MonteCarlo. The probability to get a 1 is equal to the ratio of the allowed area to the total area. Therefore the allowed area= TotalArea * Number of successful trials / total number of trials.

• this is not what he asked. This is the volume between the plane z=0 and the surface z=f(x,y) in the region where 0<f(x,y)<1 May 26 at 19:27
• @Giorgio Busoni Sorry, you are wrong. What you claim is f1[x_, y_] = If[0 <= f[x, y] <= 1, f[x,y], 0] not what I wrote. May 26 at 19:32
• yes sorry, still not what he asked for May 26 at 19:33
• r = ImplicitRegion[0 < (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + Cos[y]) && (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + Cos[y]) < 1 && x >= 0 && x <= 2*Pi && y >= 0 && y <= 2*Pi, {x, y}]; RegionMeasure[r] produces 12.7720179163377. May 27 at 4:15

Following Lukas formula for the area, you can use NIntegrate with HeavisideTheta functions to define the domain. You can add accuracy and/or precision options if required.

  NIntegrate[
Sqrt[1 + D[f[x, y], x]^2 + D[f[x, y], y]^2] HeavisideTheta[
f[x, y]] HeavisideTheta[1 - f[x, y]], {x, 0, 2 Pi}, {y, 0, 2 Pi}]
(*39.33*)

• You omitted an error message "NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. 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 39.33625175657934 and 0.5710104625871588 for the integral and error estimates." May 26 at 19:40
• As I wrote, you need to play with the integration options May 26 at 20:02

You calculate the length of a curve $$f(x)$$ between $$a$$ and $$b$$ by $$\int_a^b \sqrt{1+f'(x)^2}dx$$, see https://en.wikipedia.org/wiki/Arc_length.

The 2d generalisation would logically be $$\int_A \sqrt{1+(\partial_x f(x,y))^2+(\partial_y f(x,y))^2}dA$$. Implementing in mathematica would go like

f[x_,y_] := (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] + Cos[y])
R = ImplicitRegion[{0 < x < 2 Pi, 0 < y < 2 Pi, 0 < f[x, y] < 1}, {x, y}]
NIntegrate[Sqrt[1 + D[f[x, y], x]^2 + D[f[x, y], y]^2], {x, y} \[Element] R]


This gives an answer of about 9062.5, which seems to me like like too high of a number, so I might have made a mistake.

Edit: Second try. This answer seems more reasonable and agrees with the other solutions. Some error messages are still given about convergence about which I know very little.

integrand = Sqrt[1 + D[f[x, y], x]^2 + D[f[x, y], y]^2]
NIntegrate[If[0 < f[x, y] < 1, integrand, 0], {x, 0, 2 Pi}, {y, 0, 2 Pi}]


Edit: I thought you were interested in the area of the surface of your function. If you were interested in the are of the domain where $$0, I misunderstood, and you should follow Daniel's answer.

• Executing you code, I obtain NIntegrate[Sqrt[1 + (\!$$\*SubscriptBox[\(\[PartialD]$$, $$x$$]$$f[x, y]$$\))^2 + (\!$$\*SubscriptBox[\(\[PartialD]$$, $$y$$]$$f[x, y]$$\))^2], {x, y} \[Element] R]. May 26 at 19:48

Let us use a straightforward approach.

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


Its plot is built by

Plot3D[f,{x,0,2*Pi},{y,0,2*Pi}]


The result of

r = ImplicitRegion[0 < (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] +
Cos[y]) && (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] +
Cos[y]) < 1 && x >= 0 && x <= 2*Pi && y >= 0 && y <= 2*Pi, {x, y}]
Region[r]


shows the region of the integration. Next,

FunctionPeriod[f, x]


Pi

and

FunctionPeriod[f,y]


2*Pi

Hence, the required area approximately equlals

2*NIntegrate[Sqrt[1 + D[f, x]^2 + D[f, y]^2]*
Boole[0 < (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] +
Cos[y]) && (Cos[2 x - y] Cos[2 y] + Cos[3 y])/(Cos[2 x - y] +
Cos[y]) < 1], {x, 0, Pi}, {y, 0, 2*Pi}, PrecisionGoal -> 2,
AccuracyGoal -> 2, WorkingPrecision -> 15, Method -> "GlobalAdaptive", MaxRecursion -> 30]


39.6109795138184

and a warning "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."

Too long for a comment:

Notice: Hacky hack. I like to see what the function looks like, so...

Plot3D[ (Cos[2 y] Cos[2 x - y] + Cos[3 y])/( Cos[2 x - y] + Cos[y]), {x, 0, 2 π}, {y, 0, 2 π}, PlotRange -> {0, 1}, ClippingStyle -> None, PlotPoints -> {100, 100}, ColorFunction -> Function[{x, y, z}, Hue[z]]]

So, it's pretty clear the surface area is more than 6 (2 pi). Either use Lucas's formula, or define the region and use the surface area command built into Mathematica. For the hackiest hack, just double the surface area of the second and third objects from the top