Skip to content

Commit

Permalink
Showing 3 changed files with 35 additions and 15 deletions.
17 changes: 9 additions & 8 deletions numerics/poisson_series.hpp
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@


#pragma once

#include <map>
@@ -16,6 +16,11 @@ using quantities::Product;
using quantities::Quotient;
using quantities::Time;

// A Poisson series is the sum of terms of the form:
// aₙtⁿ aₙₖ tⁿ sin ωₖ t aₙₖ tⁿ cos ωₖ t
// Terms of the first kind are called aperiodic, terms of the second kind are
// called periodic. Poisson series form an algebra that is stable by derivation
// and integration.
template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
class PoissonSeries {
@@ -38,36 +43,32 @@ class PoissonSeries {
Value Evaluate(Time const& t) const;

private:
Polynomial const aperiodic_;
PolynomialsByAngularFrequency const periodic_;
Polynomial aperiodic_;
// All the keys in this map are positive.
PolynomialsByAngularFrequency periodic_;

template<typename V, int r, template<typename, typename, int> class E>
PoissonSeries<V, r, E> friend operator-(PoissonSeries<V, r, E> const& right);

template<typename V, int l, int r, template<typename, typename, int> class E>
PoissonSeries<V, std::max(l, r), E> friend operator+(
PoissonSeries<V, l, E> const& left,
PoissonSeries<V, r, E> const& right);

template<typename V, int l, int r, template<typename, typename, int> class E>
PoissonSeries<V, std::max(l, r), E> friend operator-(
PoissonSeries<V, l, E> const& left,
PoissonSeries<V, r, E> const& right);

template<typename Scalar,
typename V, int degree_,
template<typename, typename, int> class E>
PoissonSeries<Product<Scalar, V>, degree_, E> friend operator*(
Scalar const& left,
PoissonSeries<V, degree_, E> const& right);

template<typename Scalar,
typename V, int degree_,
template<typename, typename, int> class E>
PoissonSeries<Product<V, Scalar>, degree_, E> friend operator*(
PoissonSeries<V, degree_, E> const& left,
Scalar const& right);

template<typename Scalar,
typename V, int degree_,
template<typename, typename, int> class E>
25 changes: 22 additions & 3 deletions numerics/poisson_series_body.hpp
Original file line number Diff line number Diff line change
@@ -15,14 +15,33 @@ using quantities::Cos;
using quantities::Infinity;
using quantities::Sin;

// TODO(phl): This file should probably take advantage of emplace_hint.

template<typename Value, int degree_,
template<typename, typename, int> class Evaluator>
PoissonSeries<Value, degree_, Evaluator>::PoissonSeries(
Polynomial const& aperiodic,
PolynomialsByAngularFrequency const& periodic)
: aperiodic_(aperiodic), periodic_(periodic) {
for (auto const& [ω, _] : periodic_) {
CHECK_LT(AngularFrequency{}, ω);
: aperiodic_(aperiodic) {
// The |periodic| map may have elements with positive or negative angular
// frequencies. Normalize our member variable to only have positive angular
// frequencies.
for (auto it = periodic.crbegin(); it != periodic.crend(); ++it) {
auto const ω = it->first;
if (ω < AngularFrequency{}) {
auto const positive_it = periodic_.find(-ω);
if (positive_it == periodic_.cend()) {
periodic_.emplace(-ω,
Polynomials{/*sin=*/-it->second.sin,
/*cos=*/it->second.cos});
} else {
positive_it->second = {
/*sin=*/positive_it->second.sin - it->second.sin,
/*cos=*/positive_it->second.cos + it->second.cos};
}
} else {
periodic_.insert(*it);
}
}
}

8 changes: 4 additions & 4 deletions numerics/poisson_series_test.cpp
Original file line number Diff line number Diff line change
@@ -35,7 +35,7 @@ TEST_F(PoissonSeriesTest, VectorSpace) {
Degree1::Polynomial psa2({13, 14 / Second});
Degree1::Polynomial pca2({15, 16 / Second});

Degree1::Polynomial psb3({17, 18 / Second});
Degree1::Polynomial psb3({-17, -18 / Second});
Degree1::Polynomial pcb3({19, 20 / Second});

Degree1::Polynomials pa1{/*sin=*/psa1, /*cos=*/pca1};
@@ -47,15 +47,15 @@ TEST_F(PoissonSeriesTest, VectorSpace) {

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

Degree1 pa(pa0, {{ω1, pa1}, {ω2, pa2}});
Degree1 pb(pb0, {{ω1, pb1}, {ω3, pb3}});

EXPECT_THAT(pa.Evaluate(1 * Second),
AlmostEquals(3 + 11 * Sin(1 * Radian) + 15 * Cos(1 * Radian) +
27 * Sin(2 * Radian) + 31 * Cos(2 * Radian),
0));
0, 1));
EXPECT_THAT(pb.Evaluate(1 * Second),
AlmostEquals(7 + 19 * Sin(1 * Radian) + 23 * Cos(1 * Radian) +
35 * Sin(3 * Radian) + 39 * Cos(3 * Radian),
@@ -96,7 +96,7 @@ TEST_F(PoissonSeriesTest, VectorSpace) {
{
auto const quotient = pb / 1.5;
EXPECT_THAT(quotient.Evaluate(1 * Second),
AlmostEquals(pb.Evaluate(1 * Second) / 1.5, 32));
AlmostEquals(pb.Evaluate(1 * Second) / 1.5, 0, 32));
}
}

0 comments on commit efa11e2

Please sign in to comment.