Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improvements to Чебышёв polynomial approximants #3849

Merged
merged 3 commits into from
Jan 20, 2024
Merged

Improvements to Чебышёв polynomial approximants #3849

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
78e5ea3
Explicit control of max_degree.
pleroy Jan 20, 2024
15a75d5
Return the error estimate.
pleroy Jan 20, 2024
8c9ea2e
IWYUsing.
pleroy Jan 20, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 7 additions & 2 deletions numerics/approximation.hpp
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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