Surface integral with ImplicitRegion not working as expected?

Question

Bug introduced in 13.0 or earlier and persisting through 13.2.0 or later.

---

I'm trying to understand how ImplicitRegion works. Consider the following surface integral

$$\int_S z\,dS$$

Where $S$ is the lateral surface of the cylinder $y^2+z^2=4$ between two planes $x=y-3$ and $x=6-z$.

Using the parametrization $\mathbf{r}(x,\theta)=(x,2\cos(\theta), 2\sin(\theta))$ We have $$\int_S2\sin(\theta)|\frac{\partial \mathbf{r}}{\partial x}\times \frac{\partial \mathbf{r}}{\partial \theta}|\,dxd\theta=-8\pi$$ Notice the sign.

Now for my attempt in Mathematica

myReg2 = ImplicitRegion[y^2 + z^2 == 4, {{x, y - 3, 6 - z}, y, z}];
Integrate[z, {x, y, z} \[Element] myReg2]

> 8 \[Pi]

Even when I use my exact parametrization I get the wrong sign

myRegion = 
  ParametricRegion[{x, 2*Cos[\[Theta]], 
    2*Sin[\[Theta]]}, {{x, 2*Cos[\[Theta]] - 3, 
     6 - 2*Sin[\[Theta]]}, {\[Theta], 0, 2 \[Pi]}}];
Integrate[z, {x, y, z} \[Element] myRegion]

>8 \[Pi]

I guess I somehow have the limits confused or something, but I cant figure it out and it drives me crazy. Any input is appreciated.

Answer 1

This is a bug, but just wanted to post an answer to share what I found. Maybe it'll help someone at Wolfram figure out what's up.

Here's a plot of the region with the $z = 0$ plane shown.

Weird behavior when considering z + d

First, here's evidence that something is definitely up. I evaluated

Table[{d, Integrate[z + d, {x, y, z} \[Element] myReg2]}, {d, -1, 1, 1/3}]

and got

{{-1, -20 \[Pi]}, {-(2/3), -8 \[Pi]}, {-(1/3), 4 \[Pi]},
 {0, 8 \[Pi]}, {1/3, 28 \[Pi]}, {2/3, 40 \[Pi]}, {1, 52 \[Pi]}}

Plotted, this is

This should be linear in d. Further, we find for the surface area

Integrate[1, {x, y, z} \[Element] myReg2]

(* Out: 36 \[Pi] *)

suggesting (by linearity) that Mathematica "usually" thinks that the integral of z over the surface is (52 - 36) \[Pi] == 16 \[Pi]...except when computing it by itself, when it thinks it's 8 \[Pi].

Partitioning into component surfaces

Across the $z = 0$ plane

I also was curious what would happen when I split it up into a positive-z and negative-z surface. I got the following results:

myReg2pos = 
  ImplicitRegion[y^2 + z^2 == 4 && z > 0, {{x, y - 3, 6 - z}, y, z}];

myReg2neg = 
  ImplicitRegion[y^2 + z^2 == 4 && z < 0, {{x, y - 3, 6 - z}, y, z}]; 

posint = Integrate[z, {x, y, z} \[Element] myReg2pos]

(* 4 (18 + \[Pi]) *)

negint = Integrate[z, {x, y, z} \[Element] myReg2neg]

(* 4 (-18 + \[Pi]) *)

possurf = Integrate[1, {x, y, z} \[Element] myReg2pos]

(* 2 (-4 + 9 \[Pi]) *)

negsurf = Integrate[1, {x, y, z} \[Element] myReg2neg]

(* 8 + 18 \[Pi] *)

{possurf, negsurf} // N

(* {48.54866776, 64.54866776} *)

{posint, negint} // N

(* {84.56637061, -59.43362939} *)

Note that linearity works here, at least (whew!), and that Mathematica knows that the surface area of the lower half is higher than the surface area of the upper half. Yet, the integral of $z$ over the upper half is higher in absolute value than the integral of $z$ over the lower half. While possible in principle, given that we know how this region looks, we know something's wrong.

Along $x$

I also split the surface into three parts like this:

myReg2left = ImplicitRegion[y^2 + z^2 == 4, {{x, y - 3, -1}, y, z}];

myReg2mid = ImplicitRegion[y^2 + z^2 == 4, {{x, -1, 4}, y, z}];

myReg2right = ImplicitRegion[y^2 + z^2 == 4, {{x, 4, 6 - z}, y, z}];

{surfleft, surfmid, surfright} = 
 Integrate[1, {x, y, z} \[Element] #] & /@ {myReg2left, myReg2mid, 
   myReg2right}

(* {8 \[Pi], 20 \[Pi], 8 \[Pi]} *)

{intleft, intmid, intright} = 
 Integrate[z, {x, y, z} \[Element] #] & /@ {myReg2left, myReg2mid, 
   myReg2right}

(* {0, 0, 8 \[Pi]} *)

The problem therefore arises, perhaps unsurprisingly, when the $x \leq 6 - z$ constraint is encountered.

Answer 2

This is a bug.

Integrate[z, {x, y, z} ∈ 
  DiscretizeRegion[myReg2, MaxCellMeasure -> .01]]

-25.0531