From bfd843b249e24dd2c638997aaaf9baebdc0cb7c5 Mon Sep 17 00:00:00 2001 From: pleroy Date: Fri, 22 Apr 2022 22:00:05 +0200 Subject: [PATCH 1/3] Use Brent's method instead of bisection for computing apsides and nodes. --- physics/apsides_body.hpp | 9 ++++----- physics/ephemeris_body.hpp | 12 ++++++------ physics/kepler_orbit_body.hpp | 14 +++++++------- 3 files changed, 17 insertions(+), 18 deletions(-) diff --git a/physics/apsides_body.hpp b/physics/apsides_body.hpp index 575d642ded..c49f504e73 100644 --- a/physics/apsides_body.hpp +++ b/physics/apsides_body.hpp @@ -19,7 +19,7 @@ using geometry::Barycentre; using geometry::Instant; using geometry::Position; using geometry::Sign; -using numerics::Bisect; +using numerics::Brent; using numerics::Hermite3; using quantities::IsFinite; using quantities::Length; @@ -176,10 +176,9 @@ absl::Status ComputeNodes( node_time = Barycentre({*previous_time, time}, {z, -*previous_z}); } else { - // The normal case, find the intersection with z = 0 using bisection. - // TODO(egg): Bisection on a polynomial seems daft; we should have - // Newton's method. - node_time = Bisect( + // The normal case, find the intersection with z = 0 using Brent's + // method. + node_time = Brent( [&z_approximation](Instant const& t) { return z_approximation.Evaluate(t); }, diff --git a/physics/ephemeris_body.hpp b/physics/ephemeris_body.hpp index 6d3f6e0324..0ca4ad0bac 100644 --- a/physics/ephemeris_body.hpp +++ b/physics/ephemeris_body.hpp @@ -49,7 +49,7 @@ using integrators::ExplicitSecondOrderOrdinaryDifferentialEquation; using integrators::IntegrationProblem; using integrators::Integrator; using integrators::methods::Fine1987RKNG34; -using numerics::Bisect; +using numerics::Brent; using numerics::DoublePrecision; using numerics::Hermite3; using quantities::Abs; @@ -721,11 +721,11 @@ void Ephemeris::ComputeApsides(not_null const body1, CHECK(previous_time); // The derivative of |squared_distance| changed sign. Find its zero by - // bisection, this is the time of the apsis. Then compute the apsis and - // append it to one of the output trajectories. - Instant const apsis_time = Bisect(evaluate_square_distance_derivative, - *previous_time, - time); + // Brent's method, this is the time of the apsis. Then compute the apsis + // and append it to one of the output trajectories. + Instant const apsis_time = Brent(evaluate_square_distance_derivative, + *previous_time, + time); DegreesOfFreedom const apsis1_degrees_of_freedom = body1_trajectory->EvaluateDegreesOfFreedomLocked(apsis_time); DegreesOfFreedom const apsis2_degrees_of_freedom = diff --git a/physics/kepler_orbit_body.hpp b/physics/kepler_orbit_body.hpp index 908fd4b258..4cabaccb8c 100644 --- a/physics/kepler_orbit_body.hpp +++ b/physics/kepler_orbit_body.hpp @@ -29,7 +29,7 @@ using geometry::Sign; using geometry::Vector; using geometry::Velocity; using geometry::Wedge; -using numerics::Bisect; +using numerics::Brent; using quantities::Abs; using quantities::ArcCos; using quantities::ArcCosh; @@ -642,9 +642,9 @@ void KeplerOrbit::CompleteAnomalies(KeplerianElements& elements) { (eccentric_anomaly - e * Sin(eccentric_anomaly) * Radian); }; Angle const eccentric_anomaly = e == 0 ? *mean_anomaly - : Bisect(kepler_equation, - *mean_anomaly - e * Radian, - *mean_anomaly + e * Radian); + : Brent(kepler_equation, + *mean_anomaly - e * Radian, + *mean_anomaly + e * Radian); true_anomaly = 2 * ArcTan(Sqrt(1 + e) * Sin(eccentric_anomaly / 2), Sqrt(1 - e) * Cos(eccentric_anomaly / 2)); hyperbolic_mean_anomaly = NaN; @@ -657,9 +657,9 @@ void KeplerOrbit::CompleteAnomalies(KeplerianElements& elements) { hyperbolic_eccentric_anomaly); }; Angle const hyperbolic_eccentric_anomaly = - Bisect(hyperbolic_kepler_equation, - 0 * Radian, - *hyperbolic_mean_anomaly / (e - 1)); + Brent(hyperbolic_kepler_equation, + 0 * Radian, + *hyperbolic_mean_anomaly / (e - 1)); true_anomaly = 2 * ArcTan(Sqrt(e + 1) * Sinh(hyperbolic_eccentric_anomaly / 2), Sqrt(e - 1) * Cosh(hyperbolic_eccentric_anomaly / 2)); From e1808f71917359e8597245f693470dc66c8a2a7b Mon Sep 17 00:00:00 2001 From: pleroy Date: Fri, 22 Apr 2022 22:33:28 +0200 Subject: [PATCH 2/3] Tolerances. --- astronomy/ksp_resonance_test.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/astronomy/ksp_resonance_test.cpp b/astronomy/ksp_resonance_test.cpp index 7fc139449b..d85a3fad2d 100644 --- a/astronomy/ksp_resonance_test.cpp +++ b/astronomy/ksp_resonance_test.cpp @@ -287,13 +287,13 @@ TEST_F(KSPResonanceTest, MSVC_ONLY_TEST(Stock)) { ephemeris->t_max() - 2 * longest_joolian_period_); EXPECT_THAT(RelativeError(periods_at_mid_term.at(laythe_), expected_periods_.at(laythe_)), - IsNear(0.874_⑴)); + IsNear(0.121_⑴)); EXPECT_THAT(periods_at_mid_term.at(vall_), Eq(Infinity