# Why is regioncentroid not giving correct results?

I am trying to calculate and plot the x-axis of a volume centroid for the region between zx-plane and a function squared $$u(t,x,z)$$ from time $$t=15$$ to $$25$$. The function is a numerical solution to a particular wave equation, codes of which I will attach at the end.

The code I used to obtain the position of this volume centroid is as follows:

centroid[t_] :=
RegionPlot3D[
0 <= z <= 75 && 0 <= y <= uwave[t, x, z]^2, {z, 0, 75}, {y,
0,.2}, {x, -1, 1}] // DiscretizeGraphics // RegionCentroid;


Take one particular time-instance, for example $$t=18$$. The plot of our function squared in 3D looks like this:

and its code is

Plot3D[uwave[18, x, z]^2, {x, -1, 1}, {z, 0, 75}, PlotPoints -> 200,
PlotRange -> All, AxesLabel -> {"x", "z", "u"}]


Now clearly we would expect the centroid to have a x-axis slightly before $$x=20$$ as can be seen from the plot given. However, when I calculate it using the centroid function defined above, I find that

First[centroid[18]]


gives a value of 37.4976 which is clearly wrong. Where did I do wrong in the code?

For completeness, the numerical solution $$uwave(t,x,z)$$ is obtained by the following code which you may copy paste run, and get a result within 15 minutes.

Lx = 85;
Lz = 75;
Domain = Region[Rectangle[{-Lx, -Lz}, {Lx, Lz}]];

\[Rho]in = 5;
\[Rho]ex = 1;
w = 1;
p = 8;
\[Rho][x_] := (\[Rho]in - \[Rho]ex) (Sech[(x/w)^p])^2 + \[Rho]ex;
wx = 1;
wz = Sqrt[2];
u0saus[x_, z_] := x/wx*Exp[-(x/wx)^2 - (z/wz)^2];
VA[x_] := Sqrt[B0^2/\[Rho][x]]/Sqrt[B0^2/\[Rho][0]];
icX = u[0, x, z] == u0saus[x, z];
dicX = Derivative[1, 0, 0][u][0, x, z] == 0;
DX = DirichletCondition[u[t, x, z] == 0, True];
EqX = D[u[t, x, z], {t,
2}] - (VA[x])^2*(D[u[t, x, z], {x, 2}] +
D[u[t, x, z], {z, 2}]) == 0;
uwave = NDSolveValue[{EqX, icX, dicX, DX},
u, {t, 0, 25}, {x, z} \[Element] Domain, Method -> {
"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement",
"MeshOptions" -> {"MaxCellMeasure" -> 0.3},
"InterpolationOrder" -> {u -> 2}}}}]

• Please add definitions of B0 and u0saus (and any other missing values) to your code. Apr 28, 2022 at 19:13
• I have edited the original question to include a definition of u0saus. B0 is undefined because it will cancel out, so you may define whatever you like for B0. I have only included this thing B0 for physical completeness, a term which represents the background magnetic field. Apr 28, 2022 at 19:16
• You define u0sausX but use u0saus Apr 29, 2022 at 7:42
• Further, there is something wrong with your function "centriod". The centroid should be a 3 dim point, not a scalar. Apr 29, 2022 at 8:22
• Indeed, it is a point with 3 components as you said. I used first to extract its x-component Apr 29, 2022 at 8:26