Skip to content

Commit

Permalink
Merge pull request #3850 from pleroy/Subdivision2
Browse files Browse the repository at this point in the history 
Adaptive Чебышёв polynomial interpolant with interval subdivision
  • Loading branch information
@pleroy
pleroy authored Jan 20, 2024
2 parents 2dcfd8f + 6061d40 commit 4eb331d
Show file tree
Hide file tree
Showing 3 changed files with 162 additions and 14 deletions.
17 changes: 16 additions & 1 deletion numerics/approximation.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
#pragma once

#include <type_traits>
#include <vector>

#include "numerics/чебышёв_series.hpp"
#include "quantities/named_quantities.hpp"
Expand All @@ -16,7 +17,7 @@ using namespace principia::quantities::_named_quantities;
template<typename Argument, typename Function>
using Value = std::invoke_result_t<Function, Argument>;

// 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.
Expand All @@ -28,8 +29,22 @@ template<int max_degree, typename Argument, typename Function>
Difference<Value<Argument, Function>> const& max_error,
Difference<Value<Argument, Function>>* error_estimate = nullptr);

// 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|.
template<int max_degree, typename Argument, typename Function>
std::vector<ЧебышёвSeries<Value<Argument, Function>, Argument>>
AdaptiveЧебышёвPolynomialInterpolant(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
Difference<Value<Argument, Function>> const& max_error,
Difference<Value<Argument, Function>>* error_estimate = nullptr);

} // namespace internal

using internal::AdaptiveЧебышёвPolynomialInterpolant;
using internal::ЧебышёвPolynomialInterpolant;

} // namespace _approximation
Expand Down
56 changes: 52 additions & 4 deletions numerics/approximation_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
#include "numerics/approximation.hpp"

#include <algorithm>
#include <utility>
#include <vector>

#include "base/tags.hpp"
Expand Down Expand Up @@ -126,13 +127,60 @@ template<int max_degree, typename Argument, typename Function>
FixedVector<Value<Argument, Function>, 2> const aⱼ(
{0.5 * (f_b + f_a), 0.5 * (f_b - f_a)});

return ЧебышёвPolynomialInterpolantImplementation</*N=*/2,
max_degree,
Argument,
Function>(
return ЧебышёвPolynomialInterpolantImplementation</*N=*/2, max_degree>(
f, a, b, max_error, fₖ, aⱼ, error_estimate);
}

template<int max_degree, typename Argument, typename Function>
std::vector<ЧебышёвSeries<Value<Argument, Function>, Argument>>
AdaptiveЧебышёвPolynomialInterpolant(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
Difference<Value<Argument, Function>> const& max_error,
Difference<Value<Argument, Function>>* error_estimate) {
// Try to build an interpolation over the entire interval.
Difference<Value<Argument, Function>> full_error_estimate;
auto full_interpolant = ЧебышёвPolynomialInterpolant<max_degree>(
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<Value<Argument, Function>, 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<Value<Argument, Function>> upper_error_estimate;
Difference<Value<Argument, Function>> lower_error_estimate;
auto const midpoint =
Barycentre(std::pair(lower_bound, upper_bound), std::pair(1, 1));
auto lower_interpolants =
AdaptiveЧебышёвPolynomialInterpolant<max_degree>(
f, lower_bound, midpoint, max_error, &lower_error_estimate);
auto upper_interpolants =
AdaptiveЧебышёвPolynomialInterpolant<max_degree>(
f, midpoint, upper_bound, max_error, &upper_error_estimate);
std::vector<ЧебышёвSeries<Value<Argument, Function>, 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
Expand Down
103 changes: 94 additions & 9 deletions numerics/approximation_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -13,51 +15,134 @@ 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); };
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))));
}
}

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))));
}
}

TEST(ApproximationTest, AdaptiveSinInverse) {
auto const f = [](double const x) { return Sin(1 * Radian / x); };
double error_estimate;
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(interpolants, SizeIs(11));
EXPECT_THAT(
interpolants,
ElementsAre(AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 512, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 512, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 256, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 256, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 128, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 128, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 64, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 64, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 32, 1)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 32, 1)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 * 3 / 64, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 * 3 / 64, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 16, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 16, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 8, 1)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 8, 1)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 4, 1)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 4, 1)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(0.1 + 9.9 / 2, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8)),
AllOf(Property(&ЧебышёвSeries<double, double>::lower_bound,
AlmostEquals(0.1 + 9.9 / 2, 0)),
Property(&ЧебышёвSeries<double, double>::upper_bound,
AlmostEquals(10.0, 0)),
Property(&ЧебышёвSeries<double, double>::degree, 8))));
for (auto const& interpolant : interpolants) {
for (double x = interpolant.lower_bound();
x < interpolant.upper_bound();
x += 0.01) {
EXPECT_THAT(interpolant.Evaluate(x),
AbsoluteErrorFrom(f(x), AllOf(Lt(6.2e-7), Ge(2.7e-17))));
}
}
}

} // namespace _approximation
} // namespace numerics
} // namespace principia

0 comments on commit 4eb331d

Please sign in to comment.