From dd30c03b02fb95b88ee7010981efab3b7ade4825 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 15:30:43 +0100 Subject: [PATCH 1/4] Adaptive subdivision. --- numerics/approximation.hpp | 17 +++++- numerics/approximation_body.hpp | 55 ++++++++++++++++++-- numerics/approximation_test.cpp | 91 +++++++++++++++++++++++++++++++-- 3 files changed, 155 insertions(+), 8 deletions(-) diff --git a/numerics/approximation.hpp b/numerics/approximation.hpp index 9643eeecf8..b37c1b52a8 100644 --- a/numerics/approximation.hpp +++ b/numerics/approximation.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include #include "numerics/чебышёв_series.hpp" #include "quantities/named_quantities.hpp" @@ -16,7 +17,7 @@ using namespace principia::quantities::_named_quantities; template using Value = std::invoke_result_t; -// Returns a Чебышёв polynomial approximant of f over +// Returns a Чебышёв polynomial interpolant of f over // [lower_bound, upper_bound]. Stops if the absolute error is estimated to be // below |max_error| or if |max_degree| has been reached. If |error_estimate| // is nonnull, it receives the estimate of the L∞ error. @@ -28,8 +29,22 @@ template Difference> const& max_error, Difference>* error_estimate = nullptr); +// Returns an ordered vector of Чебышёв polynomial interpolants of f over, which +// togethor cover [lower_bound, upper_bound]. Subdivides the interval so that +// the degree of each approximant doesn't exceed |max_degree|. The final +// (estimated) absolute error is guaranteed to be below |max_error|. +template +std::vector<ЧебышёвSeries, Argument>> +AdaptiveЧебышёвPolynomialInterpolant( + Function const& f, + Argument const& lower_bound, + Argument const& upper_bound, + Difference> const& max_error, + Difference>* error_estimate = nullptr); + } // namespace internal +using internal::AdaptiveЧебышёвPolynomialInterpolant; using internal::ЧебышёвPolynomialInterpolant; } // namespace _approximation diff --git a/numerics/approximation_body.hpp b/numerics/approximation_body.hpp index a56393f5c5..2a4540757a 100644 --- a/numerics/approximation_body.hpp +++ b/numerics/approximation_body.hpp @@ -126,13 +126,60 @@ template FixedVector, 2> const aⱼ( {0.5 * (f_b + f_a), 0.5 * (f_b - f_a)}); - return ЧебышёвPolynomialInterpolantImplementation( + return ЧебышёвPolynomialInterpolantImplementation( f, a, b, max_error, fₖ, aⱼ, error_estimate); } +template +std::vector<ЧебышёвSeries, Argument>> +AdaptiveЧебышёвPolynomialInterpolant( + Function const& f, + Argument const& lower_bound, + Argument const& upper_bound, + Difference> const& max_error, + Difference>* error_estimate) { + // Try to build an interpolation over the entire interval. + Difference> full_error_estimate; + auto full_interpolant = ЧебышёвPolynomialInterpolant( + f, lower_bound, upper_bound, max_error, &full_error_estimate); + if (full_error_estimate <= max_error) { + // If the interpolant over the entire interval is within the desired error + // bound, return it. + if (error_estimate != nullptr) { + *error_estimate = full_error_estimate; + } + std::vector<ЧебышёвSeries, Argument>> + interpolants; + interpolants.emplace_back(std::move(full_interpolant)); + return interpolants; + } else { + // If the interpolant over the entire interval is not within the desired + // error bound, subdivide the interval. + Difference> upper_error_estimate; + Difference> lower_error_estimate; + auto const midpoint = + Barycentre(std::pair(lower_bound, upper_bound), std::pair(1, 1)); + auto lower_interpolants = + AdaptiveЧебышёвPolynomialInterpolant( + f, lower_bound, midpoint, max_error, &lower_error_estimate); + auto upper_interpolants = + AdaptiveЧебышёвPolynomialInterpolant( + f, midpoint, upper_bound, max_error, &upper_error_estimate); + std::vector<ЧебышёвSeries, Argument>> + all_interpolants; + std::move(lower_interpolants.begin(), + lower_interpolants.end(), + std::back_inserter(all_interpolants)); + std::move(upper_interpolants.begin(), + upper_interpolants.end(), + std::back_inserter(all_interpolants)); + if (error_estimate != nullptr) { + *error_estimate = std::max(lower_error_estimate, upper_error_estimate); + } + return all_interpolants; + } +} + } // namespace internal } // namespace _approximation } // namespace numerics diff --git a/numerics/approximation_test.cpp b/numerics/approximation_test.cpp index b87b1639c2..f698d46afe 100644 --- a/numerics/approximation_test.cpp +++ b/numerics/approximation_test.cpp @@ -2,9 +2,11 @@ #include "gmock/gmock.h" #include "gtest/gtest.h" +#include "numerics/чебышёв_series.hpp" #include "quantities/elementary_functions.hpp" #include "quantities/named_quantities.hpp" #include "quantities/si.hpp" +#include "testing_utilities/almost_equals.hpp" #include "testing_utilities/approximate_quantity.hpp" #include "testing_utilities/is_near.hpp" #include "testing_utilities/numerics_matchers.hpp" @@ -13,15 +15,20 @@ namespace principia { namespace numerics { namespace _approximation { +using ::testing::AllOf; +using ::testing::ElementsAre; +using ::testing::Ge; +using ::testing::Lt; +using ::testing::Property; +using ::testing::SizeIs; +using namespace principia::numerics::_чебышёв_series; using namespace principia::quantities::_elementary_functions; using namespace principia::quantities::_named_quantities; using namespace principia::quantities::_si; +using namespace principia::testing_utilities::_almost_equals; using namespace principia::testing_utilities::_approximate_quantity; using namespace principia::testing_utilities::_is_near; using namespace principia::testing_utilities::_numerics_matchers; -using ::testing::AllOf; -using ::testing::Ge; -using ::testing::Lt; TEST(ApproximationTest, SinInverse) { auto const f = [](double const x) { return Sin(1 * Radian / x); }; @@ -58,6 +65,84 @@ TEST(ApproximationTest, Exp) { } } +TEST(ApproximationTest, AdaptiveSinInverse) { + auto const f = [](double const x) { return Sin(1 * Radian / x); }; + double error_estimate; + auto const approximations = + AdaptiveЧебышёвPolynomialInterpolant<8>(f, + /*lower_bound=*/0.1, + /*upper_bound=*/10.0, + /*max_error=*/1e-6, + &error_estimate); + EXPECT_THAT(error_estimate, IsNear(7.1e-7_(1))); + EXPECT_THAT(approximations, SizeIs(11)); + EXPECT_THAT( + approximations, + ElementsAre(AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 512, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 512, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 256, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 256, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 128, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 128, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 64, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 64, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 32, 1)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 32, 1)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 * 3 / 64, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 * 3 / 64, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 16, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 16, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 8, 1)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 8, 1)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 4, 1)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 4, 1)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(0.1 + 9.9 / 2, 0)), + Property(&ЧебышёвSeries::degree, 8)), + AllOf(Property(&ЧебышёвSeries::lower_bound, + AlmostEquals(0.1 + 9.9 / 2, 0)), + Property(&ЧебышёвSeries::upper_bound, + AlmostEquals(10.0, 0)), + Property(&ЧебышёвSeries::degree, 8)))); + for (auto const& approximation : approximations) { + for (double x = approximation.lower_bound(); + x < approximation.upper_bound(); + x += 0.01) { + EXPECT_THAT(approximation.Evaluate(x), + AbsoluteErrorFrom(f(x), AllOf(Lt(6.2e-7), Ge(2.7e-17)))); + } + } +} + } // namespace _approximation } // namespace numerics } // namespace principia From 6ca8698f285f86ad0060ab479f79f5da188e82a0 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 15:47:54 +0100 Subject: [PATCH 2/4] Readying. --- numerics/approximation.hpp | 2 +- numerics/approximation_test.cpp | 26 +++++++++++++------------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/numerics/approximation.hpp b/numerics/approximation.hpp index b37c1b52a8..0a1ff07dc5 100644 --- a/numerics/approximation.hpp +++ b/numerics/approximation.hpp @@ -30,7 +30,7 @@ template Difference>* error_estimate = nullptr); // Returns an ordered vector of Чебышёв polynomial interpolants of f over, which -// togethor cover [lower_bound, upper_bound]. Subdivides the interval so that +// together cover [lower_bound, upper_bound]. Subdivides the interval so that // the degree of each approximant doesn't exceed |max_degree|. The final // (estimated) absolute error is guaranteed to be below |max_error|. template diff --git a/numerics/approximation_test.cpp b/numerics/approximation_test.cpp index f698d46afe..0d6b448499 100644 --- a/numerics/approximation_test.cpp +++ b/numerics/approximation_test.cpp @@ -33,17 +33,17 @@ using namespace principia::testing_utilities::_numerics_matchers; TEST(ApproximationTest, SinInverse) { auto const f = [](double const x) { return Sin(1 * Radian / x); }; double error_estimate; - auto const approximation = + auto const interpolant = ЧебышёвPolynomialInterpolant<128>(f, /*lower_bound=*/0.1, /*upper_bound=*/10.0, /*max_error=*/1e-6, &error_estimate); - EXPECT_EQ(128, approximation.degree()); + EXPECT_EQ(128, interpolant.degree()); // Didn't reach the desired accuracy. EXPECT_THAT(error_estimate, IsNear(4.9e-6_(1))); for (double x = 0.1; x < 10.0; x += 0.01) { - EXPECT_THAT(approximation.Evaluate(x), + EXPECT_THAT(interpolant.Evaluate(x), AbsoluteErrorFrom(f(x), AllOf(Lt(3.0e-6), Ge(0)))); } } @@ -51,16 +51,16 @@ TEST(ApproximationTest, SinInverse) { TEST(ApproximationTest, Exp) { auto const f = [](double const x) { return std::exp(x); }; double error_estimate; - auto const approximation = + auto const interpolant = ЧебышёвPolynomialInterpolant<128>(f, /*lower_bound=*/0.01, /*upper_bound=*/3.0, /*max_error=*/1e-6, &error_estimate); - EXPECT_EQ(16, approximation.degree()); + EXPECT_EQ(16, interpolant.degree()); EXPECT_THAT(error_estimate, IsNear(4.7e-14_(1))); for (double x = 0.01; x < 3; x += 0.01) { - EXPECT_THAT(approximation.Evaluate(x), + EXPECT_THAT(interpolant.Evaluate(x), AbsoluteErrorFrom(f(x), AllOf(Lt(7.2e-15), Ge(0)))); } } @@ -68,16 +68,16 @@ TEST(ApproximationTest, Exp) { TEST(ApproximationTest, AdaptiveSinInverse) { auto const f = [](double const x) { return Sin(1 * Radian / x); }; double error_estimate; - auto const approximations = + auto const interpolants = AdaptiveЧебышёвPolynomialInterpolant<8>(f, /*lower_bound=*/0.1, /*upper_bound=*/10.0, /*max_error=*/1e-6, &error_estimate); EXPECT_THAT(error_estimate, IsNear(7.1e-7_(1))); - EXPECT_THAT(approximations, SizeIs(11)); + EXPECT_THAT(interpolants, SizeIs(11)); EXPECT_THAT( - approximations, + interpolants, ElementsAre(AllOf(Property(&ЧебышёвSeries::lower_bound, AlmostEquals(0.1, 0)), Property(&ЧебышёвSeries::upper_bound, @@ -133,11 +133,11 @@ TEST(ApproximationTest, AdaptiveSinInverse) { Property(&ЧебышёвSeries::upper_bound, AlmostEquals(10.0, 0)), Property(&ЧебышёвSeries::degree, 8)))); - for (auto const& approximation : approximations) { - for (double x = approximation.lower_bound(); - x < approximation.upper_bound(); + for (auto const& interpolant : interpolants) { + for (double x = interpolant.lower_bound(); + x < interpolant.upper_bound(); x += 0.01) { - EXPECT_THAT(approximation.Evaluate(x), + EXPECT_THAT(interpolant.Evaluate(x), AbsoluteErrorFrom(f(x), AllOf(Lt(6.2e-7), Ge(2.7e-17)))); } } From 81d75887c923760c0c426bebb51c93bb446898ed Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 15:51:28 +0100 Subject: [PATCH 3/4] Typo. --- numerics/approximation.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/numerics/approximation.hpp b/numerics/approximation.hpp index 0a1ff07dc5..a3ccbcd7ab 100644 --- a/numerics/approximation.hpp +++ b/numerics/approximation.hpp @@ -29,7 +29,7 @@ template Difference> const& max_error, Difference>* error_estimate = nullptr); -// Returns an ordered vector of Чебышёв polynomial interpolants of f over, which +// Returns an ordered vector of Чебышёв polynomial interpolants of f, which // together cover [lower_bound, upper_bound]. Subdivides the interval so that // the degree of each approximant doesn't exceed |max_degree|. The final // (estimated) absolute error is guaranteed to be below |max_error|. From 6061d4023baa8adedbc30f5f8cc73f1df6bb2a33 Mon Sep 17 00:00:00 2001 From: pleroy Date: Sat, 20 Jan 2024 15:52:05 +0100 Subject: [PATCH 4/4] Lint. --- numerics/approximation_body.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/numerics/approximation_body.hpp b/numerics/approximation_body.hpp index 2a4540757a..6e9f73146c 100644 --- a/numerics/approximation_body.hpp +++ b/numerics/approximation_body.hpp @@ -3,6 +3,7 @@ #include "numerics/approximation.hpp" #include +#include #include #include "base/tags.hpp"