Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 734a80d024c7
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: a69de59f1bca
Choose a head ref
  • 9 commits
  • 4 files changed
  • 1 contributor

Commits on Oct 27, 2020

  1. Copy the full SHA
    2e3384f View commit details
  2. Bibliography.

    pleroy committed Oct 27, 2020
    Copy the full SHA
    eee11a5 View commit details
  3. Cleanup.

    pleroy committed Oct 27, 2020
    Copy the full SHA
    a3e55fb View commit details
  4. Use double precision.

    pleroy committed Oct 27, 2020
    Copy the full SHA
    e4bdc77 View commit details
  5. Some cleanup.

    pleroy committed Oct 27, 2020
    Copy the full SHA
    75c7c55 View commit details
  6. Tighter error bounds.

    pleroy committed Oct 27, 2020
    Copy the full SHA
    dbac23a View commit details

Commits on Oct 28, 2020

  1. Copy the full SHA
    8bc2d8c View commit details
  2. Fix tolerances.

    pleroy committed Oct 28, 2020
    Copy the full SHA
    0f35444 View commit details
  3. Merge pull request #2775 from pleroy/IntegrationByParts

    Speed-up PoissonSeries::Integrate
    pleroy authored Oct 28, 2020

    Verified

    This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
    Copy the full SHA
    a69de59 View commit details
Showing with 76 additions and 32 deletions.
  1. +26 −0 documentation/bibliography.bib
  2. BIN documentation/bibliography.pdf
  3. +5 −5 numerics/frequency_analysis_test.cpp
  4. +45 −27 numerics/poisson_series_body.hpp
26 changes: 26 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
@@ -881,6 +881,32 @@ @book{Warren2003
title = {Hacker's Delight},
}

@inbook{HuybrechsOlver2009,
author = {Huybrechs, D. and Olver, S.},
editor = {Engquist, B. and Fokas, A. and Hairer, E.},
language = {English},
publisher = {Cambridge University Press},
booktitle = {Highly Oscillatory Problems},
date = {2009},
doi = {10.1017/CBO9781139107136},
isbn = {9781139107136},
pages = {25--50},
title = {Highly oscillatory quadrature},
}

@inbook{IserlesNørsettOlver2006,
author = {Iserles, A. and Nørsett, S.P. and Olver, S.},
editor = {de Castro, A.B. and Gómez, D. and Quintela, P. and Salgado, P.},
language = {English},
publisher = {Springer Berlin Heidelberg},
booktitle = {Numerical Mathematics and Advanced Applications},
date = {2006},
doi = {10.1007/978-3-540-34288-5_6},
isbn = {978-3-540-34287-8},
pages = {97--118},
title = {Highly Oscillatory Quadrature: The Story so Far},
}

@incollection{BlanesCasasRos2001a,
author = {Blanes, S. and Casas, F. and Ros, J.},
editor = {Vulkov, Lubin and Yalamov, Plamen and Waśniewski, Jerzy},
Binary file modified documentation/bibliography.pdf
Binary file not shown.
10 changes: 5 additions & 5 deletions numerics/frequency_analysis_test.cpp
Original file line number Diff line number Diff line change
@@ -431,7 +431,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
? AllOf(Gt(6.7e-2 * Metre), Lt(7.9 * Metre))
: ω_index == 2
? AllOf(Gt(1.1e-4 * Metre), Lt(9.7e-1 * Metre))
: AllOf(Gt(2.1e-17 * Metre), Lt(6.6e-13 * Metre)))
: AllOf(Gt(1.4e-16 * Metre), Lt(2.8e-13 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
@@ -453,7 +453,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionNoSecular) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series.value()(t_min + i * (t_max - t_min) / 100),
AllOf(Ge(0), Lt(2.0e-11))));
AllOf(Ge(0), Lt(1.4e-11))));
}
}

@@ -496,8 +496,8 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
? AllOf(Gt(3.3e-2 * Metre), Lt(3.6 * Metre))
: ω_index == 3
? AllOf(Gt(7.5e-3 * Metre), Lt(5.4 * Metre))
: AllOf(Gt(1.5e-16 * Metre),
Lt(9.7e-14 * Metre)))
: AllOf(Gt(1.7e-16 * Metre),
Lt(8.7e-14 * Metre)))
<< ω_index;
}
if (ω_index == ωs.size()) {
@@ -519,7 +519,7 @@ TEST_F(FrequencyAnalysisTest, PoissonSeriesIncrementalProjectionSecular) {
EXPECT_THAT(
projection4(t_min + i * (t_max - t_min) / 100),
RelativeErrorFrom(series(t_min + i * (t_max - t_min) / 100),
AllOf(Ge(0), Lt(5.9e-15))));
AllOf(Ge(0), Lt(5.4e-15))));
}
}

72 changes: 45 additions & 27 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -49,6 +49,39 @@ constexpr int clenshaw_curtis_point_per_period = 4;
// two successive computations with increasing number of points.
constexpr double clenshaw_curtis_relative_error = 0x1p-32;

// This function computes ∫ₜ₁ᵗ²(p(t) sin ω t + q(t) cos ω t) dt.
template<typename Value,
int degree,
template<typename, typename, int> class Evaluator>
quantities::Primitive<Value, Time> AngularFrequencyIntegrate(
AngularFrequency const& ω,
PolynomialInMonomialBasis<Value, Instant, degree, Evaluator> const& p,
PolynomialInMonomialBasis<Value, Instant, degree, Evaluator> const& q,
Instant const& t1,
Instant const& t2,
double const sin_ωt1,
double const cos_ωt1,
double const sin_ωt2,
double const cos_ωt2) {
static_assert(degree >= 0);
DoublePrecision<Value> sum;
sum += q(t2) * sin_ωt2;
sum -= p(t2) * cos_ωt2;
sum -= q(t1) * sin_ωt1;
sum += p(t1) * cos_ωt1;
if constexpr (degree > 0) {
sum += AngularFrequencyIntegrate(ω,
/*p=*/-q.template Derivative<1>(),
/*q=*/p.template Derivative<1>(),
t1, t2,
sin_ωt1, cos_ωt1,
sin_ωt2, cos_ωt2);
}
return (sum.value + sum.error) / ω * Radian;
}

// This function computes ∫(p(t) sin ω t + q(t) cos ω t) dt where p and q are
// the two parts of the polynomials argument.
template<typename Value,
int aperiodic_degree, int periodic_degree,
template<typename, typename, int> class Evaluator>
@@ -293,37 +326,22 @@ quantities::Primitive<Value, Time>
PoissonSeries<Value, aperiodic_degree_, periodic_degree_, Evaluator>::
Integrate(Instant const& t1,
Instant const& t2) const {
using FirstPart = PoissonSeries<Value,
0, periodic_degree_,
Evaluator>;
using SecondPart = PoissonSeries<Variation<Value>,
0, periodic_degree_ - 1,
Evaluator>;
auto const aperiodic_primitive = aperiodic_.Primitive();
quantities::Primitive<Value, Time> result =
aperiodic_primitive(t2) - aperiodic_primitive(t1);
for (auto const& [ω, polynomials] : periodic_) {
// Integration by parts.
FirstPart const first_part(
typename FirstPart::AperiodicPolynomial({}, origin_),
{{ω,
{/*sin=*/typename FirstPart::PeriodicPolynomial(polynomials.cos),
/*cos=*/typename FirstPart::PeriodicPolynomial(-polynomials.sin)}}});
DoublePrecision<Value> sum;
sum += first_part(t2);
sum -= first_part(t1);

if constexpr (periodic_degree_ != 0) {
auto const sin_polynomial = -polynomials.cos.template Derivative<1>();
auto const cos_polynomial = polynomials.sin.template Derivative<1>();
SecondPart const second_part(
typename SecondPart::AperiodicPolynomial({}, origin_),
{{ω,
{/*sin=*/sin_polynomial,
/*cos=*/cos_polynomial}}});
sum += second_part.Integrate(t1, t2);
}
result += (sum.value + sum.error) / ω * Radian;
// This implementation follows [HO09], Theorem 1 and [INO06] equation 4.
// The trigonometric functions are computed only once as we iterate through
// the degree of the polynomials.
auto const sin_ωt1 = Sin(ω * (t1 - origin_));
auto const cos_ωt1 = Cos(ω * (t1 - origin_));
auto const sin_ωt2 = Sin(ω * (t2 - origin_));
auto const cos_ωt2 = Cos(ω * (t2 - origin_));
result += AngularFrequencyIntegrate(ω,
polynomials.sin, polynomials.cos,
t1, t2,
sin_ωt1, cos_ωt1,
sin_ωt2, cos_ωt2);
}
return result;
}