Skip to content

Commit

Permalink
Merge pull request #3849 from pleroy/Subdivision
Browse files Browse the repository at this point in the history 
Improvements to Чебышёв polynomial approximants
  • Loading branch information
@pleroy
pleroy authored Jan 20, 2024
2 parents 911d052 + 8c9ea2e commit 2dcfd8f
Show file tree
Hide file tree
Showing 4 changed files with 75 additions and 38 deletions.
9 changes: 7 additions & 2 deletions numerics/approximation.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,12 +16,17 @@ using namespace principia::quantities::_named_quantities;
template<typename Argument, typename Function>
using Value = std::invoke_result_t<Function, Argument>;

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

} // namespace internal

Expand Down
66 changes: 43 additions & 23 deletions numerics/approximation_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
#include <algorithm>
#include <vector>

#include "base/tags.hpp"
#include "geometry/barycentre_calculator.hpp"
#include "numerics/fixed_arrays.hpp"
#include "quantities/elementary_functions.hpp"
Expand All @@ -15,22 +16,39 @@ namespace numerics {
namespace _approximation {
namespace internal {

using namespace principia::base::_tags;
using namespace principia::geometry::_barycentre_calculator;
using namespace principia::numerics::_fixed_arrays;
using namespace principia::quantities::_elementary_functions;
using namespace principia::quantities::_si;

constexpr std::int64_t max_чебышёв_degree = 500;
// Compute the interpolation matrix and cache it in a static variable.
template<int N>
FixedMatrix<double, N + 1, N + 1> const& ЧебышёвInterpolationMatrix() {
static FixedMatrix<double, N + 1, N + 1> const ℐ = []() {
FixedMatrix<double, N + 1, N + 1> ℐ(uninitialized);
for (std::int64_t j = 0; j <= N; ++j) {
double const pⱼ = j == 0 || j == N ? 2 : 1;
for (std::int64_t k = 0; k <= N; ++k) {
double const pₖ = k == 0 || k == N ? 2 : 1;
ℐ(j, k) = 2 * Cos(π * j * k * Radian / N) / (pⱼ * pₖ * N);
}
}
return ℐ;
}();
return ℐ;
}

template<int N, typename Argument, typename Function>
template<int N, int max_degree, typename Argument, typename Function>
ЧебышёвSeries<Value<Argument, Function>, Argument>
ЧебышёвPolynomialInterpolantImplementation(
Function const& f,
Argument const& a,
Argument const& b,
Difference<Value<Argument, Function>> const& max_error,
FixedVector<Value<Argument, Function>, N / 2 + 1> const& previous_fₖ,
FixedVector<Value<Argument, Function>, N / 2 + 1> const& previous_aⱼ) {
FixedVector<Value<Argument, Function>, N / 2 + 1> const& previous_aⱼ,
Difference<Value<Argument, Function>>* const error_estimate) {
// This implementation follows [Boy13], section 4 and appendix A.
auto const midpoint = Barycentre(std::pair{a, b}, std::pair{0.5, 0.5});

Expand All @@ -39,7 +57,7 @@ template<int N, typename Argument, typename Function>
return 0.5 * (b - a) * Cos(π * k * Radian / N) + midpoint;
};

FixedVector<Value<Argument, Function>, N + 1> fₖ;
FixedVector<Value<Argument, Function>, N + 1> fₖ(uninitialized);

// Reuse the previous evaluations of |f|.
for (std::int64_t k = 0; k <= N / 2; ++k) {
Expand All @@ -51,32 +69,28 @@ template<int N, typename Argument, typename Function>
fₖ[k] = f(чебышёв_lobato_point(k));
}

FixedMatrix<double, N + 1, N + 1> ℐⱼₖ;
for (std::int64_t j = 0; j <= N; ++j) {
double const pⱼ = j == 0 || j == N ? 2 : 1;
for (std::int64_t k = 0; k <= N; ++k) {
double const pₖ = k == 0 || k == N ? 2 : 1;
ℐⱼₖ(j, k) = 2 * Cos(π * j * k * Radian / N) / (pⱼ * pₖ * N);
}
}

// Compute the coefficients of the Чебышёв polynomial.
FixedMatrix<double, N + 1, N + 1> const& ℐⱼₖ =
ЧебышёвInterpolationMatrix<N>();
auto const aⱼ = ℐⱼₖ * fₖ;

// Compute an upper bound for the error, based on the previous and new
// polynomials.
Difference<Value<Argument, Function>> error_estimate{};
Difference<Value<Argument, Function>> current_error_estimate{};
for (std::int64_t j = 0; j <= N / 2; ++j) {
error_estimate += Abs(previous_aⱼ[j] - aⱼ[j]);
current_error_estimate += Abs(previous_aⱼ[j] - aⱼ[j]);
}
for (std::int64_t j = N / 2 + 1; j <= N; ++j) {
error_estimate += Abs(aⱼ[j]);
current_error_estimate += Abs(aⱼ[j]);
}

if constexpr (2 * N <= max_чебышёв_degree) {
if (error_estimate > max_error) {
return ЧебышёвPolynomialInterpolantImplementation<2 * N>(
f, a, b, max_error, fₖ, aⱼ);
if constexpr (N <= max_degree) {
if (current_error_estimate > max_error) {
// Note that this recursive call overflows the stack when
// max_degree >= 256. We could allocate on the heap, but then we don't
// care about very high degree polynomials.
return ЧебышёвPolynomialInterpolantImplementation<2 * N, max_degree>(
f, a, b, max_error, fₖ, aⱼ, error_estimate);
}
}

Expand All @@ -86,31 +100,37 @@ template<int N, typename Argument, typename Function>
// 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|.
if (error_estimate != nullptr) {
*error_estimate = current_error_estimate;
}
std::vector<Value<Argument, Function>> coefficients;
std::copy(previous_aⱼ.begin(), previous_aⱼ.end(),
std::back_inserter(coefficients));
return ЧебышёвSeries<Value<Argument, Function>, Argument>(
coefficients, a, b);
}

template<typename Argument, typename Function>
template<int max_degree, typename Argument, typename Function>
ЧебышёвSeries<Value<Argument, Function>, Argument>
ЧебышёвPolynomialInterpolant(
Function const& f,
Argument const& lower_bound,
Argument const& upper_bound,
Difference<Value<Argument, Function>> const& max_error) {
Difference<Value<Argument, Function>> const& max_error,
Difference<Value<Argument, Function>>* const error_estimate) {
auto const& a = lower_bound;
auto const& b = upper_bound;
auto const f_a = f(a);
auto const f_b = f(b);
FixedVector<Value<Argument, Function>, 2> const fₖ({f_b, f_a});
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>(
f, a, b, max_error, fₖ, aⱼ);
f, a, b, max_error, fₖ, aⱼ, error_estimate);
}

} // namespace internal
Expand Down
27 changes: 19 additions & 8 deletions numerics/approximation_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include "quantities/elementary_functions.hpp"
#include "quantities/named_quantities.hpp"
#include "quantities/si.hpp"
#include "testing_utilities/approximate_quantity.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/numerics_matchers.hpp"

namespace principia {
Expand All @@ -14,19 +16,25 @@ namespace _approximation {
using namespace principia::quantities::_elementary_functions;
using namespace principia::quantities::_named_quantities;
using namespace principia::quantities::_si;
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 =
ЧебышёвPolynomialInterpolant<double>(f,
/*lower_bound=*/0.1,
/*upper_bound=*/10,
/*max_error=*/1e-6);
ЧебышёвPolynomialInterpolant<128>(f,
/*lower_bound=*/0.1,
/*upper_bound=*/10.0,
/*max_error=*/1e-6,
&error_estimate);
EXPECT_EQ(128, approximation.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),
AbsoluteErrorFrom(f(x), AllOf(Lt(3.0e-6), Ge(0))));
Expand All @@ -35,12 +43,15 @@ TEST(ApproximationTest, SinInverse) {

TEST(ApproximationTest, Exp) {
auto const f = [](double const x) { return std::exp(x); };
double error_estimate;
auto const approximation =
ЧебышёвPolynomialInterpolant<double>(f,
/*lower_bound=*/0.01,
/*upper_bound=*/3,
/*max_error=*/1e-6);
ЧебышёвPolynomialInterpolant<128>(f,
/*lower_bound=*/0.01,
/*upper_bound=*/3.0,
/*max_error=*/1e-6,
&error_estimate);
EXPECT_EQ(16, approximation.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),
AbsoluteErrorFrom(f(x), AllOf(Lt(7.2e-15), Ge(0))));
Expand Down
11 changes: 6 additions & 5 deletions physics/apsides_body.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -322,11 +322,12 @@ ComputeFirstCollision(
};

// Interpolate the height above the terrain using a Чебышёв polynomial.
auto const чебышёв_interpolant = ЧебышёвPolynomialInterpolant(
height_above_terrain_at_time,
interval.min,
interval.max,
reference_body.mean_radius() * max_error_relative_to_radius);
auto const чебышёв_interpolant =
ЧебышёвPolynomialInterpolant</*max_degree=*/128>(
height_above_terrain_at_time,
interval.min,
interval.max,
reference_body.mean_radius() * max_error_relative_to_radius);
auto const& real_roots = чебышёв_interpolant.RealRoots(
max_error_relative_to_radius);
if (real_roots.empty()) {
Expand Down

0 comments on commit 2dcfd8f

Please sign in to comment.