Game Physics Notes and Solutions

XY Particle

• Example 2.1: For a parabolic curve: r(t) = ti + t²j (Simpliy noted as: r=(t, t²)).
  
  \begin{equation} \begin{matrix} r(t) = ti + t^2j \\ \therefore{} \\ v(t) = \dot{r_t} = i + 2tj = (1, 2t) \\ \dot{s} = \vert{}v\vert{} = \sqrt{i^2 + (2tj)^2} = \sqrt{1 + 4t^2} \\ a = \dot{v} = \ddot{r} = 2j = (0, 2) \\ T(t) = \frac{v}{\vert{}v\vert{}} = \frac{i + 2tj}{\sqrt{i^2 + (2tj)^2}} = \frac{(1, 2t)}{\sqrt{1 + 4t^2}} = \frac{(1, 2t)}{\sqrt{i^2 + (2tj)^2}} \\ \kappa = \frac{v\cdot{}a^\perp}{|v|^3} = \frac{(1, 2t)\cdot{}(0, 2)^\perp}{\sqrt{1 + 4t^2}^3} = \frac{(1, 2t)\cdot{}(2, -0)}{\sqrt{1+4t^2}^3} = \frac{2 + 0}{\sqrt{1+4t^2}^3} = \frac{2}{\sqrt{1+4t^2}^3} \\ N = rotcounterclock(\frac{\pi}{2})T \\ \begin{bmatrix} cos\theta && -sin\theta \\ sin\theta && cos\theta \end{bmatrix} (@degree = 90^\circ) = \begin{bmatrix} 0 && -1 \\ 1 && 0 \\ \end{bmatrix}T \\ N(t) = \frac{(-2t, 1)}{\sqrt{i^2 + (2tj)^2}} \end{matrix} \end{equation}
• Planar Motion in Polar coordinates
  

  The position vector is represented as a unit length vector R := \(\frac{\textbf{r}}{\vert{}\textbf{r}\vert{}}\) and distance r=\(\vert{}\textbf{r}\vert{}\) and r = rR. Since R is unit length we write it as (cosθ, sinθ), θ ~ t. The unit length perpendicular vector P := (-sinθ, cosθ) obtained by \(\frac{\pi}{2}\) counter clockwise rotation of R in the plane. The moving frame is defined as {r(t); R(t), P(t)} is an alterative coordinate system. Int this system R and P need not be along the tangent and perpendicular to it. In this section there is an annoying double use of little r.

  \begin{equation} \begin{bmatrix} \dot{R} \\ \dot{P} \\ \end{bmatrix} = \begin{bmatrix} 0 && \dot{\theta} \\ -\dot{\theta} && 0 \end{bmatrix} \begin{bmatrix} R \\ P \\ \end{bmatrix} \end{equation}

  Or in otherwords: Ṙ = θ̇P and Ṗ = -θ̇R.

  This means that v = ṙ = \(\frac{d}{dt}\)(rR) = ṙR + rṘ = ṙR + rθ̇P.

  The acceleration then becomes: a = v̇ = \(\frac{d}{dt}\)(ṙR + rθ̇P) = r̈R + ṙṘ + ṙθ̇P + rθ̈P + rθ̇Ṗ = r̈R + ṙθ̇P + ṙθ̇P + rθ̈P - rθ̇θ̇R = (r̈ - rθ̇²)R + (2ṙθ̇ + rθ̈)P.

• Example 2.2: For spiral r = θ, θ ~ t.
  
  \begin{equation} \begin{matrix} \textbf{r} = \theta{}R = (\theta{}cos\theta{}, \theta{}sin\theta{}) \\ P = (-sin\theta{}, cos\theta{}) \\ \textbf{v} = \dot{\textbf{r}} = \dot{r}R + r\dot{\theta}P = \dot{\theta}R + \theta\dot{\theta}P = \dot{\theta}(cos\theta{}, sin\theta) + \dot{\theta}(-\theta{}sin\theta{}, \theta\cos\theta) = \dot{\theta}(cos\theta - \theta{}sin\theta{}, sin\theta + \theta{}cos\theta) \\ \textbf{a} = \dot{\textbf{v}} = (\ddot{\theta} - \theta{}\dot{\theta}^2)R + (2\dot{\theta}^2 + \theta\ddot{\theta})P \\ \end{matrix} \end{equation}
• Spatial motion in Cartesian coordinate
  

  In the (3d vector) spatial cooridantes we extend the planar case (i := (1, 0, 0), j := (0, 1, 0), k := (0, 0, 1)):

  \begin{equation} \begin{matrix} \textbf{r}(t) := x(t)i + y(t)j + z(t)k \\ \textbf{v}(t) := \dot{\textbf{r}} = \dot{x}(t)i + \dot{y}(t)j + \dot{z}(t)k \\ \textbf{a}(t) := \dot{\textbf{v}} = \ddot{x}(t)i + \ddot{y}(t)j + \ddot{z}(t)k \\ \end{matrix} \end{equation}

  Noting that the speed of the particle at time t is \(\vert{}\textbf{v}\vert\). If s(t) is the arclength measure along the curve, ṡ = |v|. The tangent vector stays the same T(t) = \(\frac{\textbf{v}}{\vert{}\textbf{v}\vert{}}\). But N cannot simply be the same as before as there are many (∞ on the circle normal to T, a slice plane) choices for perpendicular vector to T. We choose to stay in line with the 2D case using \(\frac{dT}{ds}\) which is \(\perp\) to T.

  \begin{equation} \begin{matrix} \frac{dT}{ds} =: \kappa(s)N(t) \end{matrix} \end{equation}

  This does not completly remove degenercy of choice as N has two varients: N and -N. If we take N κ must also be positve κ, but if -N is chosen -κ must be used to cancel it. Once we choose an N, κ is determined as well.

• Example 2.1: Curve of two pieces showing unsatisfiable N.
  
  \begin{equation} \begin{matrix} \textbf{r}(t) := (t, t^3, 0), t \lt 0 \\ \textbf{r}(t) := (t, 0, t^3), t \ge 0 \\ \end{matrix} \end{equation}

  We desire ?p := \(\textbf{r}\), \(\textbf{v}\), \(\textbf{a}\) are continuous at t = 0. ?p ≃ \(\lim_{t\to{}0}\textbf{r}(t) = \textbf{r}(0)\), \(\lim_{t\to{}0}\textbf{v}(t) = \textbf{v}(0)\), \(\lim_{t\to{}0}\textbf{a}(t) = \textbf{a}(0)\). Construct normal vectors N for each piece as a function of t. ?q := \(\lim_{t\to{}0^-}\textbf{N}(t)\) = (0, 1, 0), \(\lim_{t\to{}0^+}\textbf{N}(t)\) = (0, 0, 1). But since these one sided limits have differnet values N(t) is not continous at t=0. Changing sign of N cannot satify the equations.

  \begin{equation} \begin{matrix} !p := \\ r(t) = (t, t^3, 0) \wedge (t, 0, t^3), r(0) = (0, 0, 0) \wedge (0, 0, 0) \\ v(t) = (1, 3t^2, 0) \wedge (1, 0, 3t^2), v(0) = (1, 0, 0) \wedge (1, 0, 0) \\ a(t) = (0, 6t, 0) \wedge (0, 0, 6t^2), a(0) = (0, 0, 0) \wedge (0, 0, 0) \\ \end{matrix} \end{equation}

  Because ṡ lines up on both sides we get:

  \begin{equation} \begin{matrix} !q := \\ \dot{s} = \vert{}v\vert{} = \sqrt{1 + 9t^4} \\ T = \frac{(1, 3t^2, 0)}{\dot{s}} \wedge \frac{(1, 0, 3t^2)}{\dot{s}} \\ N = \frac{(-3t^2, 1, 0)}{\dot{s}} \wedge{} \frac{(-3t^2, 0, 1)}{\dot{s}} \\ \lim_{t\to{}0} = (0,1,0) \wedge (0,0,1) \to \bot \\ \end{matrix} \end{equation}
• Continue Spatial
  

  acceleration satifies a = s̈T + κṡ²N as before. Curvature however changes to:

  \begin{equation} \begin{matrix} \kappa := \sigma{}\frac{\vert{}\textbf{v}\times{}\textbf{a}\vert{}}{\vert{}\textbf{v}\vert{}^3} \end{matrix} \end{equation}

  where σ is a parameter chosen to be 1 or -1, if possible to make the normal vector continous. This means we can have a formula for the Normal vector as:

  \begin{equation} \begin{matrix} \textbf{N} := \frac{\sigma(\textbf{v}\times\textbf{a})\times\textbf{v}}{\vert{}\textbf{v}\times{}\textbf{a}\vert{}\vert{}\textbf{v}\vert{}} \\ \end{matrix} \end{equation}

  The tangent T and normal vector N only account for two of the three degrees of freedom in space. Question: Why not use Gram–Schmidt process for the basis? Well we use the cross product to define the binormal vector:

  \begin{equation} \begin{matrix} B = T \times N \end{matrix} \end{equation}

  Thus the coordiante system is {\(\textbf{v}\)(t);T(t);N(t);B(t)} as a moving frame of the curve. B is ⟂ to N and T and thus B ⋅ T = 0 = B ⋅ N ∀t. Also note as consequence: \(\frac{dB}{ds}\cdot{}T\) = 0 and \(\frac{dB}{ds}\cdot{}N\) = 0 (from the ⟂ equations differentiated). B is also unit (B⋅B = 1), meaning \(2\frac{dB}{ds}\cdot{}B\) = 0 = \(\frac{dB}{ds}\cdot{}B\). The derivative of B is ⟂ to T and B, making it ∥ to N.

  \begin{equation} \begin{matrix} \frac{dB}{ds} = -\tau{}N \end{matrix} \end{equation}

  For a scaler τ called the torsion of the curve (sign(-) is a common convetion). Curvature measures how much the curve wants to bend within the plane spanned by the tangent and noraml to the curve. The torsion measures the ammount the curve wants to bend out of the plane. Now lets complete the relationships with the cross product: N = B × T.

  \begin{equation} \begin{matrix} \frac{dN}{ds} = B\times\frac{dT}{ds} + \frac{dB}{ds}\times{}T = \kappa{}B\times{}N + -\tau{}N \times{} T = -\kappa{}T + \tau{}B \\ \end{matrix} \end{equation}

  These are the Frenet-Serret equations for the curve. They can be summerized in matrix notation:

  \begin{equation} \begin{bmatrix} \frac{dT}{ds} \\ \frac{dN}{ds} \\ \frac{dB}{ds} \\ \end{bmatrix} = \begin{bmatrix} 0 && \kappa && 0 \\ -\kappa && 0 && \tau \\ 0 && -\tau && 0 \end{bmatrix} \begin{bmatrix} T(s) \\ N(s) \\ B(s) \\ \end{bmatrix} \end{equation}

  For the remaining idiosyncrasies, using jerk (ȧ), v×a = κṡ³B, and v×a⋅ȧ = τκ²ṡ⁶ = τ|v×a|² we get τ:

  \begin{equation} \begin{matrix} \tau = \frac{v\times{}a\times{}\dot{a}}{\vert{}v\times{}a\vert{}^2} \\ \end{matrix} \end{equation}
• Spatial motion in cylindrical coodinates
  

  Conversion of (x,y,z) ~ Cartesian → Cylindrical (for r distance from origin = (0, 0, 0)):

  \begin{equation} \begin{matrix} x = rcos\theta \\ y = rsin\theta \\ z = z \\ \theta \in [0,2\pi) \\ \end{matrix} \end{equation}

  This gives us a unit vector R := (cosθ, sinθ, 0), ⟂ vector P := (-sinθ, cosθ, 0), and vertical unit k = (0, 0, 1). Moving frames are therfore { r (t), R (t), P (t), k }. The position of a point is:

  \begin{equation} \begin{matrix} \textbf{r} = r\textbf{R} + z\textbf{k} \\ \textbf{v} = \dot{\textbf{r}} = \dot{r}\textbf{R} + r\dot{\theta}\textbf{P} + \dot{z}\textbf{k} \\ \textbf{a} = \dot{\textbf{v}} = \ddot{r}\textbf{R} + \dot{r}\dot{\theta}\textbf{P} + \dot{r}\dot{\theta}\textbf{P} + r\ddot{\theta}\textbf{P} - r\dot{\theta}^2\textbf{R} + \ddot{z}\textbf{k} = (\ddot{r} - r\dot{\theta}^2)R + (2\dot{r}\dot{\theta} + r\ddot{\theta})P + \dot{z}\textbf{k} \\ \end{matrix} \end{equation}

  The time derivative then can be sumerized in matrix notation:

  \begin{equation} \begin{bmatrix} \dot{R} \\ \dot{P} \\ \dot{k} \\ \end{bmatrix} = \begin{bmatrix} 0 && \dot{\theta} && 0 \\ -\dot{\theta} && 0 && 0 \\ 0 && 0 && 0 \\ \end{bmatrix} \begin{bmatrix} R \\ P \\ k \\ \end{bmatrix} \end{equation}
• Exercise 2.2: For the helix (cost, sint, t) get position, velocity and acceleration vectors.   TypicallyHidden
  
  \begin{equation} \begin{matrix} ?\textbf{r}, ?\textbf{v}, ?\textbf{a} \\ Cartesian(cos(t), sin(t), t) \to Cylindrical(?r, ?\theta, ?z) \\ cos(t) = x = ?rcos?\theta \\ sin(t) = y = ?rsin?\theta \\ t = z = !z \\ \vert{}(x,y,z)\vert{} = \sqrt{cos(t)^2 + sin(t)^2 + t^2} = \sqrt{1 + t^2} \\ !r = \sqrt{cos(t)^2 + sin(t)^2} = \sqrt{1} = 1 \\ !\theta = tan^{-1}(\frac{sin(t)}{cos(t)}) \approx t \\ !\textbf{r} = rR + tk = (1, 0, t) \\ \dot{r} = 0, \dot{\theta} = 1, \dot{z} = 1 \\ !\textbf{v} = \dot{r}\textbf{R} + r\dot{\theta}\textbf{P} + \dot{z}k = 0 + 1\textbf{P} + k = (0, 1, 1) \\ \ddot{r} = 0, \ddot{\theta} = 0, \ddot{z} = 0 \\ !\textbf{a} = (\ddot{r} - r\dot{\theta}^2)R + (r\ddot{\theta} + 2\dot{r}\dot{\theta})P + \ddot{z}k = (0 - 1*1)R + (0 + 0)P + 0 = (-1, 0, 0) \\ \end{matrix} \end{equation}
• Spatial Motion in Spherical Coordinates
  

  Here a position point, r, is

  \begin{equation} \begin{matrix} \textbf{r} = \rho{}\textbf{R} \\ \textbf{R} = (cos\theta{}sin\varphi, sin\theta{}sin\varphi, cos\varphi)\ \sim\ unit \\ \textbf{P} = (-sin\theta, cos\theta, 0) \\ \textbf{Q} = \textbf{R}\times{}\textbf{P} = (-cos\theta{}cos\varphi, -sin\theta{}cos\varphi, sin\varphi) \\ \end{matrix} \end{equation}

  Making a moving frame {r(t);R(t);Q(t);R(t)}. The derivative relation is provided as (skew symmetric):

  \begin{equation} \begin{bmatrix} \dot{\textbf{P}} \\ \dot{\textbf{Q}} \\ \dot{\textbf{R}} \\ \end{bmatrix} = \begin{bmatrix} 0 && \dot{\theta}cos\varphi && -\dot{\theta}sin\varphi \\ -\dot{\theta}cos\varphi && 0 && \dot{\varphi} \\ \dot{\theta}sin\varphi && -\dot{\varphi} && 0 \\ \end{bmatrix} \begin{bmatrix} \textbf{P} \\ \textbf{Q} \\ \textbf{R} \\ \end{bmatrix} \end{equation} \begin{equation} \begin{matrix} \textbf{v} = \dot{\rho}\textbf{R} + \rho\dot{\textbf{R}} \\ \end{matrix} \end{equation}
• TODO Example 2.3: Spherical Helix (cos(t), sin(t), t)/\(\sqrt{1+t^2}\)
  
  \begin{equation} \begin{matrix} (x, y, z) = \frac{(cost, sint, t)}{\sqrt{1 + t^2}} \\ Cartesian(x,y,z) \to Spherical(?\rho, ?\theta, ?\varphi) \\ \vert{}(x,y,z)\vert{} = \sqrt{(\frac{cos(t)}{\sqrt{1 + t^2}})^2 + (\frac{sin(t)}{\sqrt{1 + t^2}})^2 + (\frac{t}{\sqrt{1 + t^2}})^2} = \sqrt{(\frac{1}{\sqrt{1 + t^2}})^2(cos(t)^2 + sin^2(t)) + (\frac{t}{\sqrt{1 + t^2}})^2} = \sqrt{(\frac{1}{\sqrt{1 + t^2}})^2 + (\frac{t}{\sqrt{1 + t^2}})^2} \& \\ \& \sqrt{\frac{1}{1 + t^2} + \frac{t^2}{1 + t^2}} = \sqrt{\frac{1 + t^2}{1 + t^2}} = 1 \\ \vert{}(x,y,z)\vert{} = 1 \therefore{}(Cartesian\to{}Spherical)(!\rho{}=1). \\ \frac{cos(t)}{\sqrt{1 + t^2}} = x = cos?\theta{}sin?\varphi) \\ \frac{sin(t)}{\sqrt{1 + t^2}} = y = sin?\theta{}sin?\varphi) \\ \frac{t}{\sqrt{1 + t^2}} = z = cos?\varphi \\ \therefore{} cos^2?\varphi{} + sin^2?\varphi = 1 \& \\ \& (\frac{t}{\sqrt{1 + t^2}})^2 + sin^2?\varphi = 1 \& \\ \& \frac{t^2}{1 + t^2} + sin^2?\varphi = 1 \& \\ \& sin^2?\varphi = 1 - \frac{t^2}{1 + t^2} \& \\ \& sin^2?\varphi = \frac{1 + t^2 - t^2}{1 + t^2} \& \\ \& sin?\varphi = \sqrt{\frac{1}{1 + t^2}} = \frac{1}{\sqrt{1 + t^2}} \\ \frac{cos(t)}{\sqrt{1 + t^2}} = x = cos(?\theta{})\frac{1}{\sqrt{1 + t^2}} \\ \frac{sin(t)}{\sqrt{1 + t^2}} = y = sin(?\theta{})\frac{1}{\sqrt{1 + t^2}} \\ \therefore{} !\theta = t \\ !\varphi = arccos(\frac{z}{\sqrt{x^2 + y^2 + z^2}}) = arccos(\frac{\frac{t}{\sqrt{1 + t^2}}}{1}) \\ \therefore{} (x,y,z) \to (\rho = 1, \theta = t, \varphi = arccos(\frac{t}{\sqrt{1 + t^2}})) \\ Cont. \end{matrix} \end{equation}
• TODO Example 2.4 page 24.   TypicallyHidden
  
• Motion about a fixed axis
  

  If we have a particle rotating about a fixed axis at a constant distance, we wish to find the position, velocity, and acceleration. Assume the axis is a line that contains the origin and unit length direction is D. We choose this in order to have D play the role of k and \(\textbf{R}\)(t) = (cosθ(t), sinθ(t), 0) is radial to that axis. Our axis are Cylindrical with vectors (ξ,η,\(\textbf{D}\)) (\(\textbf{R}\) = cosθξ + sinθη, \(\textbf{P}\) = cosθη - sinθξ, \(\textbf{k}\) = D) (ξ η are the spinning directions). Angular speed is σ(t) = θ̇(t) (or θ̇ = σ), angular velocity is w(t) = σ(t)\(\textbf{D}\), and angular accleration is α(t) = σ̇(t)\(\textbf{D}\). If (r = r₀) is the constant distance from the particle to axis and (z = h₀) is the costant height above the plane D ⋅ r = 0, then:

  \begin{equation} \begin{matrix} \textbf{r}(t) = r_0\textbf{R}(t) + h_0\textbf{D} = r_0(cos\theta{}\xi + sin\theta{}\eta)(t) + h_0\textbf{D} \\ \textbf{r}(0) = r_0\xi + h_0\textbf{D} \\ \dot{r} = \dot{r_0} = 0, \dot{z} = \dot{h_0} = 0 \\ \textbf{v}(t) = \dot{\textbf{r}}(t) = \dot{r}\textbf{R} + r\dot{\theta}\textbf{P} + \dot{z}\textbf{k} = r_0\sigma\textbf{P} = r_0\sigma\textbf{D}\times{}\textbf{R} = r_0\sigma{}\textbf{w}\times{}\textbf{r} \\ \textbf{a}(t) = (\ddot{r_0} - r_0\dot{\sigma}^2)R + (r\ddot{\sigma} + 2\dot{r_0}\dot{\theta})P + \ddot{h_0}k = (-r_0\dot{\sigma}^2,r\ddot{\theta},0) = -r_0\sigma^2R + \alpha \times\textbf{r} \\ \end{matrix} \end{equation}

  Where -r₀σ²R is the centripetal acceleration of the particle, α×\(\textbf{r}\) is the tangental acceleration.

• Motion about a Moving Axis
  

  Lets define Skew on a vector (u₁, u₂, u₃).

  \begin{equation} Skew(u) = \begin{bmatrix} 0 && -u_3 && u_2 \\ u_3 && 0 && -u_1 \\ -u_2 && u_1 && 0 \\ \end{bmatrix} \end{equation}

  Note that this matrix has property \(Skew(\textbf{u})\textbf{r} = \textbf{u}\times{}\textbf{r}\) - page(26) for the rest of the argument as it relies on a proof from chapter 10.

Particle Systems and Continuous Materials

First law

Second law

Third law

Closing

Gravitational Forces

Spring Forces

Friction and Dissipative Forces

• Friction
  

  Two objects in contact sliding across each other experience friction. In the case of a moving object sliding over a non-moving one the frictional force is tangent to the surface of the adjacent (non-moving) and opposite in direction to the moving velocity. The magnitude is ∝ the normal force between the two objects (also assumed to be independant of area of contact? and intependant to speed of objcet once that object starts to move). These assumptions argue n = 0. So the Friction becomes:

  \begin{equation} F = \begin{matrix} -c_k\frac{v}{\vert{}v\vert{}} && v \ne 0 \\ 0 && v = 0 \\ \end{matrix} \end{equation}

  cₖ is the coefficient of kinetic friction. The coefficient is the ratio of the magnitude’s of the frictional force over the normal force, cₖ = \(\frac{F_{friction}}{F_{normal}}\), but with correct units so that RHS has units force. The transition between the ≠0 and =0 force cases introduces the concept of static friction. Small initial magnitudes are not enough to overcome the friction and objects wont move when experiencing a force (when overcoming it becomes kinetic friction). The static friction coefficient is cₛ = \(\frac{max(F_{friction})}{F_{normal}}\).

• Viscosity
  

  The dissipative force in general is modeled for viscosity at n = 1.

Torque

Equilibrium

Momenta

• Linear Momentum:
  
  \begin{equation} \begin{matrix} \mathbf{p} = m\mathbf{v} \\ \end{matrix} \end{equation}

  Applied force is related to momentum by \(F = \frac{dp}{dt}\). In the p-particle case momentum becomes:

  \begin{equation} \begin{matrix} \mathbf{p} = \sum_{i=1..p}m_i\mathbf{v_i}\\ \end{matrix} \end{equation}

  And the continuum case becomes (occupying a region, R):

  \begin{equation} \begin{matrix} \mathbf{p} = \int_{R}\mathbf{v}dm = \int_R\delta\mathbf{v}dR\\ \end{matrix} \end{equation}

  δ is the mass density (dm = δdR is an infinitesimal measure of mass). Both δ and v may depend on spatial location and thus cannot be factored outside the ∫. These quantities obey the law of conservation of linear momentum, which states that if the net external force is 0, then so too is the change in linear momentum (i.e. the momentum is constant).

• Angular Momentum
  

  Linear momentum encapsulates the intuition that moving along a straight path experiencing no new forces should remain moving. Angular Momentum is a little less intuitive but follows the same philosophy. This is the tendency for an object to keep rotating about an axis. For a point mass, m, the angular momentum about the origin for th particle is

  \begin{equation} \begin{matrix} \mathbf{L} = \mathbf{r}\times\mathbf{p} = \mathbf{r}\times{}m\mathbf{v} \end{matrix} \end{equation}

  Where r is the vector from origin to the particle. L is ⟂ to r and \(\mathbf{v}\). Similarly to before if we have a system of particles p : Particle(n). The angular momentum is:

  \begin{equation} \begin{matrix} \mathbf{L} = \sum_{i..n}\mathbf{r_i}\times{}m_i\mathbf{v_i} \end{matrix} \end{equation}

  and for a continuum occupying region, R.

  \begin{equation} \begin{matrix} \mathbf{L} = \int_{R}\mathbf{r}\times{}\mathbf{v}dm = \int_{R}\delta\mathbf{r}\times{}\mathbf{v}dR \end{matrix} \end{equation}

  As force is the derivative of linear momentum, torque is the derivative of angular momentum. There also applies the law of conservation of angular momentum which states that the net external torque on a system is 0, the angular momentum (\(\textbf{L}\)) is constant.

Center of Mass

• Discrete mass in 1D
  
  • m₁ : mass
  • m₂ : mass
  • x₁ : 1DPoint
  • x₂ : 1DPoint
  • x̄ : 1DPoint ~ x₁ ≤ x̄ ≤ x₂

  Assume the line connecting m₁ and m₂ is negligible mass. Gravity is assumed to be ’downwards’. This is the setup for a seasaw. We wish to balance this system. For example if m₁ =ₘ m₂ is inhabited then x̄ = \(\frac{x_1 + x_2}{2}\) will balance the system. The force on mass mᵢ is mᵢg and the torque to x̄ is mᵢg(xᵢ - x̄), from r*F*sin(θ). Lets balance these forces by setting them to 0.

  \begin{equation} \begin{matrix} 0 = m_1g(x_1 - \bar{x}) + m_2g(x_2 - \bar{x}) = g[m_1x_1 + m_2x_2 - (m_1 + m_2)\bar{x}] \\ \bar{x} = \frac{m1x_1 + m_2x_2}{m_1 + m_2} \\ \end{matrix} \end{equation}

  This tells us the center of mass for this system. We can similarly extend to pm : zip(mass(p), 1DPoints(p)):

  \begin{equation} \begin{matrix} \sum_{p\in{}pm}^{}m(p)g(x(p) - \bar{x}) = 0 \\ \bar{x} = \sum_{p\in{}pm}^{}\frac{m(p)}{\sum_{q\in{}pm}^{}m(q)} \\ \end{matrix} \end{equation}

  Where the sum of all the masses is the total mass and the sum(mᵢxᵢ) is the moment of the system about the origin.

• Continuous mass 1D
  

  Take the endpoints of our continuum (wire) as a and b, a < b. Since the mass can vary over the interval, we use the mass density function δ(x) and the relation dm = δ(x)dx. The last thing is gravity, its infintesmal force is then gdm = δ(x)gdx and the infinitesimal torque about x̄ is (x-x̄)gdm. For x̄ to be the center of mass its torque must be 0.

  \begin{equation} \begin{matrix} \int_{a}^{b}(x-\bar{x})gdm = g\int_{a}^{b}(x-\bar{x})\delta(x)dx = 0 \\ \end{matrix} \end{equation}

  Which we solve for the center of mass at:

  \begin{equation} \begin{matrix} \bar{x} = \frac{\int_{a}^{b}x\delta{}(x)dx}{\int_a^b\delta(x)dx} \\ \bar{x} = \frac{b - a}{2}, \delta = c : Number \\ \end{matrix} \end{equation}

  The denominator is the total mass and the numerator is the moment of the system about the origin. The second equation is the case where the mass density function is independent of x.

• Discrete mass in 2D
  

  Lets extend our idea to 2D (quite naturally), start but upgrading (denoting zip as ×ₚ, ₚ being a 1..p-element-wise pair) (pm : Mass(p) ×ₚ 1DPoint(p) = {m, x}) → (pm : Mass(p) ×ₚ 2DPoint(p) = {m, x, y}). We now wish to find the center of mass at (x̄, ȳ). The gravitation force on each mass is mᵢg and the torque about (x̄, ȳ) is mᵢg(xᵢ - x̄, yᵢ - ȳ), which we try and have be 0⃗ for equilibrium.

  \begin{equation} \begin{matrix} \sum_{p∈pm}m(p)g(x(p) - \bar{x}, y(p) - \bar{y}) = (0, 0) \\ \end{matrix} \end{equation}

  Which we unsurprisingly get:

  \begin{equation} \begin{matrix} (\bar{x}, \bar{y}) = (\frac{\sum_{p\in{}pm}m(p)x(p)}{\sum_{p\in{}pm}m(p)}, \frac{\sum_{p\in{}pm}m(p)y(p)}{\sum_{p\in{}pm}m(p)}) \end{matrix} \end{equation}

  Now we get \(M_y\) the sum(m(p)x(p)) is the moment about the y-axis and Mₓ the sum(m(p)y(p)) is the moment about the x-axis. In this case our connecting line has turned into a disc.

• Continuous mass 2D
  

  Given a bounded region R that contains a plane of continous mass. Each point in the region as an associated infinitesimal mass: dm. The mass over an infinitesimal rectangle in x and y is then dA = dxdy, finally with a mass density function depending on position we get: dm = δ(x,y)dxdy. The infinitesimal torque to a location (x̄, ȳ) is (x - x̄, y - ȳ)gdm, making equilibrium achieved at:

  \begin{equation} \begin{matrix} \int\int_R(x-\bar{x}, y-\bar{y})gdm = 0 \\ (\bar{x}, \bar{y}) = \frac{\int\int_R(x,y)\delta(x,y)dxdy}{\int\int_R\delta(x,y)dxdy} = (\frac{\int\int_Rx\delta(x,y)dxdy}{\int\int_R\delta(x,y)dxdy},\frac{\int\int_Ry\delta(x,y)dxdy}{\int\int_R\delta(x,y)dxdy}) \end{matrix} \end{equation}

  The ∫∫δ is the total mass of the system, ∫∫xδdxdy is the moment of the system about the y-axis, ∫∫yδdxdy is the moment of system about the x-axis.

  • Example 2.3 - Region(y = x² ∧ y = 1), assume δ(x,y) = 1.
    

    The functions intersect at y, x² = 1 ∴ x ranges ±1. We want in y which curve is on top, to define the ranging. y = x² is below y = 1, So y ranges from x² to 1.

    \begin{equation} \begin{matrix} m = \int\int_Rdydx = \int_{-1}^1\int_{x^2}^1dxdy = \int_{-1}^1(1 - x^2)dx = \frac{4}{3} \\ M_y = \int\int_Rxdxdy = \int_{-1}^1\int_{x^2}^1xdxdy = \int_{-1}^1x(1 - x^2)dx = \frac{x^2}{2} - \frac{x^4}{4}\Big|_{-1}^1 = 0 \\ M_x = \int\int_Rydxdy = \int_{-1}^1\int_{x^2}^1ydxdy = \int_{-1}^1\frac{1 - x^4}{2}dx = \frac{4}{5} \\ (\bar{x}, \bar{y}) = \frac{(M_y, M_x)}{m} = (0, \frac{3}{5}) \\ \end{matrix} \end{equation}
  • Mass distributed along a curve, such as mass varying arclength planar wire
    

    Assume the curve is continously differentiable and that it is specified parametrically by (x(t),y(t)) for t ∈ [a, b]. In terms of arc length s, the mass density is δ̄(s) (in terms of curve paramter its δ(t)). Infinitesimal mass at the postion correspoding to t is distributed over infinitesimal arc length ds. Thusly, dm = δ(t)ds, ds = √(ẋ² + ẏ²) (dot derivative W.R.T t). With L = total length of the curve.

    \begin{equation} \begin{matrix} m = \int_0^L\bar{\delta}(s)ds = \int_a^b\delta(t)\sqrt{\dot{x}^2+\dot{y}^2}dt \\ M_y = \int_0^Lx\bar{\delta}(s)ds = \int_a^bx(t)\delta(t)\sqrt{\dot{x}^2+\dot{y}^2}dt \\ M_x = \int_0^Ly\bar{\delta}(s)ds = \int_a^by(t)\delta(t)\sqrt{\dot{x}^2+\dot{y}^2}dt \\ (\bar{x}, \bar{y}) = \frac{(M_y, M_x)}{m} \\ \end{matrix} \end{equation}
  • Example 2.4: Center of mass for a wire, δ(x,y) = 1, Region(x² + y² = 1 ∧ y ≥ 0)
    

    Here as the circle is effectly cut in half at the x-axis, x ranges from (-1, 1) and y ranges from (0, 1).

    x = var('x')
    y = var('y')
    eq1 = x^2 + y^2 == 1
    eq2 = y >= 0
    #; solve([eq1, eq2], x, y)
    #; The curve is a circle so we parameterize as (cos(t), sin(t))
    t = var('t')
    v = vector([cos(t), sin(t)]) # over [0,π]
    vdt = diff(v, t) # (-sin(t), cos(t))
    m = integral(sqrt(vdt[0]^2 + vdt[1]^2), t, 0, pi) #; π
    Mx = integral(v[1]*sqrt(vdt[0]^2 + vdt[1]^2), t, 0, pi) #; 2
    My = integral(v[0]*sqrt(vdt[0]^2 + vdt[1]^2), t, 0, pi) #; 0
    CoM = vector([My, Mx]) / m #; (0, 2/π)
    

• Discrete mass in 3D
  

  Now just like before we promote (pm : Mass(p) ×ₚ Point2D(p) = {m, x, y}) → (pm : Mass(p) ×ₚ Point3D(p) = {m, x, y, z}). The logic continues, and the torque at (x̄, ȳ, z̄) is now: mᵢg(xᵢ-x̄, yᵢ-ȳ, zᵢ-z̄) which we want equal to 0⃗.

  \begin{equation} \begin{matrix} \sum_{p∈pm}m(p)g(x(p) - \bar{x}, y(p) - \bar{y}, z(p) - \bar{z}) = (0, 0, 0) \\ (\bar{x}, \bar{y}, \bar{z}) = (\frac{\sum_{p\in{}pm}m(p)x(p)}{\sum_{p\in{}pm}m(p)}, \frac{\sum_{p\in{}pm}m(p)y(p)}{\sum_{p\in{}pm}m(p)}, \frac{\sum_{p\in{}pm}m(p)z(p)}{\sum_{p\in{}pm}m(p)}) \\ \end{matrix} \end{equation}

  Now our moments are \(M_{yz}\) = sum(mᵢxᵢ) (yz-plane moment), \(M_{xz}\) = sum(mᵢyᵢ) (xz-plane moment), and \(M_{xy}\) = sum(mᵢzᵢ) (xy-plane moment).

• Continuous mass in 3D
  

  Now instead of a region, we work with a bounded volume V, with infinitesimal mass dm at (x,y,z) being a cube dxdydz = dV. The mass density function is also upgraded to δ(x,y,z). The total torque (which we want to be 0⃗) is:

  \begin{equation} \begin{matrix} \int\int\int_V(x - \bar{x}, y - \bar{y}, z - \bar{z})gdm = 0 \\ \end{matrix} \end{equation}
  • Example 2.5
    

    Compute the center of mass of the solid hemisphere x² + y² + z² ≤ 1, with z ≥ 0, and δ(x,y,z) ≡ 1.

    We want to calculate the volume so we use the volume of a sphere formula. Now remember this is a full sphere and we want the hemisphere. I get a bit confused sometimes as I see the first and second equations for volume of sphere, I use the second here. The ρ²sin(<angle>) should work if <angle> is the polar angle.

    \begin{equation} \begin{matrix} V = \rho^2sin(\theta)d\theta{}d\varphi{}d\rho{} \\ dV = \rho^2sin(\varphi)d\theta{}d\varphi{}d\rho{} \\ \int_0^R\int_0^{2\pi}\int_0^\pi\rho^2sin(\varphi)d\theta{}d\varphi{}d\rho{} \\ \end{matrix} \end{equation}

    x, y, z = var('x, y, z')
    eq1 = x^2 + y^2 + z^2 <= 1
    eq2 = z >= 0
    #; But we want to work with a sphere so spherical coords
    ρ = var('ρ') #; ρ² ≤ 1 ∴ ρ ∈ [0,1]
    θ = var('θ') #; θ ∈ [0, 2π]
    φ = var('φ') #; φ is cut in half (upper half is 0 to π/2, non-negative values) so ∈ [0, π/2]
    eqx = ρ*sin(φ)*cos(θ)
    eqy = ρ*sin(φ)*sin(θ)
    eqz = ρ*cos(φ)
    eqm = ρ^2*sin(φ)
    m = integral(integral(integral(eqm, θ, 0, 2*π), ρ, 0, 1), φ, 0, π/2) #; 2/3*π
    Mxy = integral(integral(integral(eqz*eqm, θ, 0, 2*π), ρ, 0, 1), φ, 0, π/2) #; π/4 
    Mxz = integral(integral(integral(eqy*eqm, θ, 0, 2*π), ρ, 0, 1), φ, 0, π/2) #; 0
    Myz = integral(integral(integral(eqx*eqm, θ, 0, 2*π), ρ, 0, 1), φ, 0, π/2) #; 0
    z̄ = Mxy / m
    

    From this we know (x̄, ȳ, z̄) = \(\frac{(M_{yz}, M_{xz}, M_{xy})}{m}\) = (0, 0, \(\frac{3}{8}\)).

  • Surface mass
    

    In the case of a bounded surface S, the infinitesimal mass dm at (x,y,z) is distributed in an infinitesimal surface area dS. This of-course depends on the definition of the surface, for example if the surface is z = f(x,y) then:

    \begin{equation} \begin{matrix} dS = \sqrt{1 + (\frac{\partial{}f}{\partial{}x})^2 + (\frac{\partial{}f}{\partial{}y})^2}dxdy \end{matrix} \end{equation}

    If it is instead defined parametrically: P(u, v) = (x(u, v), y(u, v), z(u, v)) then:

    \begin{equation} dS = \begin{vmatrix} \frac{\partial{}P}{\partial{}u} \times{} \frac{\partial{}P}{\partial{}v} \end{vmatrix} dudv \end{equation}

    The density of the distrobution is δ a function of (x,y,z), then center of mass is then:

    \begin{equation} \begin{matrix} (\bar{x},\bar{y},\bar{z}) = \frac{\int\int_S\delta(x,y,z)dS}{\int\int_SdS} \end{matrix} \end{equation}

    Where the denominator is the total mass and the numerators contains the moments as before (they use exclusionary names remember).

  • Example 2.6
    

    Center of mass of hemishpere: x² + y² + z² = 1, with z ≥ 0, δ ≡ 1.

    x, y, z = var('x y z')
    eq1 = x^2 + y^2 + z^2 == 1
    eq2 = z >= 0
    ρ = 1 #; From eq1
    θ, φ = var('θ φ') #; φ ranges [0,π/2] from eq2
    eqm = ρ^2*sin(φ)
    m = integral(integral(eqm, θ, 0, 2*π), φ, 0, π/2) #; ρ is 1 so left out, 2π
    eqyz_x = ρ*sin(θ)*cos(φ)
    eqxz_y = ρ*sin(θ)*sin(φ)
    eqxy_z = ρ*cos(φ)
    M_yz = integral(integral(eqm * eqyz_x, θ, 0, 2*π), φ, 0, π/2) #; 0
    M_xz = integral(integral(eqm * eqxz_y, θ, 0, 2*π), φ, 0, π/2) #; 0
    M_xy = integral(integral(eqm * eqxy_z, θ, 0, 2*π), φ, 0, π/2) #; π
    (x̄, ȳ, z̄) = result = vector([M_yz, M_xz, M_xy]) / m #; (0, 0, 1/2)
    

  • Curve Mass
    

    If we have a continuous material that is a wire along a curve in 3D. Assume this curve is continually differentiable and parametrically defined: (x(t),y(t),z(t)) for t ∈ [a,b]. For arc length, s, the mass is δ̄(s) and in terms of curve parameter, t, it is δ(t). The infinitesimal mass is over infinitesimal arc length ds, dm = δ(t)ds with ds = \(\sqrt{\dot{x}^2+\dot{y}^2+\dot{z}^2}dt\) for parametric curves in space (dot for derivative respecting t). Total mass is given as:

    \begin{equation} \begin{matrix} m = \int_0^L\bar{\delta}(s)ds = \int_a^b\delta(t)\sqrt{\dot{x}^2+\dot{y}^2+\dot{z}^2}dt \\ \end{matrix} \end{equation}

    Each moment is substituted (with exculsionary names α) as:

    \begin{equation} \begin{matrix} M_{\neg\alpha} = \int_0^L\alpha\bar{\delta}(s)ds = \int_a^b\alpha(t)\delta(t)\sqrt{\dot{x}^2+\dot{y}^2+\dot{z}^2}dt \\ \end{matrix} \end{equation}

    And as before the center of mass is the moments divided by the total mass.

  • Example 2.7
    

    Compute center of mass of constant density (δ ≡ 1) wire in the shape of a helix: (x(t), y(t), z(t)) = (cos(t),sin(t),t²/2), t∈[0,2π].

    t = var('t') #; ∈ [0, 2π]
    f_x = cos(t)
    f_y = sin(t)
    f_z = t^2 / 2
    eqm = sqrt(diff(f_x, t)^2 + diff(f_y, t)^2 + diff(f_z, t)^2)
    m = integral(eqm, t, 0, 2*π) #; ≈ 21.2562941482091
    M_yz = numerical_approx(integral(eqm * f_x, t, 0, 2*π))
    M_xz = numerical_approx(integral(eqm * f_y, t, 0, 2*π))
    M_xy = numerical_approx(integral(eqm * f_z, t, 0, 2*π))
    (x̄, ȳ, z̄) = result = vector([M_yz, M_xz, M_xy]) / numerical_approx(m) 
    

    From this we get the mass is: \(\pi \sqrt{4 \, \pi^{2} + 1} + \frac{1}{2} \, \operatorname{arsinh}\left(2 \, \pi\right)\) ≈ 21.2562941482091. And the center of mass as: (0.01820556170283781, -0.27415670682952065, 9.390631601721266).

• Moments and Products of Inertia
  

  The moment of inertia is the measure of the rotational inertia of a body about an axis. The more difficult to set rotation about an axis the stronger the moment is about that axis (intuitively, if the mass is large or distance from axis is large). For a single particle emperical results show mr² (r ~ distance to axis) as the moment.

  • Moments of inertia in 1D
    

    With Particle1D := Mass(p) ×ₚ 1DPoint(p). If we have ps : Particle1D(p), total mass is \(m = \sum_{p\in{}ps}m(p)\). We then define a useful quantity for moment about x=0 (origin)

    \begin{equation} \begin{matrix} M_0 = \sum_{p\in{}ps}m(p)x(p) \end{matrix} \end{equation}

    These equations are special cases of

    \begin{equation} \begin{matrix} \sum_{p\in{}ps}m(p)x(p)^k \end{matrix} \end{equation}

    k = 0 for mass k = 1 for moment. Now what is k = 2? This is the case of moment of inertia with respect to the origin, I₀, of a real line. The moment of inertia with respect to the center of mass (x̄) is then:

    \begin{equation} \begin{matrix} I = \sum_{p\in{}ps}m(p)(x(p) - \bar{x})^2 = I_0 - m\bar{x}^2 \end{matrix} \end{equation}

    This translates simply for a continuous mass (with density function δ over [a,b]) as total mass m = ∫ₐᵇδ(x)dx and M₀ = ∫ₐᵇxδ(x)dx, x̄ = \(\frac{M_0}{m}\). The moment of inertia becomes (with respect to origin then center of mass):

    \begin{equation} \begin{matrix} I_0 = \int_a^bx^2\delta(x)dx \\ I = \int_a^b(x - \bar{x})^2\delta(x)dx = I_0 - m\bar{x}^2 \\ \end{matrix} \end{equation}
  • Moments of inertia in 2D
    

    With Particle2D := Mass(p) ×ₚ 2DPoint(p). Assume we have ps : Particle2D(p). Avoiding obvious restating, we define the moment of inertia with respect to the origin and with respect to center of mass (x̄, ȳ) as:

    \begin{equation} \begin{matrix} I_0 = \sum_{p\in{}ps}m(p)(x(p)^2 + y(p)^2) \\ I = \sum_{p\in{}ps}m(p)\vert{}(x(p), y(p)) - (\bar{x}, \bar{y})\vert{}^2 = I_0 - m(\bar{x}^2 + \bar{y}^2). \\ \end{matrix} \end{equation}

    For a continous mass function δ(x,y) over region, R.

    \begin{equation} \begin{matrix} I_0 = \int\int_R(x^2 + y^2)\delta(x,y)dxdy \\ I = \int\int_R((x - \bar{x})^2 + (y - \bar{y}^2))\delta(x,y)dxdy = I_0 - m(\bar{x}^2 + \bar{y}^2)\\ \end{matrix} \end{equation}
  • Moment of inertia in 3D
    

    Here we will mix up the nomenclature by not talking about inertia about a point but instead a line. Given an origin \(\mathcal{O}\), Particle3D := Mass(p) ×ₚ 3DPoint(p), ps : Particle3D(p) and ∀3DPoint(x,y,z) ≃ ∃r (For any 3DPoint we have a ’r’ we could refer to the whole point as relative to \(\mathcal{O}\)). The moments of inertia are then defined as (about x, y, and z axes):

    \begin{equation} \begin{matrix} I_{xx} = \sum_{p∈ps}m(p)(y(p)^2 + z(p)^2) && I_{yy} = \sum_{p∈ps}m(p)(x(p)^2 + z(p)^2) && I_{zz} = \sum_{p∈ps}m(p)(x(p)^2 + y(p)^2) \\ \end{matrix} \end{equation}

    For the moment of inertia of a line, L, through \(\mathcal{O}\) given parametrically as \(\mathcal{O}\) + \(t\mathbf{D}\) with unit length direction vector \(\mathbf{D}\) = (d₁, d₂, d₃) is the sum of m(p)r(p)² where r is distance from r(p) to L.

    \begin{equation} \begin{matrix} I_L = \sum_{p\in{}ps}m(p)(\vert{}r(p)\vert{}^2 - (\mathbf{D}\cdot{}r(p))^2) = \\ d_1^2\sum_{p\in{}ps}m(p)(y(p)^2 + z(p)^2) + d_2^2\sum_{p\in{}ps}m(p)(x(p)^2 + z(p)^2) + d_3^2\sum_{p\in{}ps}m(p)(x(p)^2 + y(p)^2) \\ - 2d_1d_2\sum_{p\in{}ps}m(p)x(p)y(p) - 2d_1d_3\sum_{p\in{}ps}m(p)x(p)z(p) - 2d_2d_3\sum_{p\in{}ps}m(p)y(p)z(p) \\ = d_1^2I_{xx} + d_2^2I_{yy} + d_3^2I_{zz} - 2d_1d_2I_{xy} - 2d_1d_3I_{xz} - 2d_2d_3I_{yz}\\ \end{matrix} \end{equation}

    Where we define the products of inertia as follows

    \begin{equation} \begin{matrix} I_{xy} = \sum_{p\in{}ps}m(p)x(p)y(p) & I_{xz} = \sum_{p\in{}ps}m(p)x(p)z(p) & I_{yz} = \sum_{p\in{}ps}m(p)y(p)z(p) \\ \end{matrix} \end{equation}

    This allows for a neat matrix formalism:

    \begin{equation} I_L = \mathbf{D}^T \begin{bmatrix} I_{xx} && -I_{xy} && -I_{xz} \\ -I_{xy} && I_{yy} && -I_{yz} \\ -I_{xz} && -I_{yz} && I_{zz} \\ \end{bmatrix} \mathbf{D} =: \mathbf{D}^T\mathtt{J}\mathbf{D} \end{equation}

    Which defines a matrix \(\mathtt{J}\), the inertia tensor or mass matrix. For the continuous case these values become (over a region R, the book doesn’t put in δ(x,y,z) but do I assume its part of these integrals?):

    \begin{equation} \begin{matrix} I_{xx} = \int_R(y^2 + z^2)dm && I_{yy} = \int_R(x^2 + z^2)dm && I_{zz} = \int_R(x^2 + y^2)dm \\ I_{xy} = \int_R(xy)dm && I_{xz} = \int_R(xz)dm && I_{yz} = \int_R(yz)dm \\ \end{matrix} \end{equation}

    For a single particle of mass m, located at r about \(\mathcal{O}\) rotating about a fixed axis, v = w × r (w is angular velocity), a = \(-r_0\sigma^2\mathbf{R}+\mathbf{\alpha}\times{}\mathbf{r}\) (σ is angular speed and α is angular acceleration). Angular momentum is then:

    \begin{equation} \mathbf{L} = \mathbf{r}\times{}m\mathbf{v} = \mathtt{J}\mathbf{w} \\ ,\mathtt{J} = m(\vert\mathbf{r}\vert^2I - \mathbf{r}\mathbf{r}^T) = m \begin{bmatrix} y^2 + z^2 && -xy && -xz \\ -xy && x^2 + z^2 && -yz \\ -xz && -yz && y^2 + z^2 \\ \end{bmatrix} \end{equation}

    Angular momentum L = \(\mathtt{J}\mathbf{w}\) is analogous to linear momentum p = \(m\mathbf{v}\), torque also follows as:

    \begin{equation} \mathbf{\tau} = \mathbf{r}\times{}m\mathbf{a} = \mathtt{J}\mathbf{\alpha}. \begin{matrix} \end{matrix} \end{equation}

    which is analogous with F=ma. If we choose to factor \(\mathtt{J}\) into eigen-decomposition: \(\mathtt{J}\) = RMRᵀ, M = diag(μ₁, μ₂, μ₃) (eigenvalues of \(\mathtt{J}\)) and R = [U₁,U₂,U₃] (eigenvectors of \(\mathtt{J}\)). Then μᵢ are the principle moments of inertia and Uᵢ are the principle directions of inertia. Torque is the time derivative of angular momentum so lets compute the time derivative of a vector in a moving frame for a rigid body:

    \begin{equation} \begin{matrix} \mathbf{\tau} = \frac{d\mathbf{L}}{dt} = \frac{D\mathbf{L}}{Dt} + \mathbf{w}\times{}\mathbf{L} = \mathtt{J}\frac{d\mathbf{w}}{dt}+\mathbf{w}\times{}(\mathtt{J}\mathbf{w})\\ \end{matrix} \end{equation}

    This may be viewed as the equation of motion for a rigid body. When τ and w are represented in body-coordinate system (coordiante axis directions ar ethe principle directions of inertia):

    \begin{equation} \begin{matrix} \mathbf{\tau} = M\frac{d\textbf{w}}{dt}+\textbf{w}\times{}(M\textbf{w}) \\ \end{matrix} \end{equation}

    Where M is the diagonal matrix as before. This is the Euler equation of motion. (This is achieved by substituting the decompostion in and projecting onto the basis of \(\mathtt{J}\) eigenvectors Rᵀτ and Rᵀw, use the identity (Ra) × (Rb) = R(a×b), the cross product of rotated vectors is the rotation of the cross product of the vectors).

  • Example 2.8
    
    • Compute the inertial tensor for a solid triangle of constant mass density δ ≡ 1, with vertices: (0,0), (a,0), (0,b). So a right triangle of edges: a, b, \(\sqrt{a^2 + b^2}\), we wish to solve for \(\mathtt{J}\). Since the coords are (x,y) (z has no ranging order and thus is 0). The line conecting (a, 0) and (0, b), we want to use parametric form and eliminate t[(a,b)=(a,0),(c,d)=(0,b)]: x = a-ta → 1 - t = x/a → t = -x/a + 1, y = tb → y = b(-x/a + 1). x ∈ [0..a] and y ∈ [0..b(-x/a + 1)]. Mass = Mass Density * Area and for a triangle, area = \(\frac{ab\times{}sin(\gamma)}{2}\) = \(\frac{ab}{2}\) = Mass (as δ = 1).
      

      a, b = var('a b')
      z = 0
      x, y = var('x y')
      I_xx = integral(integral(y^2 + z^2, y, 0, b*(-x/a + 1)), x, 0, a)
      I_yy = integral(integral(x^2 + z^2, y, 0, b*(-x/a + 1)), x, 0, a)
      I_zz = integral(integral(x^2 + y^2, y, 0, b*(-x/a + 1)), x, 0, a)
      I_xy = integral(integral(x*y, y, 0, b*(-x/a + 1)), x, 0, a)
      I_xz = integral(integral(x*z, y, 0, b*(-x/a + 1)), x, 0, a) 
      I_yz = integral(integral(y*z, y, 0, b*(-x/a + 1)), x, 0, a)
      

      This gives (\(I_{xx}\), \(I_{yy}\), \(I_{zz}\), \(I_{xy}\), \(I_{xz}\), \(I_{yz}\)) = (\(\frac{1}{12} \, a b^{3}\), \(\frac{1}{12} \, a^{3} b\), \(\frac{a^{6} b + a^{4} b^{3}}{12 \, a^{3}}\), \(\frac{1}{24} \, a^{2} b^{2}\), 0, 0). Where for completeness we would write in terms with m = \(\frac{ab}{2}\).

  • Example 2.9
    

    Box centered at \(\mathcal{O}\) with dimensions 2a > 0, 2b > 0, and 2c > 0. Vertices are (±a, ±b, ±c). δ ≡ 1.

    1. Compute the inertia tensor for the 8 vertices (each mass is 1 at discrete points):

    m = 8 #; 8 points each mass 1
    a, b, c = var('a b c')
    xs = (-a, a)
    ys = (-b, b)
    zs = (-c, c)
    #; The mass m(p) ≡ 1
    I_xx = sum(1 * (y^2 + z^2) for y in ys for z in zs for x in xs) #; 8*b^2 + 8*c^2 
    I_yy = sum(1 * (x^2 + z^2) for y in ys for z in zs for x in xs) #; 8*a^2 + 8*c^2
    I_zz = sum(1 * (x^2 + y^2) for y in ys for z in zs for x in xs) #; 8*a^2 + 8*b^2
    I_xy = sum(1 * x*y for y in ys for z in zs for x in xs) #; 0
    I_xz = sum(1 * x*z for y in ys for z in zs for x in xs) #; 0
    I_yz = sum(1 * y*z for y in ys for z in zs for x in xs) #; 0
    

    1. Compute the inertia tensor treating the box as a wire frame with wires between vertices:

    We calculate the moment of inertia for each edge then sum them. For this we integrate over 1 dimension and sum the iterative combinations of other dimensions points.

    a, b, c = var('a b c')
    x, y, z = var('x y z')
    def δ(x,y,z):
         return 1
    xs = (-a, a)
    ys = (-b, b)
    zs = (-c, c)
    #; mass is sum of length of edges
    #; Formula is 4*(l×w×h)
    m = 8*(a + b + c)
    #; Sub is with something like expr.subs(z=c)
    eqx = y^2 + z^2
    eqy = x^2 + z^2
    eqz = x^2 + y^2
    I_xx_x = sum(integral(eqx, x, -a, a).subs(y=yb,z=zc) for yb in ys for zc in zs)
    I_xx_y = sum(integral(eqx, y, -b, b).subs(x=xa,z=zc) for xa in xs for zc in zs)
    I_xx_z = sum(integral(eqx, z, -c, c).subs(x=xa,y=yb) for xa in xs for yb in ys)
    I_xx = I_xx_x + I_xx_y + I_xx_z
    #; 8/3*b^3 + 8*b^2*c + 8*b*c^2 + 8/3*c^3 + 8*(b^2 + c^2)*a
    #; we test with
    #; bool(I_xx == (m/3)*(((b+c)^3 + 3*a*(b^2 + c^2))/(a + b + c))) = True against the book
    def getSum(eq):
         return sum(integral(eq, x, -a, a).subs(y=yb,z=zc) for yb in ys for zc in zs) \
                + sum(integral(eq, y, -b, b).subs(x=xa,z=zc) for xa in xs for zc in zs) \
                + sum(integral(eq, z, -c, c).subs(x=xa,y=yb) for xa in xs for yb in ys)
    I_yy = getSum(eqy) #; 8/3*a^3 + 8*a^2*c + 8*a*c^2 + 8/3*c^3 + 8*(a^2 + c^2)*b
    I_zz = getSum(eqz) #; 8/3*a^3 + 8*a^2*b + 8*a*b^2 + 8/3*b^3 + 8*(a^2 + b^2)*c
    eqxy = x*y
    eqxz = x*z
    eqyz = y*z
    I_xy = getSum(eqxy) #; 0
    I_xz = getSum(eqxz) #; 0
    I_yz = getSum(eqyz) #; 0
    

    1. Compute the inertia tensor for the box as a hollow body (surface mass):

    Now two dimensions range at a time!

    a, b, c, x, y, z = var('a b c x y z')
    xs = (-a, a)
    ys = (-b, b)
    zs = (-c, c)
    #; The surfce area formula for a rect prism is: 2(lw + wh + lh)
    l,w,h = 2*a, 2*b, 2*c #; mass density is 1
    m = 2*(l*w + w*h + l*h)
    eqx = y^2 + z^2
    eqy = x^2 + z^2
    eqz = x^2 + y^2
    eqxy = x*y
    eqxz = x*z
    eqyz = y*z
    I_xx = sum(integral(integral(eqx, x, -a, a), y, -b, b).subs(z=zc) for zc in zs) \
           + sum(integral(integral(eqx, x, -a, a), z, -c, c).subs(y=yb) for yb in ys) \
           + sum(integral(integral(eqx, y, -b, b), z, -c, c).subs(x=xa) for xa in xs)
    #; 8/3*b^3*c + 8/3*b*c^3 + 8/3*(b^3 + 3*b*c^2)*a + 8/3*(3*b^2*c + c^3)*a
    #; bool(I_xx == (m/3)*((a*(b+c)^3 + b*c*(b^2 + c^2))/(a*b + a*c + b*c))) => True!
    def getSum(eq):
        return sum(integral(integral(eq, x, -a, a), y, -b, b).subs(z=zc) for zc in zs) \
           + sum(integral(integral(eq, x, -a, a), z, -c, c).subs(y=yb) for yb in ys) \
           + sum(integral(integral(eq, y, -b, b), z, -c, c).subs(x=xa) for xa in xs)
    I_yy = getSum(eqy) #; 8/3*a^3*c + 8/3*a*c^3 + 8/3*(a^3 + 3*a*c^2)*b + 8/3*(3*a^2*c + c^3)*b
    I_zz = getSum(eqz) #; 8/3*a^3*b + 8/3*a*b^3 + 8/3*(a^3 + 3*a*b^2)*c + 8/3*(3*a^2*b + b^3)*c
    I_xy = getSum(eqxy) #; 0
    I_xz = getSum(eqxz) #; 0 
    I_yz = getSum(eqyz) #; 0
    

    1. Compute the inertia tensor for the box treated as a solid (volume mass):

    Now the entire region is being integrated over!

    a, b, c, x, y, z = var('a b c x y z')
    #; mass is just the volume as mass density is 1
    l,w,h = 2*a, 2*b, 2*c
    m = l * w * h
    eqx = y^2 + z^2
    eqy = x^2 + z^2
    eqz = x^2 + y^2
    eqxy = x*y
    eqxz = x*z
    eqyz = y*z
    def getInt(eq):
        return integral(integral(integral(eq, x, -a, a), y, -b, b), z, -c, c)
    
    I_xx = getInt(eqx) #; 8/3*(b^3*c + b*c^3)*a
    #; bool(I_xx == (m*(b^2 + c^2))/(3)) => True
    I_yy = getInt(eqy) #; 8/3*(a^3*c + a*c^3)*b
    I_zz = getInt(eqz) #; 8/3*(a^3*b + a*b^3)*c
    I_xy = getInt(eqxy) #; 0
    I_xz = getInt(eqxz) #; 0
    I_yz = getInt(eqyz) #; 0
    

  • TODO Exercise 2.12   TypicallyHidden
    
  • TODO Exercise 2.13   TypicallyHidden
    
  • TODO Exercise 2.14   TypicallyHidden
    
  • TODO Exercise 2.15   TypicallyHidden
    
• Mass an Inertia Tensor of a solid Polyhedron
  

  This problem becomes more difficult for non-constant mass density, but here we assume δ ≡ 1. We will be using the Mirtich algorithm for constant mass polyhedra. It uses the divergence theorem for reducing volume to surface integrals (Green’s theorem, the 2D-divergence analog is used reduce the planar integrals further to line integrals). The full picture is: Divergence Theorem => Projection => Green’s Theorem, then propagate the values backwards.

  • Reduction of Volume Integrals
    

    The mass, center of mass, and inertial tensor require computing volume integrals of the type:

    \begin{equation} \begin{matrix} \int_Vp(x,y,z)dV \end{matrix} \end{equation}

    V is the volumetric region and dV is the infinitesimal measure of the region. p(x,y,z) is the polynomial function selected from {1, x, y, z, x², y², z², xy, xz, yz}. Here V is the region bounded by a simple polyhedron. Let’s apply the divergence theorem:

    \begin{equation} \begin{matrix} \int_{V}p(x,y,z)dV = \int_{V}\nabla{}\cdot{}FdV = \int_{S}N\cdot{}FdS \end{matrix} \end{equation}

    S is the boundary of the polyhedron, a union of triangular face and dS is a infinitesimal measure of surface area. F(x,y,z) is chosen as ∇⋅F = p, N is the outward pointing unit surface normals. Here are the choosen F’s from the paper:

    p	F	p	F
    1	(x,0,0)	y²	(0, y³/3, 0)
    x	(x²/2, 0, 0)	z²	(0, 0, z³/3)
    y	(0, y²/2, 0)	xy	(x²y/2, 0, 0)
    z	(0, 0, z²/2)	xz	(0, 0, z²x/2)
    x²	(x³/3, 0, 0)	yz	(0, y²z/2, 0)

    Computational effort is now solcing ∫ₛN⋅FdS. The boundary S is just the union of polyhedral face \(\mathcal{F}\). The unit length Normal for a given face \(\mathcal{F}\) is \(N_{\mathcal{F}}=(\hat{\eta}_x, \hat{\eta}_y, \hat{\eta}_y)\). This surface intergral is then decomposed into:

    \begin{equation} \begin{matrix} \int_S\mathbf{N}\cdot{}\mathbf{F}ds = \sum_{\mathcal{F}\in{}S}\int_{\mathcal{F}}\mathbf{N_\mathcal{F}}\cdot{}\mathbf{F}dS \\ \end{matrix} \end{equation}

    The integrals we care about are then:

    \begin{equation} \begin{matrix} \int_{V}dV = \sum_{\mathcal{F}\in{}S}\hat{\eta}_x\int_{\mathcal{F}}xdS && \int_{V}y^2dV = \frac{1}{3}\sum_{\mathcal{F}\in{}S}\hat{\eta}_y\int_{\mathcal{F}}y^3dS \\ \int_{V}xdV = \frac{1}{2}\sum_{\mathcal{F}\in{}S}\hat{\eta}_x\int_{\mathcal{F}}x^2dS && \int_{V}z^2dV = \frac{1}{3}\sum_{\mathcal{F}\in{}S}\hat{\eta}_z\int_{\mathcal{F}}z^3dS \\ \int_{V}ydV = \frac{1}{2}\sum_{\mathcal{F}\in{}S}\hat{\eta}_y\int_{\mathcal{F}}y^2dS && \int_{V}xydV = \frac{1}{2}\sum_{\mathcal{F}\in{}S}\hat{\eta}_x\int_{\mathcal{F}}x^2ydS \\ \int_{V}zdV = \frac{1}{2}\sum_{\mathcal{F}\in{}S}\hat{\eta}_z\int_{\mathcal{F}}z^2dS && \int_{V}y^2zdV = \frac{1}{2}\sum_{\mathcal{F}\in{}S}\hat{\eta}_y\int_{\mathcal{F}}y^2zdS \\ \int_{V}x^2dV = \frac{1}{3}\sum_{\mathcal{F}\in{}S}\hat{\eta}_x\int_{\mathcal{F}}x^3dS && \int_{V}xzdV = \frac{1}{2}\sum_{\mathcal{F}\in{}S}\hat{\eta}_z\int_{\mathcal{F}}z^2xdS \\ \end{matrix} \end{equation}

    or more conveniently (q is a \(\int_\mathcal{F}\) body from above, \(\mathcal{l}\) := x | y | z):

    \begin{equation} \begin{matrix} \hat{\eta}_\mathcal{l}\int_{\mathcal{F}}q(x,y,z)dS \end{matrix} \end{equation}
  • Computation by Reduction to Line integrals.
    

    When q = 1, we don’t need to compute the surface integral: we use this case as motivation. \(\int_{\mathcal{F}}dS\) is the area A of the polygonal face. The area may be computed as a 3D quantity:

    \begin{equation} \begin{matrix} A = \int_{\mathcal{F}}ds = \frac{1}{2}\mathbf{N}_\mathcal{F}\cdot{}\sum_{i=0}^{n-1}\mathbf{P}_i\times{}\mathbf{P}_{i+1} \\ \end{matrix} \end{equation}

    Where \(\mathbf{P}_i\) = (xᵢ, yᵢ, zᵢ) are the polygon vertices, 0 ≤ i ≤ n - 1, where the vertices are ordered counterclockwise relative to the face normal \(\mathbf{N}_\mathcal{F}\), additionally indexes are modular i.e. Pᵢ ≡ Pₙ. This formula is from [Arv91] and is derived using Stoke’s theorem, but not optimal. Instead our goal is to project the polygon onto a coordinate plane and compute the area in 2D (adjusting to account for projection). The plane of the polynomial face is \(\hat{\eta}_xx\) + \(\hat{\eta}_yy\) + \(\hat{\eta}_zz\) + w = 0, w = \(-\mathbf{N}_\mathcal{F}\cdot{}\mathbf{P}_0\). If \(\hat{\eta}_z\) ≠ 0 then we can project into xy-plane, z = f(x,y) = \(\frac{-(w + \hat{\eta}_xx + \hat{\eta}_yy)}{\hat{\eta_z}}\). The infinitesimal surface area in terms of independent x and y becomes:

    \begin{equation} \begin{matrix} dS = \sqrt{1 + (\frac{\partial{}f}{\partial{}x})^2 + (\frac{\partial{}f}{\partial{}y})^2}\ dxdy \\ dS = \sqrt{1 + (\frac{-\hat{\eta}_x}{\hat{\eta}_z})^2 + (\frac{-\hat{\eta}_y}{\hat{\eta}_z})^2}\ dxdy \\ dS = \frac{1}{\vert{}\hat{\eta}_z\vert{}}dxdy \\ \vert{}\mathbf{N}_\mathcal{F}\vert{} = 1 \\ \end{matrix} \end{equation}

    The surface integral now is:

    \begin{equation} \begin{matrix} \int_\mathcal{F}dS = \frac{1}{\vert{}\hat{\eta}_z\vert{}}\int_Rdxdy\\ \end{matrix} \end{equation}

    R is the bounding region by the projected polygon. The vertices of the polygon are \(\mathbf{Q}_i\) = (xᵢ, yᵢ) using Green’s theorem:

    \begin{equation} \begin{matrix} \int_Rp(x,y)dxdy = \int_R\nabla\cdot\mathbf{G}dxdy = \int_L\mathbf{M}\cdot\mathbf{G}ds \end{matrix} \end{equation}

    L is the boundary of the region R and ds is the infinitesimal measure of arc length. G(x,y) is chosen such that ∇ ⋅ G = p. The vector M is the ourward-unit-length curve normals. In the special case p(x,y) = 1 and L is a polygon, as such many such G can work but we choose G(x,y) = (x,0). The boundary is the union of edges, \(\mathcal{E}\), and outward unit normal being \(\mathbf{M}_\mathcal{E}\), decomposing the area integral to (the second equation to see similarity above):

    \begin{equation} \begin{matrix} \int_L\mathbf{M}\cdot\mathbf{G}ds=\sum_\mathcal{E}\int_\mathcal{E}\mathbf{M}_\mathcal{E}\cdot\mathbf{G}ds\\ \int_S\mathbf{N}\cdot{}\mathbf{F}ds = \sum_{\mathcal{F}\in{}S}\int_{\mathcal{F}}\mathbf{N_\mathcal{F}}\cdot{}\mathbf{F}dS \\ \end{matrix} \end{equation}

    If edge \(\mathcal{E}\) is connecting Qᵢ to Qᵢ₊₁ and has length Lᵢ = |Qᵢ₊₁ - Qᵢ|, then the edge is parameterized by (x(s), y(s)) = (1 - \(\frac{s}{L_i}\))Qᵢ + (\(\frac{s}{L_i}\))Qᵢ₊₁ for s ∈ [0..Lᵢ]. The edge normals are: Mᵢ = sign(\(\hat{\eta}_z)\frac{(y_{i+1} - y_{i}, x_i - x_{i+1})}{L_i}\).

    \begin{equation} \begin{matrix} \int_L\mathbf{M}\cdot\mathbf{G}ds=\sum_{i=0}^{n-1}\int_0^{L_i}M_i\cdot(x(s),0)ds = \frac{sign(\hat{\eta}_z)}{2}\sum_{i=0}^{n-1}x_i(y_{i+1}-y_{i-1})\\ \end{matrix} \end{equation}

    There is a possible ill-conditioning here when \(\hat{\eta}_z\) ≈ 0, so we can choose to project on another coordinate plane, but these can also be ill conditioned (but not both as one component of \(\mathbf{N}_\mathcal{F}\) ≥ \(\frac{1}{\sqrt{3}}\), as its normal, thus its project will be conditioned well). Polygonal faces are computed (R’s are the respective regions of projection of the polygonal face):

    \begin{equation} A = \begin{Bmatrix} \frac{\int_{R_{xy}}dxdy}{\vert{}\hat{\eta}_z\vert{}} = \frac{1}{2\hat{\eta}_z}\sum_{i=0}^{n-1}x_i(y_{i+1}-y_{i-1}), && \vert{}\hat{\eta}_z\vert{} = max(\vert{}\hat{\eta}_x\vert{},\vert{}\hat{\eta}_y\vert{},\vert{}\hat{\eta}_z\vert{}) \\ \frac{\int_{R_{yz}}dydz}{\vert{}\hat{\eta}_x\vert{}} = \frac{1}{2\hat{\eta}_x}\sum_{i=0}^{n-1}y_i(z_{i+1}-z_{i-1}), && \vert{}\hat{\eta}_x\vert{} = max(\vert{}\hat{\eta}_x\vert{},\vert{}\hat{\eta}_y\vert{},\vert{}\hat{\eta}_z\vert{}) \\ \frac{\int_{R_{zx}}dzdx}{\vert{}\hat{\eta}_y\vert{}} = \frac{1}{2\hat{\eta}_y}\sum_{i=0}^{n-1}z_i(x_{i+1}-x_{i-1}), && \vert{}\hat{\eta}_y\vert{} = max(\vert{}\hat{\eta}_x\vert{},\vert{}\hat{\eta}_y\vert{},\vert{}\hat{\eta}_z\vert{}) \\ \end{Bmatrix} \end{equation}

    The construction of ∫ₛq(x,y,z)dS is handled in a similar manner (now with integrand variables) we replace the dependant variable using the appropriate formulation of the plane equation:

    \begin{equation} \int_Sq(x,y,z)dS = \begin{Bmatrix} \frac{1}{\vert{}\hat{\eta}_z\vert{}}\int_{R_{xy}}q(x,y,-(\frac{w + \hat{\eta}_xx + \hat{\eta}_yy}{\hat{\eta}_z})), && \vert{}\hat{\eta}_z\vert{} = max(\vert{}\hat{\eta}_x\vert{},\vert{}\hat{\eta}_y\vert{},\vert{}\hat{\eta}_z\vert{}) \\ \frac{1}{\vert{}\hat{\eta}_x\vert{}}\int_{R_{yz}}q(-(\frac{w + \hat{\eta}_yy + \hat{\eta}_zz}{\hat{\eta}_x}),y,z), && \vert{}\hat{\eta}_x\vert{} = max(\vert{}\hat{\eta}_x\vert{},\vert{}\hat{\eta}_y\vert{},\vert{}\hat{\eta}_z\vert{}) \\ \frac{1}{\vert{}\hat{\eta}_y\vert{}}\int_{R_{zx}}q(x,-(\frac{w + \hat{\eta}_zz + \hat{\eta}_xx}{\hat{\eta}_y}),z), && \vert{}\hat{\eta}_y\vert{} = max(\vert{}\hat{\eta}_x\vert{},\vert{}\hat{\eta}_y\vert{},\vert{}\hat{\eta}_z\vert{}) \\ \end{Bmatrix} \end{equation}

    Each integrand is a polynomial of two variables, G must be chosen such that ∇ ⋅ G is that polynomial. When α is the max term of (x,y,z) for η the polynomial will be in terms of ¬α, thus:

    \begin{equation} \int_Sq(x,y,z)dS = \begin{matrix} \frac{1}{\hat{\eta}_z}\sum_{i=0}^{n-1}(y_{i+1}-y_i,x_i-x_{i+1})\cdot{}\int_0^1G((1-t)Q_i + tQ_{i+1})dt \\ \end{matrix} \end{equation}

    To see the generalized (α, β, γ) permutations [(x,y,z),(y,z,x),(z,x,y)] for these equations see page: 70.

  • Computation by Direct parameterization of Triangles
    

    The formule greatly simplify when using triangle faces for the polyhedron. Let the triangular face be counterclockwise oriented and have vertices \(\mathbf{P}_i\) = (xᵢ, yᵢ, zᵢ), 0 ≤ i ≤ 2. Two edges are (for 1 ≤ i ≤ 2):

    \begin{equation} \begin{matrix} \mathbf{E_i} = \mathbf{P_i} - \mathbf{P_0} = (x_i - x_0, y_i - y_0, z_i - z_0) = (\alpha_i, \beta_i, \gamma_i) \\ \end{matrix} \end{equation}

    Then the parameterized face is (u ≥ 0, v ≥ 0, u + v ≤ 1)

    \begin{equation} \begin{matrix} \mathbf{P}(u, v) = \mathbf{P}_0 + u\mathbf{E}_1 + v\mathbf{E}_2 \\ = (x_0 + \alpha_1u + \alpha_2v, y_0 + \beta_1u + \beta_2v, z_0 + \gamma_1u + \gamma_2v) \\ = (x(u,v), y(u,v), z(u,v))\\ \end{matrix} \end{equation}

    Infinitesimal measure of surface area is:

    \begin{equation} dS = \begin{vmatrix} \frac{\partial{}\mathbf{P}}{\partial{}u} \times{} \frac{\partial{}\mathbf{P}}{\partial{}v} \end{vmatrix}dudv = \begin{vmatrix} \mathbf{E}_1 \times \mathbf{E}_2 \end{vmatrix}dudv \end{equation}

    The outer-pointing unit-length face normal is:

    \begin{equation} \mathbf{N}_\mathcal{F} = \begin{matrix} \frac{\mathbf{E}_1\times{}\mathbf{E}_2}{\vert{}\mathbf{E}_1\times{}\mathbf{E}_2\vert{}} = \frac{(\beta_1\gamma_2 - \beta_2\gamma_1, \alpha_2\gamma_1 - \alpha_1\gamma_2, \alpha_1\beta_2 - \alpha_2\beta_1)}{\vert{}\mathbf{E}_1\times{}\mathbf{E}_2\vert{}} = \frac{(\delta_0, \delta_1, \delta_2)}{\vert{}\mathbf{E}_1\times{}\mathbf{E}_2\vert{}} \end{matrix} \end{equation}

    From here we implement a sample program to calculate these in full (Start from page 73, top). Here I rework some magic source code into a fortran version:

    module Polyhedral
      implicit none
      type Point
          real(kind=8) :: x,y,z
      contains
          ! we name the pass self here (nopass is also an attr)
          procedure, private, pass(self) :: point_add_point
          procedure, private, pass(self) :: point_sub_point
          procedure, private, pass(self) :: point_cross_point
          procedure, private, pass(self) :: point_dot_point
          generic, public :: operator(+) => point_add_point
          generic, public :: operator(-) => point_sub_point
          generic, public :: operator(.cross.) => point_cross_point
          generic, public :: operator(.dot.) => point_dot_point
      end type Point
    
    
      ! You would want a function -> convert to matrix for this type as well
      type InertiaTensor
          real(kind=8) :: xx, yy, zz, xy, xz, yz
      end type InertiaTensor
    
      integer, parameter :: long = selected_real_kind(9, 99) ! 8 on my machine
      ! So I use 8 directly, this is bad practice make a module of kinds or something
      real(kind=8), parameter :: OneDiv6 = 1._long/6._long
      real(kind=8), parameter :: OneDiv24 = 1._long/24._long
      real(kind=8), parameter :: OneDiv60 = 1._long/60._long
      real(kind=8), parameter :: OneDiv120 = 1._long/120._long
    contains
      pure function point_add_point(self, o) result(c)
          class(Point), intent(in) :: self
          type(Point), intent(in) :: o
          type(Point) :: c
          c%x = self%x + o%x
          c%y = self%y + o%y
          c%z = self%z + o%z
      end
      pure function point_sub_point(self, o) result(c)
          class(Point), intent(in) :: self
          type(Point), intent(in) :: o
          type(Point) :: c
          c%x = self%x - o%x
          c%y = self%y - o%y
          c%z = self%z - o%z
      end function point_sub_point
    
      ! RHR cross product
      pure function point_cross_point(self, o) result(c)
          class(point), intent(in) :: self
          type(point), intent(in) :: o
          type(point) :: c
          c%x = self%y * o%z - self%z * o%y
          c%y = self%z * o%x - self%x * o%z
          c%z = self%x * o%y - self%y * o%x
      end function point_cross_point
    
      pure function point_dot_point(self, o) result(c)
          class(point), intent(in) :: self
          type(point), intent(in) :: o
          real(kind=8) :: c
          c = self%x * o%x + self%y * o%y + self%z * o%z
      end function point_dot_point
    
    
    
      subroutine Subexpressions(w0, w1, w2, f1, f2, f3, g0, g1, g2)
          ! Fortran is a pass by reference, opposed to C's pass by value
          real(kind=8), intent(inout) :: w0, w1, w2, f1, f2, f3, g0, g1, g2
          real(kind=8) :: temp0, temp1, temp2
          temp0 = w0 + w1
          f1 = temp0 + w2
          temp1 = w0 * w0
          temp2 = temp1 + w1 * temp0
          f2 = temp2 + w2 * f1
          f3 = w0 * temp1 + w1 * temp2 + w2 * f2
          g0 = f2 + w0 * (f1 + w0)
          g1 = f2 + w1 * (f1 + w1)
          g2 = f2 + w2 * (f1 + w2)
      end subroutine Subexpressions
    
    
      subroutine Compute(p, tmax, index, mass, cm, inertia)
          type(Point), intent(in), dimension(:) :: p
          integer, intent(in) :: tmax
          integer, intent(in) :: index(:)
          real(kind=8), intent(inout) :: mass
          type(Point), intent(inout) :: cm
          type(InertiaTensor), intent(inout) :: inertia
          ! logical :: bodyCoords ! I forgo this but its simple
          ! Array temporaries, we have bound freedom in fortran
          ! Could do 0:9
          real(kind=8) :: integral(10);
          ! The temporaries
          integer :: t;
          integer :: i0, i1, i2
          type(Point) :: v0, v1, v2, v1mv0, v2mv0, N
          real(kind=8) :: f1x, f2x, f3x, &
            & g0x, g1x, g2x, f1y, f2y, f3y, g0y, g1y, g2y, &
            & f1z, f2z, f3z, g0z, g1z, g2z
          ! Order: 1, x, y, z, x^2, y^2, z^2, xy, yz, zx
          data integral(1:10)/10*0./ ! Initialize as 10 0.'s
          ! Loops are inclusive in Fortran
          do t = 0, (tmax - 1)
             ! Get the vertices of the triangle
             i0 = index(3 * t + 1)
             i1 = index(3 * t + 2)
             i2 = index(3 * t + 3)
             v0 = p(i0)
             v1 = p(i1)
             v2 = p(i2)
             ! Get the cross product of edges
             ! And compute the Normals
             v1mv0 = v1 - v0
             v2mv0 = v2 - v0
             N = v1mv0 .cross. v2mv0
             ! Compute integral terms
             call Subexpressions(v0%x, v1%x, v2%x, &
             & f1x, f2x, f3x, g0x, g1x, g2x)
             call Subexpressions(v0%y, v1%y, v2%y, &
             & f1y, f2y, f3y, g0y, g1y, g2y)
             call Subexpressions(v0%z, v1%z, v2%z, &
             & f1z, f2z, f3z, g0z, g1z, g2z)
             ! Update the integral
             integral(1) = integral(1) + N%x*f1x
             integral(2) = integral(2) + N%x*f2x
             integral(3) = integral(3) + N%y*f2y
             integral(4) = integral(4) + N%z*f2z
             integral(5) = integral(5) + N%x*f3x
             integral(6) = integral(6) + N%y*f3y
             integral(7) = integral(7) + N%z*f3z
             integral(8) = integral(8) + N%x*(v0%y*g0x + v1%y*g1x + v2%y*g2x)
             integral(9) = integral(9) + N%y*(v0%z*g0y + v1%z*g1y + v2%z*g2y)
             integral(10) = integral(10) + N%z*(v0%x*g0z + v1%x*g1z + v2%x*g2z)
          end do
          ! Update by weight
          integral(1) = integral(1) * OneDiv6
          integral(2) = integral(2) * OneDiv24
          integral(3) = integral(3) * OneDiv24
          integral(4) = integral(4) * OneDiv24
          integral(5) = integral(5) * OneDiv60
          integral(6) = integral(6) * OneDiv60
          integral(7) = integral(7) * OneDiv60
          integral(8) = integral(8) * OneDiv120
          integral(9) = integral(9) * OneDiv120
          integral(10) = integral(10) * OneDiv120
          ! get Mass
          mass = integral(1)
          cm%x = integral(2) / mass
          cm%y = integral(3) / mass
          cm%z = integral(4) / mass
          ! Set the inertia tensor
          ! Relative to world origin
          inertia%xx = integral(6) + integral(7)
          inertia%xy = -integral(8)
          inertia%xz = -integral(10)
          inertia%yy = integral(5) + integral(7)
          inertia%yz = -integral(9)
          inertia%zz = integral(5) + integral(6)
          ! Relative to center of mass
          ! This only should happen if we had a com view
          inertia%xx = inertia%xx - mass * (cm%y * cm%y + cm%z * cm%z)
          inertia%yy = inertia%yy - mass * (cm%z * cm%z + cm%x * cm%x)
          inertia%zz = inertia%zz - mass * (cm%x * cm%x + cm%y * cm%y)
          inertia%xy = inertia%xy + mass * cm%x * cm%y
          inertia%yz = inertia%yz + mass * cm%y * cm%z
          inertia%xz = inertia%xz + mass * cm%z * cm%x
      end subroutine Compute
    
    end module Polyhedral
    
    
    program Run
      use Polyhedral
      implicit none
      ! Example of using
    end program Run
    

Energy

• Work and Kinetic Energy
  

  Given p :: Particle ~ {mass = m, 1D = {D :: UnitVector}}, F :: Force ~ {Constant}, θ :: Angle ~ {Between(F,D)}. If F is appiled to p while p is moving along D over a distance L = |x₁ - x₀| in a time interval. The work done by F on p is

  \begin{equation} \begin{matrix} W = (\vert{}F\vert{}cos(\theta))L = (F\cdot{}D)L \end{matrix} \end{equation}

  Force in the direction of the D is (F⋅D)D, and work is this component times length traveled. If the path the particle takes is any x(t) :: Curve ~ {Smooth}, then we use infinitesimals and calculus. The direction updates to be the arclength instance direction, D = dx/ds, making infinitesimal work:

  \begin{equation} \begin{matrix} dW = (F\cdot{}D)ds = (F\cdot\frac{dx}{ds})ds \end{matrix} \end{equation}

  and total work:

  \begin{equation} \begin{matrix} W = \int_0^LF\cdot{}\frac{dx}{ds}ds \end{matrix} \end{equation}

  We typically want work as a function of time (over interval [t₀..t₁], and v(t) :: Velocity) instead of arclength so we update the equations to:

  \begin{equation} \begin{matrix} dW = (F\cdot{}dx)dt \\ W = \int_{t_0}^{t_1}F\cdot\frac{dx}{dt}dt = \int_{t_0}^{t_1}F\cdot{}vdt \end{matrix} \end{equation}

  It is also to remember that this is the work done by a singular force and the particle can be acted on by other forces as well and thus have postion not dependant on a single F. If F is a net force on a particle which controls its position and if mass is consant then we can use Newtons Second law, F = ma:

  \begin{equation} \begin{matrix} W = \frac{m}{2}(\vert{}v(t_1)\vert{}^2 - \vert{}v(t_0)\vert{}^2) \end{matrix} \end{equation}

  The quantity W is the energy required to change the particle’s velocity from v(t₀) to v(t₁). If v(t₀) ≡ 0 then W is called kinetic energy, which is the particle from rest case. The general kinetic energy for moving particle is:

  \begin{equation} \begin{matrix} T = \frac{m}{2}\vert{}\mathbf{v}\vert{}^2 \end{matrix} \end{equation}

  Kinetic energy is specified with respsect to the inertial frame of reference, the same sued to specify the postion x(t) of the particle. Kinetic energy is additive in a system of particles (summation for discrete and integration for continous materials).

• Conservative forces and potential energy
  

  Finding the force on a particle traveling along a path x(t), t ∈ [t₀, t₁] is defined previously. The force in general is dependant on the path taken, x(t), but when it does not matter the path taken the force is said to be conservative.

  • Example (gravity):
    

    Two endpoints for a particle are (x₀, y₀, z₀) and (x₁, y₁, z₁) along path (x(t), y(t), z(t)), t ∈ [t₀, t₁]. The particle is subject to gravity: F = -mgk (i, j, k are cardinal unit vectors), the work of gravity is:

    # using Pkg;
    # Pkg.add("Symbolics");
    # Pkg.add("SymbolicNumericIntegration");
    using Symbolics;
    using SymbolicNumericIntegration;
    using Latexify;
    # W = integrate(F⋅v, t, t₀, t₁)
    # k ⋅ (ẋ(t), ẏ(t), ż(t)) ≡ ż(t)
    # W = integrate(-mgk⋅(ẋ(t), ẏ(t), ż(t)), t, t₀, t₁)
    # W = -m*g*[z(t), t₀, t₁]
    @variables t₀ t₁ m g k t ż(..) z(..);
    Work = -m * g * (z(t₁) - z(t₀));
    print(latexify(Work))
    

    \begin{equation} - g m \left( z\left( \mathtt{t{_1}} \right) - z\left( \mathtt{t{_0}} \right) \right) \end{equation}

    Here the work is independant of the path it is always a C * ’Difference in height of endpoints’. Hence, gravity is a conservative force.

  • Example (spring):
    

    One end of a spsring is attached to (0,0,0) and the other is attached to a particle varing on time x(t), t ∈ [t₀, t₁]. The spring constant is c > 0 and unstreached length is L. The force by the spring is then F = \(-c(x-\frac{Lx}{\vert{}x\vert{}})\) in the direction restoring to L. We can make l = \(\frac{Lx}{\vert{}x\vert{}}\) as a constant vector length L in direction of x. We use the line integral style work integrals.

    # x, ẋ, t, t₀, t₁, c, l, l̇ = var('x ẋ t t₀ t₁ c l l̇')
    L, c, t0, t1, t, dly, dlz = var('L c t0 t1 t dly dlz')
    o = vector([0, 0, 0])
    x = function('x', nargs=1)(t) # Atleast one more direction is supported here
    x = vector([x, 0, 0]) # y and z not changing
    l = L * x / abs(x)
    #xv = vector([x(t)])
    dx = derivative(x, t)
    l2 = l.dot_product(l)
    l2 # autoprint
    # dl ⟂ l because l⋅l = L² ∴ ||l|| = √(l⋅l) = √(L²) = L
    # Const Mag vector is ⟂ its derivative is used here
    # We know it is in y and z directions
    dl = vector([0, dly, dlz])
    l.dot_product(dl) # = 0
    # d/dt[a(t) ⋅ b(t)] = ȧ(t)⋅b(t) + a(t)⋅ḃ(t) = if a = b = y then 2[y(t)⋅ẏ(t)]
    y = x - l
    dy = dx - dl
    y.dot_product(dy) == 1/2 * derivative(y.dot_product(y), t)
    W = -c * integral(derivative(1/2 * abs(x - l) ** 2, t), t, t0, t1)
    # W == (-c/2)*(abs(x(t1) - l(t1))**2 - abs(x(t0) - l(t0))**2)
    W
    

    L^2
    0
    -(L*x(t)/sqrt(x(t)^2) - x(t))*diff(x(t), t) == -(L*x(t)^2*diff(x(t), t)/(x(t)^2)^(3/2) - L*diff(x(t), t)/sqrt(x(t)^2) + diff(x(t), t))*(L*x(t)/sqrt(x(t)^2) - x(t))
    1/2*(x(t0)^2 - x(t1)^2 - 2*L*sqrt(x(t0)^2) + 2*L*sqrt(x(t1)^2))*c
    

    Which means that work only depends on the endpoints and noth the path connecting them, so springs are a conservative force.

  • Example (Non-conservative)
    

    A particle moves from (0,0,0) to (1,1,0) and is subject to F = (y, -x, 0). This is a non-conservative force. The work done by the force on the particle traveling across the line segment (t, t, 0) [∴ x = t and y = t], t ∈ [0, 1]. Secondly we consider a arc (t,t²,0)

    x, y, z = var('x y z')
    t0, t1 = (0, 1)
    s = vector([0, 0, 0])
    e = vector([1, 1, 0])
    line = vector([t, t, 0])
    F = vector([y, -x, 0]).subs(x=t, y=t)
    # W = integral(F ⋅ v, t, 0, 1)
    # W1 = integral(F.dot_product(derivative(pos, t)), t0, t1)
    #W1
    W1 = integral(F.dot_product(derivative(line, t)), t, t0, t1)
    # Whereas if we had it taking a arc
    arc = vector([t, t**2, 0])
    # Define the force substitution properly
    F = vector([y, -x, 0]).subs(x=t, y=t**2)
    W2 = integral(F.dot_product(derivative(arc, t)), t, t0, t1)
    W1 != W2
    

    0 != (-1/3)
    

    As we can see the path taken by the particle matters so we can say the force is non-conservative. A major example of a non-conservative force is viscosity.

  • Continuing
    

    Consider particles, ps :: Particle(p) ~ {mass = m, pos = (x,y,z)}. These positions are the reference system conservative force Fᵢ are applied to the particle in a general system, where the pos and work in the general system is transferred to the reference system. The work done by transferring the forces is referred to as the potential energy that the general system has with repset to the reference system. This is denoted with V:

    \begin{equation} \begin{matrix} V = -\sum_{i=1}^{p}\int_{0}^{1}F_{i}\cdot{}v_{i}dt \end{matrix} \end{equation}

    The paths implictly connect (xᵢ, yᵢ, zᵢ) at t = 0 to (xᵢ, yᵢ, zᵢ) at t = 1 as the forces are conservative so any smooth paths parameterized by t ∈ [0,1] will work. V is only dependant on the general system and the referneces are assumed as constants (i.e. a function of 3p variables).

    \begin{equation} \begin{matrix} V = V(x_1, y_1, y_1, ..., x_p, y_p, z_p) \\ \end{matrix} \end{equation}

    This leads to the derivative and force formulations of (use an exactness test, pg 83 for conservative force checking: ∇×F = 0)

    \begin{equation} \begin{matrix} dV = \sum_{i=1}^{p}(\frac{\partial{}V}{\partial{}x_i}dx_i + \frac{\partial{}V}{\partial{}y_i}dy_i + \frac{\partial{}V}{\partial{}z_i}dz_i) \\ dV = -\sum_{i=1}^{p}F_{i}\cdot{}d\mathbf{x}_{i} = -\sum_{i=1}^{p}(F_{x_i}dx_{i} + F_{y_i}dy_{i} + F_{z_i}dz_{i}) \\ (F_{x_i}, F_{y_i}, F_{z_i}) = (-\frac{\partial{}V}{\partial{}x_i},-\frac{\partial{}V}{\partial{}y_i},-\frac{\partial{}V}{\partial{}z_i}) \\ F_i = -\nabla_iV \\ \end{matrix} \end{equation}
  • Example (gravity conservation, spring conservation, and non-conservation)
    

    (F₁, F₂, F₃) = (0, 0, -mg)

    # Gravity
    m, g, x, y, z = var('m g x y z')
    F = vector([0, 0, -m*g])
    # Exactness conditions
    derivative(F[0], y) == derivative(F[0], x)
    derivative(F[0], z) == derivative(F[2], x)
    derivative(F[1], z) == derivative(F[2], y)
    # So the partial derivatives for each is 0
    # Spring
    c = var('c')
    F = -c * vector([x, y, z])
    derivative(F[0], y) == derivative(F[1], x)
    derivative(F[0], z) == derivative(F[2], x)
    derivative(F[1], z) == derivative(F[2], y)
    # Non-conservative
    F = vector([y, -x, 0])
    derivative(F[0], y) == derivative(F[1], x)
    derivative(F[0], z) == derivative(F[2], x)
    derivative(F[1], z) == derivative(F[2], y)
    

    0 == 0
    0 == 0
    0 == 0
    0 == 0
    0 == 0
    0 == 0
    1 == -1
    0 == 0
    0 == 0
    

    Having a conservative force means the total energy is conserved E = T + V, \frac{dE}{dt} = 0. Page 84.