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: 57a213fca36b
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: f651d9e37d71
Choose a head ref
  • 11 commits
  • 3 files changed
  • 2 contributors

Commits on Jul 1, 2020

  1. Copy the full SHA
    2dfb881 View commit details
  2. Copy the full SHA
    38ad9e5 View commit details
  3. Polynomial conversion.

    pleroy committed Jul 1, 2020
    Copy the full SHA
    ecdf7b9 View commit details
  4. Something that compiles.

    pleroy committed Jul 1, 2020
    Copy the full SHA
    8893e53 View commit details

Commits on Jul 2, 2020

  1. Move code out.

    pleroy committed Jul 2, 2020
    Copy the full SHA
    7f23d30 View commit details
  2. Copy the full SHA
    8eb9d71 View commit details
  3. Test the primitive.

    pleroy committed Jul 2, 2020
    Copy the full SHA
    c3b9de2 View commit details
  4. Copy the full SHA
    0c2ef8d View commit details
  5. Cleanup.

    pleroy committed Jul 2, 2020
    Copy the full SHA
    8bb6420 View commit details
  6. More cleanup.

    pleroy committed Jul 2, 2020
    Copy the full SHA
    55fd8aa View commit details
  7. Merge pull request #2623 from pleroy/PoissonPrimitive

    Primitive of a Poisson series
    pleroy authored Jul 2, 2020
    Copy the full SHA
    f651d9e View commit details
Showing with 106 additions and 12 deletions.
  1. +4 −0 numerics/poisson_series.hpp
  2. +49 −0 numerics/poisson_series_body.hpp
  3. +53 −12 numerics/poisson_series_test.cpp
4 changes: 4 additions & 0 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -43,6 +43,10 @@ class PoissonSeries {

Value Evaluate(Time const& t) const;

// The constant term of the result is zero.
PoissonSeries<quantities::Primitive<Value, Time>, degree_ + 1, Evaluator>
Primitive() const;

private:
Polynomial aperiodic_;
// All the keys in this map are positive.
49 changes: 49 additions & 0 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -16,8 +16,40 @@ namespace internal_poisson_series {

using quantities::Cos;
using quantities::Infinity;
using quantities::Primitive;
using quantities::Sin;

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
typename PoissonSeries<Primitive<Value, Time>,
degree_ + 1,
Evaluator>::Polynomials
AngularFrequencyPrimitive(
AngularFrequency const& ω,
typename PoissonSeries<Value, degree_, Evaluator>::Polynomials const&
polynomials) {
using Argument = PoissonSeries<Value, degree_, Evaluator>;
using Result = PoissonSeries<Primitive<Value, Time>, degree_ + 1, Evaluator>;

// Integration by parts.
typename Result::Polynomials const first_part{
/*sin=*/Result::Polynomial(polynomials.cos / ω * Radian),
/*cos=*/Result::Polynomial(-polynomials.sin / ω * Radian)};
if constexpr (degree_ == 0) {
return first_part;
} else {
auto const sin_polynomial = -polynomials.cos.Derivative<1>() / ω * Radian;
auto const cos_polynomial = polynomials.sin.Derivative<1>() / ω * Radian;
auto const second_part =
AngularFrequencyPrimitive<Value, degree_ - 1, Evaluator>(
ω,
{/*sin=*/sin_polynomial,
/*cos=*/cos_polynomial});
return {/*sin=*/first_part.sin + second_part.sin,
/*cos=*/first_part.cos + second_part.cos};
}
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
@@ -58,6 +90,23 @@ Value PoissonSeries<Value, degree_, Evaluator>::Evaluate(Time const& t) const {
return result;
}

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<quantities::Primitive<Value, Time>, degree_ + 1, Evaluator>
PoissonSeries<Value, degree_, Evaluator>::Primitive() const {
using Result =
PoissonSeries<quantities::Primitive<Value, Time>, degree_ + 1, Evaluator>;
typename Result::PolynomialsByAngularFrequency periodic;
typename Result::Polynomial const aperiodic = aperiodic_.Primitive();
for (auto const& [ω, polynomials] : periodic_) {
periodic.emplace_hint(
periodic.cend(),
ω,
AngularFrequencyPrimitive<Value, degree_, Evaluator>(ω, polynomials));
}
return Result{aperiodic, std::move(periodic)};
}

template<typename Value, int rdegree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, rdegree_, Evaluator> operator+(
65 changes: 53 additions & 12 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -16,6 +16,7 @@ namespace numerics {
using quantities::AngularFrequency;
using quantities::Cos;
using quantities::Sin;
using quantities::Time;
using quantities::si::Radian;
using quantities::si::Second;
using testing_utilities::AlmostEquals;
@@ -24,7 +25,11 @@ class PoissonSeriesTest : public ::testing::Test {
protected:
using Degree1 = PoissonSeries<double, 1, HornerEvaluator>;

PoissonSeriesTest() {
PoissonSeriesTest()
: ω0_(0 * Radian / Second),
ω1_(1 * Radian / Second),
ω2_(2 * Radian / Second),
ω3_(-3 * Radian / Second) {
Degree1::Polynomial pa0({0, 0 / Second});
Degree1::Polynomial psa0({100, 200 / Second});
Degree1::Polynomial pca0({1, 2 / Second});
@@ -50,19 +55,19 @@ class PoissonSeriesTest : public ::testing::Test {

Degree1::Polynomials pscb3{/*sin=*/psb3, /*cos=*/pcb3};

AngularFrequency ω0 = 0 * Radian / Second;
AngularFrequency ω1 = 1 * Radian / Second;
AngularFrequency ω2 = 2 * Radian / Second;
AngularFrequency ω3 = -3 * Radian / Second;

pa_ = std::make_unique<Degree1>(pa0,
Degree1::PolynomialsByAngularFrequency{
0, psca0}, {ω1, psca1}, {ω2, psca2}});
pb_ = std::make_unique<Degree1>(pb0,
Degree1::PolynomialsByAngularFrequency{
1, pscb1}, {ω3, pscb3}});
pa_ = std::make_unique<Degree1>(
pa0,
Degree1::PolynomialsByAngularFrequency{
{ω0_, psca0}, {ω1_, psca1}, {ω2_, psca2}});
pb_ = std::make_unique<Degree1>(
pb0,
Degree1::PolynomialsByAngularFrequency{{ω1_, pscb1}, {ω3_, pscb3}});
}

AngularFrequency ω0_;
AngularFrequency ω1_;
AngularFrequency ω2_;
AngularFrequency ω3_;
std::unique_ptr<Degree1> pa_;
std::unique_ptr<Degree1> pb_;
};
@@ -125,5 +130,41 @@ TEST_F(PoissonSeriesTest, Algebra) {
AlmostEquals(pa_->Evaluate(1 * Second) * pb_->Evaluate(1 * Second), 6));
}

TEST_F(PoissonSeriesTest, Primitive) {
auto const actual_primitive = pb_->Primitive();

// The primitive was computed using Mathematica.
auto const expected_primitive = [=](Time const& t){
auto const a0 = 3;
auto const a1 = 4 / Second;
auto const b0 = 9;
auto const b1 = 10 / Second;
auto const c0 = 11;
auto const c1 = 12 / Second;
auto const d0 = -17;
auto const d1 = -18 / Second;
auto const e0 = 19;
auto const e1 = 20 / Second;
return a0 * t + (a1 * t * t) / 2 +
(c1 * Cos(ω1_ * t) * Radian * Radian) / (ω1_ * ω1_) -
(b0 * Cos(ω1_ * t) * Radian) / ω1_ -
(b1 * t * Cos(ω1_ * t) * Radian) / ω1_ +
(e1 * Cos(ω3_ * t) * Radian * Radian) / (ω3_ * ω3_) -
(d0 * Cos(ω3_ * t) * Radian) / ω3_ -
(d1 * t * Cos(ω3_ * t) * Radian) / ω3_ +
(b1 * Sin(ω1_ * t) * Radian * Radian) / (ω1_ * ω1_) +
(c0 * Sin(ω1_ * t) * Radian) / ω1_ +
(c1 * t * Sin(ω1_ * t) * Radian) / ω1_ +
(d1 * Sin(ω3_ * t) * Radian * Radian) / (ω3_ * ω3_) +
(e0 * Sin(ω3_ * t) * Radian) / ω3_ +
(e1 * t * Sin(ω3_ * t) * Radian) / ω3_;
};

for (int i = -10; i < 10; ++i) {
EXPECT_THAT(actual_primitive.Evaluate(i * Second),
AlmostEquals(expected_primitive(i * Second), 0, 6));
}
}

} // namespace numerics
} // namespace principia