Merge pull request #3621 from pleroy/Jacobian2

A more elegant computation of the Jacobian

physics/ephemeris.hpp

@@ -200,7 +200,7 @@ class Ephemeris {
 
   // Returns the Jacobian of the acceleration field exerted on the given |body|
   // by the rest of the system.
-  Jacobian<Frame> ComputeJacobianOnMassiveBody(
+  JacobianOfAcceleration<Frame> ComputeJacobianOnMassiveBody(
       not_null<MassiveBody const*> body,
       Instant const& t) const EXCLUDES(lock_);
 
@@ -325,7 +325,7 @@ class Ephemeris {
       std::size_t b2_begin,
       std::size_t b2_end,
       std::vector<Position<Frame>> const& positions,
-      std::vector<Jacobian<Frame>>& jacobians);
+      std::vector<JacobianOfAcceleration<Frame>>& jacobians);
 
   // Computes the accelerations between one body, |body1| (with index |b1| in
   // the |positions| and |accelerations| arrays) and the bodies |bodies2| (with

physics/ephemeris_body.hpp

@@ -451,12 +451,12 @@ absl::Status Ephemeris<Frame>::FlowWithFixedStep(
 }
 
 template<typename Frame>
-Jacobian<Frame> Ephemeris<Frame>::ComputeJacobianOnMassiveBody(
+JacobianOfAcceleration<Frame> Ephemeris<Frame>::ComputeJacobianOnMassiveBody(
     not_null<MassiveBody const*> body,
     Instant const& t) const {
   // NOTE(phl): This doesn't take high-order geopotential into account.
   std::vector<Position<Frame>> positions;
-  std::vector<Jacobian<Frame>> jacobians(bodies_.size());
+  std::vector<JacobianOfAcceleration<Frame>> jacobians(bodies_.size());
   int b1 = -1;
 
   // Evaluate the |positions|.  Locking is necessary to be able to call the
@@ -1094,40 +1094,30 @@ void Ephemeris<Frame>::ComputeJacobianByMassiveBodyOnMassiveBodies(
     std::size_t b2_begin,
     std::size_t b2_end,
     std::vector<Position<Frame>> const& positions,
-    std::vector<Jacobian<Frame>>& jacobians) {
+    std::vector<JacobianOfAcceleration<Frame>>& jacobians) {
   Position<Frame> const& position_of_b1 = positions[b1];
-  Jacobian<Frame>& jacobian_on_b1 = jacobians[b1];
+  JacobianOfAcceleration<Frame>& jacobian_on_b1 = jacobians[b1];
   GravitationalParameter const& μ1 = body1.gravitational_parameter();
   for (std::size_t b2 = b2_begin; b2 < b2_end; ++b2) {
-    Jacobian<Frame>& jacobian_on_b2 = jacobians[b2];
+    JacobianOfAcceleration<Frame>& jacobian_on_b2 = jacobians[b2];
     MassiveBody const& body2 = *bodies2[b2];
     GravitationalParameter const& μ2 = body2.gravitational_parameter();
 
     // A vector from the center of |b2| to the center of |b1|.
     Displacement<Frame> const Δq = position_of_b1 - positions[b2];
-    R3Element<Length> const& Δq_coordinates = Δq.coordinates();
-    Length const& Δq_x = Δq_coordinates.x;
-    Length const& Δq_y = Δq_coordinates.y;
-    Length const& Δq_z = Δq_coordinates.z;
-
-    Length const Δq_norm = Δq.Norm();
-    Square<Length> const Δq_norm² = Δq.Norm²();
-    auto const one_over_Δq_norm⁵ = 1 / (Δq_norm * Δq_norm² * Δq_norm²);
-
-    // Note that the matrix below is symmetric by construction (and that the
-    // parentheses matters for that purpose).
-    SymmetricBilinearForm<Square<Length>, Frame, Vector> const numerator(
-        {{+2 * Pow<2>(Δq_x) - Pow<2>(Δq_y) - Pow<2>(Δq_z),
-          3 * (Δq_x * Δq_y),
-          3 * (Δq_x * Δq_z)},
-         {3 * (Δq_y * Δq_x),
-          -Pow<2>(Δq_x) + 2 * Pow<2>(Δq_y) - Pow<2>(Δq_z),
-          3 * (Δq_y * Δq_z)},
-         {3 * (Δq_z * Δq_x),
-          3 * (Δq_z * Δq_y),
-          -Pow<2>(Δq_x) - Pow<2>(Δq_y) + 2 * Pow<2>(Δq_z)}});
-    jacobians[b2] -= μ1 * numerator * one_over_Δq_norm⁵;
-    jacobians[b1] += μ2 * numerator * one_over_Δq_norm⁵;
+
+    Square<Length> const Δq² = Δq.Norm²();
+    Length const Δq_norm = Sqrt(Δq²);
+    Cube<Length> const Δq_norm³ = Δq² * Δq_norm;
+    auto const Δq_norm⁵ = Δq_norm³ * Δq²;
+
+    auto const form = -InnerProductForm<Frame, Vector>() / Δq_norm³ +
+                      3 * SymmetricSquare(Δq) / Δq_norm⁵;
+
+    // Note that the Jacobian is independent from the sign of Δq, hence the +
+    // signs.
+    jacobian_on_b2 += μ1 * form;
+    jacobian_on_b1 += μ2 * form;
   }
 }
 

physics/tensors.hpp

@@ -20,12 +20,13 @@ using InertiaTensor = SymmetricBilinearForm<MomentOfInertia, Frame, Bivector>;
 
 // The type of the Jacobian of the gravitational acceleration field.
 template<typename Frame>
-using Jacobian = SymmetricBilinearForm<Inverse<Square<Time>>, Frame, Vector>;
+using JacobianOfAcceleration =
+    SymmetricBilinearForm<Inverse<Square<Time>>, Frame, Vector>;
 
 }  // namespace internal
 
 using internal::InertiaTensor;
-using internal::Jacobian;
+using internal::JacobianOfAcceleration;
 
 }  // namespace _tensors
 }  // namespace physics