Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 585223bd99aa
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: b5dd83ac1c51
Choose a head ref
  • 2 commits
  • 3 files changed
  • 1 contributor

Commits on Feb 18, 2020

  1. Copy the full SHA
    a398933 View commit details
  2. Merge pull request #2473 from pleroy/Diagonalize

    Deal with a degenerate case when diagonalizing a form
    pleroy authored Feb 18, 2020
    Copy the full SHA
    b5dd83a View commit details
Showing with 38 additions and 3 deletions.
  1. +2 −1 geometry/symmetric_bilinear_form.hpp
  2. +14 −2 geometry/symmetric_bilinear_form_body.hpp
  3. +22 −0 geometry/symmetric_bilinear_form_test.cpp
3 changes: 2 additions & 1 deletion geometry/symmetric_bilinear_form.hpp
Original file line number Diff line number Diff line change
@@ -66,7 +66,8 @@ class SymmetricBilinearForm {

// The eigensystem for a form is described by (1) the form in its eigenbasis,
// which gives the eigenvalues; and (2) a rotation from the current basis to
// the eigenbasis, which gives the eigenvectors.
// the eigenbasis, which gives the eigenvectors. The eigenvalues are in
// increasing order.
template<typename Eigenframe>
struct Eigensystem {
SymmetricBilinearForm<Scalar, Eigenframe, Multivector> form;
16 changes: 14 additions & 2 deletions geometry/symmetric_bilinear_form_body.hpp
Original file line number Diff line number Diff line change
@@ -129,6 +129,7 @@ template<typename Eigenframe>
typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
template Eigensystem<Eigenframe>
SymmetricBilinearForm<Scalar, Frame, Multivector>::Diagonalize() const {
Scalar const zero;
R3x3Matrix<Scalar> const& A = matrix_;
auto const I = R3x3Matrix<double>::Identity();

@@ -138,20 +139,31 @@ typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
Scalar const q = A.Trace() / 3;
R3x3Matrix<Scalar> const A_minus_qI = A - q * I;
Scalar const p = Sqrt((A_minus_qI * A_minus_qI).Trace() / 6);

if (p == zero) {
// A is very close to q * I.
SymmetricBilinearForm<Scalar, Eigenframe, Multivector> form(
R3x3Matrix<Scalar>({q, zero, zero},
{zero, q, zero},
{zero, zero, q}));
return {form, Rotation<Eigenframe, Frame>::Identity()};
}

R3x3Matrix<double> const B = A_minus_qI / p;
double const det_B = B.Determinant();
Angle const θ = ArcCos(det_B * 0.5);
Angle const θ = ArcCos(std::clamp(det_B, -2.0, 2.0) * 0.5);
double const β₀ = 2 * Cos(θ / 3);
double const β₁ = 2 * Cos((θ + 2 * π * Radian) / 3);
double const β₂ = 2 * Cos((θ + 4 * π * Radian) / 3);
std::array<Scalar, 3> αs = {p * β₀ + q, p * β₁ + q, p * β₂ + q};
// We expect αs[1] <= αs[2] <= αs[0] here, but sorting ensures that we are
// correct irrespective of numerical errors.
std::sort(αs.begin(), αs.end());
Scalar const& α₀ = αs[0];
Scalar const& α₁ = αs[1];
Scalar const& α₂ = αs[2];

// The form in its diagonal basis.
Scalar const zero;
SymmetricBilinearForm<Scalar, Eigenframe, Multivector> form(
R3x3Matrix<Scalar>({α₀, zero, zero},
{zero, α₁, zero},
22 changes: 22 additions & 0 deletions geometry/symmetric_bilinear_form_test.cpp
Original file line number Diff line number Diff line change
@@ -359,6 +359,28 @@ TEST_F(SymmetricBilinearFormTest, Diagonalize) {
VanishesBefore(1, 2),
AlmostEquals(1, 0)));
}

// A degenerate case.
{
auto const f = MakeSymmetricBilinearForm<World>(R3x3Matrix<double>(
{2, 0, 0},
{0, 2, 0},
{0, 0, 2}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

Vector<double, Eigenworld> const e₀({1, 0, 0});
Vector<double, Eigenworld> const e₁({0, 1, 0});
Vector<double, Eigenworld> const e₂({0, 0, 1});
EXPECT_THAT(f_eigensystem.form(e₀, e₀), AlmostEquals(2 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₀, e₁), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₀, e₂), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₁, e₀), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₁, e₁), AlmostEquals(2 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₁, e₂), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₀), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₁), AlmostEquals(0 * Metre, 0));
EXPECT_THAT(f_eigensystem.form(e₂, e₂), AlmostEquals(2 * Metre, 0));
}
}

} // namespace internal_symmetric_bilinear_form