feat(ClassicalMechanics): inertia tensor about a point and the parallel-axis theorem

Physlib/ClassicalMechanics/RigidBody/API-map.yaml

@@ -6,7 +6,9 @@ Overview: |
     The key data structure is `RigidBody d`, a rigid body in `d` dimensions specified by
     its mass distribution (a linear map from smooth functions on `Space d` to `ℝ`). The API
     provides the body's total mass, centre of mass and inertia tensor, and works out the
-    solid sphere as a concrete example.
+    solid sphere as a concrete example. The inertia tensor is also available about an arbitrary
+    point, and the parallel-axis theorem expresses it in terms of the inertia tensor about the
+    centre of mass.
 
     A rigid body *in motion*, `RigidBodyMotion d`, additionally records the centre-of-mass
     trajectory and the time-dependent orientation. From these the API builds the
@@ -49,6 +51,14 @@ Requirements:
     done: true
     location: Physlib/ClassicalMechanics/RigidBody/Basic.lean (RigidBody.inertiaTensor_symmetric)
 
+  - description: >
+      The API allows the calculation of the inertia tensor about an arbitrary point, and contains
+      the parallel-axis theorem relating it to the inertia tensor about the centre of mass.
+    done: true
+    location: >
+      Physlib/ClassicalMechanics/RigidBody/Basic.lean
+      (RigidBody.inertiaTensorAbout, RigidBody.inertiaTensorAbout_eq_centerOfMass_add_pointMass)
+
   - description: >
       The API provides the solid sphere as a worked example, with its total mass and centre
       of mass computed.

Physlib/ClassicalMechanics/RigidBody/Basic.lean

@@ -19,14 +19,19 @@ In this module we will define the basic properties of a rigid body, including
 - mass
 - center of mass
 - inertia tensor
+- inertia tensor about an arbitrary point
+
+The moments of the last are taken about the given point rather than about the origin of the
+reference frame. The parallel-axis theorem expresses it in terms of the inertia tensor about the
+centre of mass.
 
 ## References
 - Landau and Lifshitz, Mechanics, page 100, Section 32
 -/
 
 @[expose] public section
 
-open Manifold
+open Manifold InnerProductSpace
 
 TODO "The definition of a rigid body is currently defined via linear maps
   from the space of smooth functions to ℝ. When possible, it should be change
@@ -91,6 +96,87 @@ lemma rho_coord_sub_centerOfMass {d : ℕ} (R : RigidBody d) (h : R.mass ≠ 0)
   rw [hsplit, map_sub, map_smul, R.rho_one, hcoord, smul_eq_mul]
   ring
 
+/-- The inertia tensor of the rigid body about the point `p`, that is with the moments taken
+about `p` rather than about the origin of the reference frame. Taking `p` to be the origin
+recovers `RigidBody.inertiaTensor`, see `RigidBody.inertiaTensorAbout_zero`. -/
+noncomputable def inertiaTensorAbout {d : ℕ} (R : RigidBody d) (p : Space d) :
+    Matrix (Fin d) (Fin d) ℝ := fun i j =>
+  R.ρ (cmap (fun y => (if i = j then 1 else 0) * ∑ k : Fin d, (y k - p k) ^ 2
+    - (y i - p i) * (y j - p j)) (by fun_prop))
+
+/-- The `(i, j)` entry of the inertia tensor about `p` as a moment of the mass distribution. -/
+lemma inertiaTensorAbout_apply {d : ℕ} (R : RigidBody d) (p : Space d) (i j : Fin d) :
+    R.inertiaTensorAbout p i j = R.ρ (cmap (fun y => (if i = j then 1 else 0) *
+      ∑ k : Fin d, (y k - p k) ^ 2 - (y i - p i) * (y j - p j)) (by fun_prop)) := rfl
+
+/-- The inertia tensor about the origin is the inertia tensor. -/
+@[simp]
+lemma inertiaTensorAbout_zero {d : ℕ} (R : RigidBody d) :
+    R.inertiaTensorAbout 0 = R.inertiaTensor := by
+  ext i j
+  simp [inertiaTensorAbout, inertiaTensor, cmap]
+
+/-- The inertia tensor about any point is symmetric. -/
+lemma inertiaTensorAbout_symmetric {d : ℕ} (R : RigidBody d) (p : Space d) (i j : Fin d) :
+    R.inertiaTensorAbout p i j = R.inertiaTensorAbout p j i := by
+  simp only [inertiaTensorAbout, eq_comm, mul_comm]
+
+/-- The integrand of the inertia tensor about `p` splits into the integrand about the centre of
+mass, terms linear in the centred coordinates `y − c`, and a constant built from the
+displacement `c − p`. -/
+lemma inertiaTensorAbout_integrand_split {d : ℕ} (R : RigidBody d) (p : Space d) (i j : Fin d) :
+    cmap (fun y => (if i = j then 1 else 0) * ∑ k : Fin d, (y k - p k) ^ 2
+        - (y i - p i) * (y j - p j)) (by fun_prop)
+      = cmap (fun y => (if i = j then 1 else 0) * ∑ k : Fin d, (y k - R.centerOfMass k) ^ 2
+            - (y i - R.centerOfMass i) * (y j - R.centerOfMass j)) (by fun_prop)
+        + ∑ k : Fin d, (2 * (if i = j then 1 else 0) * (R.centerOfMass k - p k)) •
+            cmap (fun y => y k - R.centerOfMass k) (by fun_prop)
+        - (R.centerOfMass j - p j) • cmap (fun y => y i - R.centerOfMass i) (by fun_prop)
+        - (R.centerOfMass i - p i) • cmap (fun y => y j - R.centerOfMass j) (by fun_prop)
+        + ((if i = j then 1 else 0) * ∑ k : Fin d, (R.centerOfMass k - p k) ^ 2
+            - (R.centerOfMass i - p i) * (R.centerOfMass j - p j)) •
+          (1 : C^⊤⟮𝓘(ℝ, Space d), Space d; 𝓘(ℝ, ℝ), ℝ⟯) := by
+  ext y
+  simp only [cmap_apply, ContMDiffMap.coe_add, ContMDiffMap.coe_sub, ContMDiffMap.coe_smul,
+    ContMDiffMap.coe_one, Pi.add_apply, Pi.sub_apply, Pi.smul_apply, Pi.one_apply]
+  rw [← ContMDiffMap.coeFnAddMonoidHom_apply, map_sum, Finset.sum_apply]
+  simp only [ContMDiffMap.coeFnAddMonoidHom_apply, ContMDiffMap.coe_smul, Pi.smul_apply,
+    cmap_apply, smul_eq_mul]
+  have hlin (a : ℝ) : ∑ k : Fin d, a * (R.centerOfMass k - p k) * (y k - R.centerOfMass k)
+      = a * ∑ k : Fin d, (R.centerOfMass k - p k) * (y k - R.centerOfMass k) := by
+    rw [Finset.mul_sum]
+    exact Finset.sum_congr rfl fun k _ => by ring
+  have hsum : ∑ k : Fin d, (y k - p k) ^ 2
+      = ∑ k : Fin d, ((y k - R.centerOfMass k) ^ 2
+        + 2 * (R.centerOfMass k - p k) * (y k - R.centerOfMass k)
+        + (R.centerOfMass k - p k) ^ 2) := Finset.sum_congr rfl fun k _ => by ring
+  rw [hsum, Finset.sum_add_distrib, Finset.sum_add_distrib, hlin, hlin]
+  ring
+
+/-- **The parallel-axis theorem**: the inertia tensor about a point `p` is the inertia tensor
+about the centre of mass `c` plus the inertia tensor about `p` of a point particle carrying the
+total mass of the body and sitting at `c`, that is `M (|c − p|² 1 − (c − p) ⊗ (c − p))`. -/
+theorem inertiaTensorAbout_eq_centerOfMass_add_pointMass {d : ℕ} (R : RigidBody d) (h : R.mass ≠ 0)
+    (p : Space d) :
+    R.inertiaTensorAbout p = R.inertiaTensorAbout R.centerOfMass
+      + R.mass • (⟪R.centerOfMass - p, R.centerOfMass - p⟫_ℝ • (1 : Matrix (Fin d) (Fin d) ℝ)
+        - Matrix.vecMulVec ⇑(R.centerOfMass - p) ⇑(R.centerOfMass - p)) := by
+  ext i j
+  rw [inertiaTensorAbout_apply, inertiaTensorAbout_integrand_split R p i j, map_add, map_sub,
+    map_sub, map_add, map_sum, ← inertiaTensorAbout_apply]
+  simp only [map_smul, R.rho_coord_sub_centerOfMass h, smul_eq_mul, mul_zero,
+    Finset.sum_const_zero, R.rho_one, Matrix.add_apply, Matrix.smul_apply, Matrix.sub_apply,
+    Matrix.one_apply, Matrix.vecMulVec_apply, Space.inner_eq_sum, Space.sub_apply, pow_two]
+  ring
+
+/-- The inertia tensor about the origin is the inertia tensor about the centre of mass `c` plus
+`M (|c|² 1 − c ⊗ c)`. -/
+lemma inertiaTensor_eq_centerOfMass_add_pointMass {d : ℕ} (R : RigidBody d) (h : R.mass ≠ 0) :
+    R.inertiaTensor = R.inertiaTensorAbout R.centerOfMass
+      + R.mass • (⟪R.centerOfMass, R.centerOfMass⟫_ℝ • (1 : Matrix (Fin d) (Fin d) ℝ)
+        - Matrix.vecMulVec ⇑R.centerOfMass ⇑R.centerOfMass) := by
+  simpa using R.inertiaTensorAbout_eq_centerOfMass_add_pointMass h 0
+
 /-- One can describe the motion of rigid body with a fixed (inertial) coordinate system (X,Y,Z)
     and a moving system (x₁,x₂,x₃) rigidly attached to the body. -/
 informal_definition coordinate_system where
@@ -157,12 +243,6 @@ informal_lemma kinetic_energy_decomposition where
   tag := "LL32-TK"
   deps := [``RigidBody]
 
-/-- If I_O is the inertia tensor about a point O, then the inertia tensor about a parallel point O'
-    displaced by a is I_{O'} = I_O + M(|a|² 1 − a ⊗ a). This is the parallel-axis theorem. -/
-informal_lemma parallel_axis_theorem where
-  tag := "LL32-PA"
-  deps := [``RigidBody]
-
 /-- Because the inertia tensor is real symmetric, there exists an orthonormal basis of principal
     axes in which it is diagonal. The corresponding directions are the principal axes of inertia. -/
 informal_definition principal_axes_of_inertia where