Some improvements to the computations of collisions

mathematica/mathematica.hpp

@@ -231,6 +231,10 @@ template<typename Tuple,
 std::string ToMathematica(Tuple const& tuple,
                           OptionalExpressIn express_in = std::nullopt);
 
+template<typename Scalar, typename OptionalExpressIn = std::nullopt_t>
+std::string ToMathematica(UnboundedMatrix<Scalar> const& matrix,
+                          OptionalExpressIn express_in = std::nullopt);
+
 template<typename Scalar, typename OptionalExpressIn = std::nullopt_t>
 std::string ToMathematica(UnboundedLowerTriangularMatrix<Scalar> const& matrix,
                           OptionalExpressIn express_in = std::nullopt);

mathematica/mathematica_body.hpp

@@ -460,6 +460,22 @@ std::string ToMathematica(Tuple const& tuple, OptionalExpressIn express_in) {
   return RawApply("List", expressions);
 }
 
+template<typename Scalar, typename OptionalExpressIn>
+std::string ToMathematica(UnboundedMatrix<Scalar> const& matrix,
+                          OptionalExpressIn express_in) {
+  std::vector<std::string> rows;
+  rows.reserve(matrix.rows());
+  for (int i = 0; i < matrix.rows(); ++i) {
+    std::vector<std::string> row;
+    row.reserve(matrix.columns());
+    for (int j = 0; j < matrix.columns(); ++j) {
+      row.push_back(ToMathematica(matrix(i, j), express_in));
+    }
+    rows.push_back(RawApply("List", row));
+  }
+  return RawApply("List", rows);
+}
+
 template<typename Scalar, typename OptionalExpressIn>
 std::string ToMathematica(UnboundedLowerTriangularMatrix<Scalar> const& matrix,
                           OptionalExpressIn express_in) {

numerics/approximation_body.hpp

@@ -79,8 +79,16 @@ template<int N, typename Argument, typename Function>
           f, a, b, max_error, fₖ, aⱼ);
     }
   }
+
+  // Unlike [Boy13], section 4, we return the polynomial of the lower degree
+  // that is within the |max_error| bound (that is, the one of degree N / 2).
+  // Even though we computed the polynomial of degree N, returning it would
+  // impose an unnecessary cost on the client (e.g., more costly evaluation). If
+  // a client wants a more precise approximation, they just need to give a
+  // smaller |max_error|.
   std::vector<Value<Argument, Function>> coefficients;
-  std::copy(aⱼ.begin(), aⱼ.end(), std::back_inserter(coefficients));
+  std::copy(previous_aⱼ.begin(), previous_aⱼ.end(),
+            std::back_inserter(coefficients));
   return ЧебышёвSeries<Value<Argument, Function>, Argument>(
       coefficients, a, b);
 }

numerics/approximation_test.cpp

@@ -26,10 +26,10 @@ TEST(ApproximationTest, SinInverse) {
                                            /*lower_bound=*/0.1,
                                            /*upper_bound=*/10,
                                            /*max_error=*/1e-6);
-  EXPECT_EQ(256, approximation.degree());
+  EXPECT_EQ(128, approximation.degree());
   for (double x = 0.1; x < 10.0; x += 0.01) {
     EXPECT_THAT(approximation.Evaluate(x),
-                AbsoluteErrorFrom(f(x), AllOf(Lt(1.1e-13), Ge(0))));
+                AbsoluteErrorFrom(f(x), AllOf(Lt(3.0e-6), Ge(0))));
   }
 }
 
@@ -40,10 +40,10 @@ TEST(ApproximationTest, Exp) {
                                            /*lower_bound=*/0.01,
                                            /*upper_bound=*/3,
                                            /*max_error=*/1e-6);
-  EXPECT_EQ(32, approximation.degree());
+  EXPECT_EQ(16, approximation.degree());
   for (double x = 0.01; x < 3; x += 0.01) {
     EXPECT_THAT(approximation.Evaluate(x),
-                AbsoluteErrorFrom(f(x), AllOf(Lt(1.8e-14), Ge(0))));
+                AbsoluteErrorFrom(f(x), AllOf(Lt(7.2e-15), Ge(0))));
   }
 }
 

numerics/matrix_computations_body.hpp

@@ -6,8 +6,10 @@
 #include <limits>
 #include <utility>
 
+#include "absl/container/btree_set.h"
 #include "base/tags.hpp"
 #include "numerics/fixed_arrays.hpp"
+#include "numerics/root_finders.hpp"
 #include "numerics/transposed_view.hpp"
 #include "numerics/unbounded_arrays.hpp"
 #include "quantities/elementary_functions.hpp"
@@ -20,6 +22,7 @@ namespace internal {
 
 using namespace principia::base::_tags;
 using namespace principia::numerics::_fixed_arrays;
+using namespace principia::numerics::_root_finders;
 using namespace principia::numerics::_transposed_view;
 using namespace principia::numerics::_unbounded_arrays;
 using namespace principia::quantities::_elementary_functions;
@@ -335,6 +338,21 @@ void PostMultiply(Matrix& A, JacobiRotation const& J) {
   }
 }
 
+template<typename Matrix>
+absl::btree_set<typename Matrix::Scalar> Compute2By2Eigenvalues(
+    BlockView<Matrix> const& block) {
+  static constexpr typename Matrix::Scalar zero{};
+  // TODO(phl): Would |SymmetricSchurDecomposition2By2| work to shoot one zero
+  // (even though the |block| is not symmetric)?
+  auto const& a = block(0, 0);
+  auto const& b = block(0, 1);
+  auto const& c = block(1, 0);
+  auto const& d = block(1, 1);
+  auto const solutions = SolveQuadraticEquation(zero, a * d - b * c, -a - d, 1);
+  return absl::btree_set<typename Matrix::Scalar>(solutions.begin(),
+                                                  solutions.end());
+}
+
 // [GV13] algorithm 7.5.1.
 template<typename Scalar, typename Matrix>
 void FrancisQRStep(Matrix& H) {
@@ -469,6 +487,7 @@ template<typename Scalar>
 struct RealSchurDecompositionGenerator<UnboundedMatrix<Scalar>> {
   struct Result {
     UnboundedMatrix<Scalar> T;
+    absl::btree_set<double> real_eigenvalues;
   };
 };
 
@@ -477,6 +496,7 @@ struct RealSchurDecompositionGenerator<
     FixedMatrix<Scalar, dimension, dimension>> {
   struct Result {
     FixedMatrix<Scalar, dimension, dimension> T;
+    absl::btree_set<double> real_eigenvalues;
   };
 };
 
@@ -855,8 +875,40 @@ RealSchurDecomposition(Matrix const& A, double const ε) {
     FrancisQRStep<Scalar>(H₂₂);
   }
 
-  typename RealSchurDecompositionGenerator<Matrix>::Result result{.T = H};
-  return result;
+  // Find the real eigenvalues.  Note that they may be part of a 2×2 block which
+  // happens to have real roots.
+  absl::btree_set<Scalar> real_eigenvalues;
+  for (int i = 0; i < H.rows();) {
+    int first_in_block = 0;
+    if (i == H.rows() - 1) {
+      if (H(i, i - 1) == 0) {
+        real_eigenvalues.insert(H(i, i));
+      }
+      break;
+    } else if (H(i + 1, i) == 0) {
+      real_eigenvalues.insert(H(i, i));
+      ++i;
+      continue;
+    } else {
+      first_in_block = i;
+    }
+    auto const block = BlockView<Matrix>{.matrix = H,
+                                         .first_row = first_in_block,
+                                         .last_row = first_in_block + 1,
+                                         .first_column = first_in_block,
+                                         .last_column = first_in_block + 1};
+    auto const block_real_eigenvalues = Compute2By2Eigenvalues(block);
+    real_eigenvalues.insert(block_real_eigenvalues.begin(),
+                            block_real_eigenvalues.end());
+
+    // The 2×2 block is processed on its first index, so the next index is
+    // skipped.
+    i += 2;
+  }
+
+  return typename RealSchurDecompositionGenerator<Matrix>::Result{
+      .T = std::move(H),
+      .real_eigenvalues = std::move(real_eigenvalues)};
 }
 
 template<typename Matrix>

numerics/matrix_computations_test.cpp

@@ -151,10 +151,11 @@ TYPED_TEST(MatrixComputationsTest, RealSchurDecomposition) {
                     8, -9, -2,  4});
   auto s4 = RealSchurDecomposition(m4, 1e-6);
   // Only check the real eigenvalues.
-  EXPECT_THAT(s4.T(2, 2),
-              RelativeErrorFrom(8.8004352424313246181, IsNear(6.0e-7_(1))));
-  EXPECT_THAT(s4.T(3, 3),
-              RelativeErrorFrom(6.2103405225078473234, IsNear(8.4e-7_(1))));
+  EXPECT_THAT(
+      s4.real_eigenvalues,
+      ElementsAre(
+          RelativeErrorFrom(6.2103405225078473234, IsNear(8.4e-7_(1))),
+          RelativeErrorFrom(8.8004352424313246181, IsNear(6.0e-7_(1)))));
 }
 
 TYPED_TEST(MatrixComputationsTest, ClassicalJacobi) {

numerics/чебышёв_series.hpp

@@ -98,7 +98,7 @@ class ЧебышёвSeries final {
       Argument const& argument) const;
 
   // Returns the Frobenius companion matrix suitable for the Чебышёв basis.
-  UnboundedMatrix<Value> FrobeniusCompanionMatrix() const;
+  UnboundedMatrix<double> FrobeniusCompanionMatrix() const;
 
   void WriteToMessage(not_null<serialization::ЧебышёвSeries*> message) const;
   static ЧебышёвSeries ReadFromMessage(

numerics/чебышёв_series_body.hpp

@@ -259,10 +259,10 @@ Derivative<Value, Argument> ЧебышёвSeries<Value, Argument>::EvaluateDeriv
 }
 
 template<typename Value, typename Argument>
-UnboundedMatrix<Value>
+UnboundedMatrix<double>
 ЧебышёвSeries<Value, Argument>::FrobeniusCompanionMatrix() const {
   int const N = degree();
-  UnboundedMatrix<Value> A(N, N, uninitialized);
+  UnboundedMatrix<double> A(N, N, uninitialized);
 
   // j = 1.
   for (int k = 1; k <= N; ++k) {

numerics/чебышёв_series_test.cpp

@@ -5,7 +5,9 @@
 #include "geometry/instant.hpp"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
+#include "numerics/matrix_computations.hpp"
 #include "numerics/unbounded_arrays.hpp"
+#include "quantities/elementary_functions.hpp"
 #include "quantities/named_quantities.hpp"
 #include "quantities/quantities.hpp"
 #include "quantities/si.hpp"
@@ -14,11 +16,14 @@
 namespace principia {
 namespace numerics {
 
+using ::testing::ElementsAre;
 using namespace principia::astronomy::_frames;
 using namespace principia::geometry::_grassmann;
 using namespace principia::geometry::_instant;
+using namespace principia::numerics::_matrix_computations;
 using namespace principia::numerics::_unbounded_arrays;
 using namespace principia::numerics::_чебышёв_series;
+using namespace principia::quantities::_elementary_functions;
 using namespace principia::quantities::_named_quantities;
 using namespace principia::quantities::_quantities;
 using namespace principia::quantities::_si;
@@ -127,14 +132,17 @@ TEST_F(ЧебышёвSeriesTest, T2Double) {
 TEST_F(ЧебышёвSeriesTest, FrobeniusCompanionMatrix) {
   ЧебышёвSeries<double, Instant> series({-2, 3, 5, 6}, t_min_, t_max_);
   auto const matrix = series.FrobeniusCompanionMatrix();
-  // Checked with Mathematica that the eigenvalues of this matrix are the zeroes
-  // of the above series.
   EXPECT_THAT(matrix,
               AlmostEquals(UnboundedMatrix<double>(
                   {0.0,             1.0,         0.0,
                    1.0 / 2.0,       0.0,   1.0 / 2.0,
                    1.0 / 6.0, 1.0 / 4.0, -5.0 / 12.0}),
                   0));
+  auto const matrix_schur_decomposition = RealSchurDecomposition(matrix, 1e-16);
+  EXPECT_THAT(matrix_schur_decomposition.real_eigenvalues,
+              ElementsAre(AlmostEquals((1.0 - Sqrt(337.0)) / 24.0, 4),
+                          AlmostEquals(-0.5, 1),
+                          AlmostEquals((1.0 + Sqrt(337.0)) / 24.0, 2)));
 }
 
 TEST_F(ЧебышёвSeriesTest, X6Vector) {