Center of Mass

Code

centerOfMass:=(
mass = NIntegrate[f[x,y]Boole[regionFunction[x,y]],{x,xLeft,xRight},{y,yLeft,yRight}];
xMoment = NIntegrate[y f[x,y]Boole[regionFunction[x,y]],{x,xLeft,xRight},{y,yLeft,yRight}];
yMoment = NIntegrate[x f[x,y]Boole[regionFunction[x,y]],{x,xLeft,xRight},{y,yLeft,yRight}];
{yMoment/mass,xMoment/mass}
)

(*Given a point in the xy-plane, this function "tilts" it an angle theta relative to the support3D vector computed in initialComputations.*)
tilt[theta_,point_]:=(
relPoint = point-support3D;
N[point+(Cos[theta]-1)*(relPoint.unit)*unit+Sin[theta]*(relPoint.unit)*zUnit]
)

(*This function tilts all points in a collection of graphics data.*)
tiltGraphics[theta_,g_]:=(
g/.{x_Real,y_Real,0|0.}->tilt[theta,{x,y,0}]
)

graphics[theta_]:=(
Show[Graphics3D[{Opacity[0.75],tiltGraphics[theta,graphicsData]}],Graphics3D[supportGraphic],ViewPoint->{8,4,4},PlotRange->plotRange,Boxed->False]
)

ySupportPlane:=(
xMin=NMinimize[{x,regionFunction[x,support[[2]]],xLeft≤x≤xRight},x][[1]];
xMax=NMaximize[{x,regionFunction[x,support[[2]]],xLeft≤x≤xRight},x][[1]];
{Red,Polygon[{{xMin,support[[2]],0},{xMin,support[[2]],zMin},{xMax,support[[2]],zMin},{xMax,support[[2]],0}}]}
)
xSupportPlane:=(
yMin=NMinimize[{y,regionFunction[support[[1]],y],yLeft≤y≤yRight},y][[1]];
yMax=NMaximize[{y,regionFunction[support[[1]],y],yLeft≤y≤yRight},y][[1]];
{Red,Polygon[{{support[[1]],yMin,0},{support[[1]],yMin,zMin},{support[[1]],yMax,zMin},{support[[1]],yMax,0}}]}
)

initialComputations:=(
zMin=plotRange[[3]][[1]];

laminaPlot=DensityPlot[f[x,y],{x,xLeft,xRight},{y,yLeft,yRight},RegionFunction->(regionFunction[#1,#2]&),ColorFunction->colorFunction];

(*Embed the plot in the xy-plane of a 3D coordinate system.*)
laminaData = laminaPlot[[1]]/.{x_Real,y_Real}->{x,y,0};

cm = Append[centerOfMass,0];

support3D = Switch[balanceOn,"x",{support[[1]],cm[[2]],0},"y",{cm[[1]],support[[2]],0},_,Append[support,0]];

supportGraphic=Switch[balanceOn,"x",xSupportPlane,"y",ySupportPlane,_,{Red,Arrow[{Append[support,plotRange[[3]][[1]]],support3D}]}];

xLine = If[showXLine,ParametricPlot3D[{cm[[1]],t,0},{t,yLeft,yRight},PlotStyle->Green,RegionFunction->(regionFunction[cm[[1]],#2]&)][[1]],{}];
yLine = If[showYLine,ParametricPlot3D[{t,cm[[2]],0},{t,xLeft,xRight},PlotStyle->Green,RegionFunction->(regionFunction[#1,cm[[2]]]&)][[1]],{}];
cmPoint=If[showCenterOfMass,{Green,PointSize[0.02],Point[cm]},{}];

graphicsData = Graphics3D[{laminaData,cmPoint,xLine,yLine}][[1]];

unit = Normalize[cm-support3D];
zUnit = {0,0,1};

(*By considering torque, we find that the differential equation modeling this situation is
theta'' = C cos(theta)
where the constant C is the force acting on the lamina (assumed 1 in the calculation below) times the distance between the support point and the center of mass, and theta is the angle between the xy-plane and the vector between the support point and the center of mass.  We have Mathematica numerically compute solutions for equally spaced time values from 0 to tMax, stopping if our lamina ever hits vertical.*)
thetaSolution=NDSolve[{theta''[t]==-Norm[cm-support3D]*Cos[theta[t]],theta[0]==0,theta'[0]==0},theta[t],{t,0,tMax}];
lastTheta=If[stopAtVertical,-Pi/2,-Pi];
thetas=TakeWhile[Table[Evaluate[theta[t]/.thetaSolution][[1]],{t,0,tMax,.1}],#>lastTheta&];
)

singleFrame[theta_] := (
initialComputations;
Image[graphics@theta]
)

animation:=(
initialComputations;
Image[graphics@#]&/@thetas
)

Graphics

Change the right-hande side of f[x_,y_] to change the density function (constant density 1 by default).
Change the regionFunction to change the shape of the lamina.
xLeft, xRight, yLeft, and yRight control the bounds on the lamina.
plotRange determines how large the overall plotting box should be.

support is the point on which we’re balancing.
If balanceOn = “x” (with quotes), then we will balance on the plane x = support[[1]] (i.e., the first coordinate of support).
If balanceOn = “y” (still with quotes), then we will balance on the plane y = support[[2]] (i.e., the second coordinate of support).
For any other value of balanceOn, we will balance on the single point support.

showXLine and showYLine determine whether to draw the lines on the lamina corresponding to the components of the center of mass.
showCenterOfMass determines whether to draw the center of mass itself (which by definition is at the intersection of the two lines shown by showXLine and showYLine).

stopAtVertical determines whether to stop the animation when the lamina reaches vertical.  If False, then the lamina will appear to swing through to the other side of its support.

f[x_,y_]:=1;
regionFunction[x_,y_]:=(y^2+x^2≤1);
xLeft=-1;
xRight=1;
yLeft=0;
yRight=1;
plotRange={{-1.5,1.5},{-.5,1.5},{-1.5,1.5}};


support = {0,0.5};
balanceOn="y";

showXLine=False;
showYLine=True;
showCenterOfMass=False;

stopAtVertical=True;

(*Display parameters*)
colorFunction=Automatic; 
tMax = 10;

frames = animation;
ListAnimate[animation,12]

(*The following settings are used to generate the Key Concepts section of the webpage:
f[x_,y_]:=1;
regionFunction[x_,y_]:=(y^2≤x^3-2x+4);
xLeft=-4;
xRight=4;
yLeft=-4;
yRight=4;
plotRange={{-5,5},{-5,5},{-5,5}};

support = {2,-.5};

Additionally, the color function (RGBColor[1-#/2,1-#,1-#/2]&) was used to generate the animation in Visualizing the Example with non-constant density, so that denser parts of the lamina were darker.  By default (i.e., with a colorFunction of Automatic), Mathematica makes denser parts lighter instead.*)