2
$\begingroup$

Given a partial differential equation, we can NDSolve for its solution. In my case I got a function of three variables $u(t,x,z)$. I wish to study the case when we restrict $z=0$. Most importantly, I wish to animate how $u(t,x,0)$ changes with time $t$ when we plot the function against $x$, and at each instant of time $t$ add, along with the graph, a modified centroid point which is defined to be the usual centroid point of the function $(u(t,x,0))^2$. While I am able to plot the function in animation, I find difficulty in plotting the point. The code is as follows:

\[CapitalOmega] = Region[Rectangle[{-10, -10}, {10, 10}]]
IC11 = u[0, x, z] == x*Exp[-(x/s)^2 - (z/3)^2]
IC12 = Derivative[1, 0, 0][u][0, x, z] == 0
s = 1
\[Rho]0 = 3
\[Rho]max = 1
\[Rho][x_] := (\[Rho]0 - \[Rho]max) (Sech[x/s])^2 + \[Rho]max
CA[x_] := 1/(B0/Sqrt[\[Rho][0] \[Mu]0])*B0/Sqrt[\[Rho][x] \[Mu]0]
uwave1 = 
 NDSolveValue[{1/(CA[x])^2 D[u[t, x, z], {t, 2}] - 
     D[u[t, x, z], {x, 2}] - D[u[t, x, z], {z, 2}] == 0, 
   u[0, x, z] == x*Exp[-(x/s)^2 - (z/3)^2], 
   Derivative[1, 0, 0][u][0, x, z] == 0, 
   DirichletCondition[u[t, x, z] == 0, True]}, 
  u, {t, 0, 4 \[Pi]}, {x, z} \[Element] \[CapitalOmega], Method -> {
    "PDEDiscretization" -> {"MethodOfLines",
      "SpatialDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {"MaxCellMeasure" -> 0.2}, 
        "InterpolationOrder" -> {u -> 2}}}}];
Animate[Plot[(uwave1[t, x, 0])^2, {x, 0, 10}, 
  PlotRange -> {0, 0.2}], {t, 0, 4 \[Pi]}, AnimationRunning -> False]

where in the last part I plotted the animation. How should I proceed from here to add the centroid point into the animation? To clarify, the region is the region under the curve $u[t,x,0]^2$ at any specific time $t$ when plotted against x-axis.

Improve this question
$\endgroup$
2
  • $\begingroup$ centroid point only define on a regions instead of functions. $\endgroup$
    – cvgmt
    Commented Nov 26, 2021 at 10:56
  • $\begingroup$ The function seems to be symmetric uwave1[t,x,0]^2==uwave1[t,-x,0]^2, that's why centroid[t]==0 $\endgroup$
    – Ulrich Neumann
    Commented Nov 26, 2021 at 12:56

2 Answers 2

Reset to default
2
$\begingroup$

Edit

Using Polygon is faster.

center[t_] := Module[{pts},
   pts = Cases[Plot[uwave1[t, x, 0]^2, {x, 0, 10}, PlotRange -> All], 
     Line[a_] :> a, Infinity];
   RegionCentroid@
    Polygon@Join[{{pts[[1, 1, 1]], 0}}, 
      pts[[1]], {{pts[[1, -1, 1]], 0}}]];
ani = Animate[
  Plot[(uwave1[t, x, 0])^2, {x, 0, 10}, PlotRange -> {0, 0.2}, 
   Epilog -> {PointSize[Large], Red, Point[center[t]]}, 
   PerformanceGoal -> "Speed"], {t, 0, 4 π}, 
  AnimationRunning -> False, AnimationRate -> .05]

enter image description here

Original

To define the centroid of a plane region $\mathcal{R}$ we need to use double integral.

Assume that $\chi(\mathcal{R})$ is the characteristic function of the region $\mathcal{R}$,then the centroid $(\bar{x},\bar{y})$ is $$\bar{x}=\frac{\iint x\cdot\chi(\mathcal{R})}{\iint \chi(\mathcal{R})},\bar{y}=\frac{\iint y\cdot\chi(\mathcal{R})}{\iint \chi(\mathcal{R})}$$

For a positive function $f(x)$,it can be reduce to one-dimensional integral by

$$\bar{x}=\frac{1}{A}\int_{a}^{b}x f(x)\,\mathrm{d}x,\bar{y}=\frac{1}{A}\int_{a}^{b}\frac{1}{2}[f[x)]^2\,\mathrm{d}x$$ $$A=\int_a^b f(x)\,\mathrm{d}x$$

So we can define

χ[t_?NumericQ, x_?NumericQ, y_?NumericQ] := 
  Boole[0 <= x <= 10 && 0 <= y <= uwave1[t, x, 0]^2];
centroid[t_?NumericQ] := 
  NIntegrate[{x, y} χ[t, x, y], {x, 0, 10}, {y, 0, .2}]/
  NIntegrate[χ[t, x, y], {x, 0, 10}, {y, 0, .2}];

centroid2[
   t_?NumericQ] := {NIntegrate[x*uwave1[t, x, 0]^2, {x, 0, 10}]/
   NIntegrate[uwave1[t, x, 0]^2, {x, 0, 10}], 
   NIntegrate[1/2 uwave1[t, x, 0]^4, {x, 0, 10}]/
   NIntegrate[uwave1[t, x, 0]^2, {x, 0, 10}]};

but both of them are too slower to play the animation.

Here we discrete the region under the graph and use RegionCentroid to defined the cnetroid1[t].

centroid1[t_] := 
  RegionPlot[
     0 <= x <= 10 && 0 <= y <= uwave1[t, x, 0]^2, {x, 0, 10}, {y, 
      0, .2}] // DiscretizeGraphics // RegionCentroid;
Animate[Plot[(uwave1[t, x, 0])^2, {x, 0, 10}, PlotRange -> {0, 0.2}, 
  Epilog -> {PointSize[Large], Red, Point[centroid1[t]]}, 
  PerformanceGoal -> "Speed"], {t, 0, 4 π}, 
 AnimationRunning -> False, AnimationRate -> .05]
Improve this answer
$\endgroup$
5
  • $\begingroup$ The second part in your definiton centroid2[] , what kind of centroid does it describe? $\endgroup$
    – Ulrich Neumann
    Commented Nov 26, 2021 at 16:11
  • $\begingroup$ @UlrichNeumann The definaion of the centroid of a region is a double integral,it can be reduce to a one dimensional integral say centroid2[] $\endgroup$
    – cvgmt
    Commented Nov 26, 2021 at 23:59
  • $\begingroup$ Thank you, all agreed , I got it! $\endgroup$
    – Ulrich Neumann
    Commented Nov 27, 2021 at 13:08
  • $\begingroup$ Clever idea to use Polygon $\endgroup$
    – Ulrich Neumann
    Commented Nov 27, 2021 at 13:12
  • $\begingroup$ Could you please explain or refer me to an explanation as to what pts[1,1,1] and similar expressions mean inPolygon@Join[{{pts[[1, 1, 1]], 0}}, pts[[1]], {{pts[[1, -1, 1]], 0}}]]? $\endgroup$
    – Rescy_
    Commented Dec 5, 2021 at 6:42
2
$\begingroup$

If you are looking for the centroid in the range 0<x<10 proceed as follows.

The centroid xs of a function f[x] is definded xs Integrate[f[x],x]==Integrate[x f[x],x]

xs[t_?NumericQ]:= cp[t_?NumericQ] := NIntegrate[x uwave1[t, x, 0]^2 , {x, 0, 10} ]/ NIntegrate[  uwave1[t, x, 0]^2 , {x, 0, 10}]

Evaluation takes some time

Manipulate[
Plot[(uwave1[t, x, 0])^2, {x, 0, 10}, PlotRange -> {0, .5},GridLines -> {{xs[t]}, None}], {{t, 0}, 0, 4 Pi, Appearance -> "Labeled"}]

enter image description here

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.