Moments of Inertia:  3D Integration

(copyright Kieran Mullen 1996)

You've seen in class that moments of inertia can be calculated by doing an integral over  the distribution of mass in an object.  The idea is pretty simple - you just sum up the  bits of mass that comprise the object, weighting them by the square of their distance from the axis of rotation.
      Start by evaluating the expression below to fire up the Mathematica kernel.

2+2

Cubes and Parallelpipeds:

First we consider the case of a cube.  We'd like to calculate its mass.  The simple way is to multiply its volume,  V, by its density, ρ.  Since the volume of a cube with sides of length L  is (L^3), we have that:

                    (Subscript[M, cube  ])=  ρ (L^3)

where we have treated the density ρ as a constant.  (We will assume throughout this exercise that the density is constant.)  However, we may consider the single cube as being made up of many little cubes:

We can view the big cube as a "sum" of little cubes.  To get the total mass, we have to sum over each individual mini-cube.  How do we do that?  

Well, we can see that each mini- cube has an x,  y and z coordinate.  Let's denote a given mini-cube by its coordinates,  {x,y,z}.  Each mini-cube has a mass, (Subscript[Δm, {x,y,z}]) which is given by the density times the volume of the minicube, so that 

              Δ(Subscript[m, {x,y,z}])=(Subscript[ρ   ΔV, {x,y,z}])
              
 In order to get the total mass we must sum up all the masses of the mini-cubes.  
 
 (Subscript[(∑_(x=(Subscript[x, i]))^x), f])   (Subscript[(∑_(y=(Subscript[y, i]))^y), f])   (Subscript[(∑_(z=(Subscript[z, i]))^z), f])    Δ(Subscript[m, {x,y,z}]) =
 
                                       (Subscript[(∑_(x=(Subscript[x, i]))^x), f])   (Subscript[(∑_(y=(Subscript[y, i]))^y), f])   (Subscript[(∑_(z=(Subscript[z, i]))^z), f])  ρ   (Subscript[ΔV, {x,y,z}])
                                       
 The volume of the mini-cube is given by the length of its sides,
 
                                                       (Subscript[ΔV, {x,y,z}]) =  Δx  Δy  Δz
 
 where we have assumed that they are all the same size.  If we take the limit that the sides of the mini-cubes become very small, we have to sum (integrate) over them in the x direction, the y direction, and the z direction, so that we have three integrals or a triple integral.
 
    M=(∫_(-(L/2))^(L/2))(∫_(-(L/2))^(L/2))(∫_(-(L/2))^(L/2))  (Subscript[ρ, 0])  dzdydx

We do this triple integral in Mathematica with a single "Integrate" command.   Evaluate this command below:

(∫_(-(L/2))^(L/2))(∫_(-(L/2))^(L/2))(∫_(-(L/2))^(L/2))  (Subscript[ρ, 0])  \[DifferentialD]z\[DifferentialD]y\[DifferentialD]x

We are integrating a constant, "ρ", over the volume of a cube.  How do we know that it is a cube?  The limits of the integrals specify the extent of the object.  For example, x has a minium possible value of -L/2, and a maximum value of L/2, and the same is true for y and x. Convince yourself that this will just cover the volume of the cube.   Note that we get the correct answer, which is reassuring!

We must also calculate the moment of inertia.  If we are rotating about the z-axis, then the distance of any point (x,y,z)  from the z-axis, is just x^2 + y^2.  (Make sure you understand this!).  We integrate over the same range as before:

(∫_(-(L/2))^(L/2))(∫_(-(L/2))^(L/2))(∫_(-(L/2))^(L/2))  (Subscript[ρ, 0]) ((x^2)+(y^2))\[DifferentialD]z\[DifferentialD]y\[DifferentialD]x

Ok.  We have the mass and the mment of inertia in terms of rho and L.  However, it is customary to write the moment of inertia in terms of the mass of the object.  We can do this by factoring out the expression for the mass from the expression for the moment of inertia.  This can be done easily in Mathematica, if you remember a simple shorthand.  The output printed above as "Out[3]" can be referred to as "%3".  Thus to divide Out[3] by Out[2] we simply write "%3/%2".  The moment of inertia of the cube is simply:

((M %3)/(%2))

(1) Below is the expression for calculating the mass of a rectangular solid (a parallelpiped)  height H, width W and length L.  Evaluate it.  Next, copy it, and alter the expression to calculate the moment of inertia of the rectangular solid.  Finally, express your answer in   terms of the mass of the object.

(∫_(-(L/2))^(L/2))(∫_(-(W/2))^(W/2))(∫_(-(H/2))^(H/2))  (Subscript[ρ, 0])  \[DifferentialD]z\[DifferentialD]y\[DifferentialD]x

Spheres:

How do w calculate the  calculate the mass or moment of inertia of  a sphere.? We can use the same trick as above.  We can break our sphere into a bunch of mini-cubes:

How do we specify that our mini-cubes add up to  a sphere and not a cube or a rectangular parallelpiped?  The trick here is understanding the limits of integration.  A sphere is described by the equation: 

                          (x^2)+(y^2)+(z^2) ≤ (R^2)
                           
We have integrals over x, y, and z, but they are not independent.  Let us imagine that we specify z first.  The variable z can range anywhere from -R to +R and still be inside the sphere.  Once we have specified z, we have limited our choice in y.  The variable y has a minimum value of -(Sqrt[(R^2)-(z^2)]) (which it can take if x=0).  It can take on a maximum  value of  (Sqrt[(R^2)-(z^2)]) (which again requires x=0).

   Once we have specified the value of y, we have specified the possible range of x.  If z and y are specified then x lies between -(Sqrt[(R^2)-(y^2)-(z^2)]) and +(Sqrt[(R^2)-(y^2)-(z^2)]).  If we choose any values of x outside this range, we would be outside the sphere.  Thus the  range of x, y and z, are determined by the equation of the sphere.

(2) Below are Mathematica expressions for the mass and moment of inertia of a sphere.  Evaluate them, and then show they produce the standard result for a sphere, 

                                     I=2/5 M R^2
                                     
There are two versions of the integral below.  The first versions may take too long to evaluate on your PC. (Different versions of Mathematica on different computers run slower or faster).  If you execute them and nothing happens after a few minutes (Mathematica keeps chugging along with no result)  then click on the darkened bracket, and then  select "Abort Evaluation" from the "Kernel" menu above.   Use the versions below that have the explicit function Integrate in them.

(Subscript[M, sphere])=(∫_(-R)^R)(∫_(-(Sqrt[(R^2)-(z^2)]))^(Sqrt[(R^2)-(z^2)]))(∫_(-(Sqrt[(R^2)-(z^2)-(y^2)]))^(Sqrt[(R^2)-(z^2)-(y^2)]))  (Subscript[ρ , 0])  \[DifferentialD]x\[DifferentialD]y\[DifferentialD]z

(Subscript[I, sphere])=(∫_(-R)^R)(∫_(-(Sqrt[(R^2)-(z^2)]))^(Sqrt[(R^2)-(z^2)]))(∫_(-(Sqrt[(R^2)-(z^2)-(y^2)]))^(Sqrt[(R^2)-(z^2)-(y^2)]))  (Subscript[ρ , 0]) ((x^2)+(y^2))  \[DifferentialD]x\[DifferentialD]y\[DifferentialD]z

Alternate version:  The "Assumptions" option in the second Integrate command tells  Mathematica to assume that R>z throughout the integration.   After all, Mathematica  doesn't know that we are trying to integrate over the inside of a sphere!  By telling it this, the integrals remain real (not complex) and everything proceeds smoothly.

(Subscript[M, sphere])=
	Integrate[
	 Integrate[
	    Integrate[(Subscript[ρ, 0  ]),{x, -(Sqrt[(R^2)-(z^2)-(y^2)]),  (Sqrt[(R^2)-(z^2)-(y^2)])}],
		                           {y,- (Sqrt[(R^2)-(z^2)]),(Sqrt[(R^2)-(z^2)])},    
		                                               Assumptions->{ ((R^2)-(z^2))>0}],
	                   {z,-R,R}]

(Subscript[I, sphere])=
	Integrate[
	 Integrate[
	    Integrate[(Subscript[ρ, 0])((x^2)+(y^2)),{x, -(Sqrt[(R^2)-(z^2)-(y^2)]),  (Sqrt[(R^2)-(z^2)-(y^2)])}],
		                           {y,- (Sqrt[(R^2)-(z^2)]),(Sqrt[(R^2)-(z^2)])},    
		                                               Assumptions->{ ((R^2)-(z^2))>0}],
	                   {z,-R,R}]

Ellipsoids:

Let's do something more interesting that a sphere.  We could have looked that answer up in the textbook!  Consider the case of a rotating ellipsoid.  The surface of a  general ellipsoid in 3D is described by the equation:

                                  ((((x/a)))^2)+((((y/b)))^2)+((((z/c)))^2)   = 1

which means that the points on the interior obey:

                                 ((((x/a)))^2)+((((y/b)))^2)+((((z/c)))^2)   ≤ 1

where  a, b and c are positive real numbers.  We'd like to evaluate the moment of inertia of an ellipsoid, too.  Here's a plot of what one looks like:

x=a Sin[θ] Cos[φ];
y=b Sin[θ] Sin[φ];
z=c Cos[θ];
a=1;
b=1;
c=3;
ParametricPlot3D[{x,y,z},{θ,0,π},{φ,0,2 π},BoxRatios->{a,b,c}]
Clear[x,y,z,a,b,c];

The Clear command is needed because below we will use the letters a, b, and c below as variables.  We have to tell Mathematica to "forget" the  values we assigned to them here in this calculation.  This can be a source of problems if you are not careful - you may at times forget that you have defined a quantity earlier on and be puzzled when Mathematica uses that definition at a later stage of calculation.

(3)  Show (on a separate piece of paper, by direct calculation) that the limits of integration on the  below expression are correct for evaluating the  total mass of the ellipsoid.  I used a different form for integration here - the Integrate command rather than "  ∫ ".   I had to do so because Mathematica got confused.  It needed to be told that certain constants were positive.  (That's the purpose of the Assumptions -> { ...}  option in the second Integrate command.)

Integrate[
	 Integrate[
	    Integrate[(Subscript[ρ, 0  ]),{x, -a (Sqrt[1-((((z/c)))^2)-((((y/b)))^2)]), a (Sqrt[1-((((z/c)))^2)-((((y/b)))^2)])}],
		                          {y,-b (Sqrt[1-((((z/c)))^2)]),b (Sqrt[1-((((z/c)))^2)])},    Assumptions->{b (Sqrt[1-((z^2)/(c^2))])>0&&(b^2) (1-((z^2)/(c^2)))>0}],
	                   {z,-c,c}]

Evaluate the above expression.  When you do, you will be a bit disappointed.  It is too messy!  It is messy because Mathematica is confused.  Mathematica does not know that Sqrt[b^2]=b.  (It doesn't know that b is a positive number).  In order to force it to simplify this expression, we can use the function PowerExpand[expression].   Insert the "%number" of the output above.

PowerExpand[]

Note that if  a=b=c, we recover the result for a sphere.

(4)  Below is the expression for the moment of inertia of the ellipsoid.  Evaluate it,  and force simplification via the PowerExpand command.   Determine the moment of inertia of an ellipsoid in the form of the total mass, M, times some function of a, b  and c.  Don't get confused - the big square root thingies are the limits of integration.  They are not part of the integrand.

(∫_(-c)^c)(∫_(-b (Sqrt[1-((((z/c)))^2)]))^(b (Sqrt[1-((((z/c)))^2)])))(∫_(-a (Sqrt[1-((((z/c)))^2)-((((y/b)))^2)]))^(a (Sqrt[1-((((z/c)))^2)-((((y/b)))^2)])))    (Subscript[ρ, 0])   ((x^2)+(y^2))\[DifferentialD]x\[DifferentialD]y\[DifferentialD]z

Other cases:

(5)  Use Mathematica to calculate the moment of inertia of a cube about one of its edges.
(Hint:  all you need to do is change the limits of integration in the integral we did above.)
Verify that this result agrees with that given by the parallel axis theorem.

(6)   Determine the moment of inertia of a pyramid of height L, with a square
LxL base, when rotated about a line passing through its apex and perpendicular to the base.
(See the figure below.)   You must use Mathematica to to do the integrals.

a=.9;
n=5;
Show[Graphics3D[
		Table[Cuboid[{i,j,-k},{i+a,j+a,-k+a}],{k,n,0,-1},{j,-k,k},{i,-k,k}]
		]]

If you were curious how I made the figures of the minicubes, the command for the sphere  and pyramid pictures are given below:

a=.8;
n=4;
Show[Graphics3D[
		Table[Cuboid[{i,j,k},{i+a,j+a,k+a}],{i,-n,n},{j,-Floor[Sqrt[n^2-i^2]],Floor[Sqrt[n^2-i^2]]},{k,-Floor[Sqrt[n^2-i^2-j^2]],Floor[Sqrt[n^2-i^2-j^2]]}]
		]]

a=.9;
n=5;
Show[Graphics3D[
		Table[Cuboid[{i,j,-k},{i+a,j+a,-k+a}],{k,n,0,-1},{j,-k,k},{i,-k,k}]
		]]