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

Commits on Sep 15, 2020

  1. Copy the full SHA
    c5c3fe8 View commit details

Commits on Sep 17, 2020

  1. Copy the full SHA
    7400b36 View commit details
  2. Use the generated tables.

    pleroy committed Sep 17, 2020
    Copy the full SHA
    1a3284e View commit details
  3. More tests.

    pleroy committed Sep 17, 2020
    Copy the full SHA
    4ceaf05 View commit details
  4. Cleanup.

    pleroy committed Sep 17, 2020
    Copy the full SHA
    bf7f796 View commit details

Commits on Sep 18, 2020

  1. Merge pull request #2721 from pleroy/Gauss

    Gauss-Legendre quadrature
    pleroy authored Sep 18, 2020
    Copy the full SHA
    61f8ccb View commit details
Showing with 142 additions and 0 deletions.
  1. +1 −0 numerics/numerics.vcxproj
  2. +3 −0 numerics/numerics.vcxproj.filters
  3. +7 −0 numerics/quadrature.hpp
  4. +31 −0 numerics/quadrature_body.hpp
  5. +100 −0 numerics/quadrature_test.cpp
1 change: 1 addition & 0 deletions numerics/numerics.vcxproj
Original file line number Diff line number Diff line change
@@ -94,6 +94,7 @@
<ClCompile Include="poisson_series_test.cpp" />
<ClCompile Include="polynomial_evaluators_test.cpp" />
<ClCompile Include="polynomial_test.cpp" />
<ClCompile Include="quadrature_test.cpp" />
<ClCompile Include="root_finders_test.cpp" />
<ClCompile Include="unbounded_arrays_test.cpp" />
<ClCompile Include="чебышёв_series_test.cpp" />
3 changes: 3 additions & 0 deletions numerics/numerics.vcxproj.filters
Original file line number Diff line number Diff line change
@@ -258,6 +258,9 @@
<ClCompile Include="poisson_series_basis_test.cpp">
<Filter>Test Files</Filter>
</ClCompile>
<ClCompile Include="quadrature_test.cpp">
<Filter>Source Files</Filter>
</ClCompile>
</ItemGroup>
<ItemGroup>
<Text Include="xgscd.proto.txt">
7 changes: 7 additions & 0 deletions numerics/quadrature.hpp
Original file line number Diff line number Diff line change
@@ -12,6 +12,12 @@ namespace internal_quadrature {

using quantities::Primitive;

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> GaussLegendre(
Function const& function,
Argument const& lower_bound,
Argument const& upper_bound);

template<typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(
Function const& function,
@@ -21,6 +27,7 @@ Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(

} // namespace internal_quadrature

using internal_quadrature::GaussLegendre;
using internal_quadrature::Midpoint;

} // namespace quadrature
31 changes: 31 additions & 0 deletions numerics/quadrature_body.hpp
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

#pragma once

#include "numerics/gauss_legendre_weights.mathematica.h"
#include "numerics/legendre_roots.mathematica.h"
#include "numerics/quadrature.hpp"

namespace principia {
@@ -10,6 +12,35 @@ namespace internal_quadrature {

using quantities::Difference;

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> Gauss(
Function const& function,
Argument const& lower_bound,
Argument const& upper_bound,
double const* const nodes,
double const* const weights) {
Difference<Argument> half_width = (upper_bound - lower_bound) / 2;
std::invoke_result_t<Function, Argument> result{};
for (int i = 0; i < points; ++i) {
Argument const scaled_node = lower_bound + half_width * (nodes[i] + 1);
// TODO(phl): Consider compensated summation.
result += weights[i] * function(scaled_node);
}
return result * half_width;
}

template<int points, typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> GaussLegendre(
Function const& function,
Argument const& lower_bound,
Argument const& upper_bound) {
return Gauss<points>(function,
lower_bound,
upper_bound,
LegendreRoots[points],
GaussLegendreWeights[points]);
}

template<typename Argument, typename Function>
Primitive<std::invoke_result_t<Function, Argument>, Argument> Midpoint(
Function const& function,
100 changes: 100 additions & 0 deletions numerics/quadrature_test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,100 @@

#include "numerics/quadrature.hpp"

#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "quantities/elementary_functions.hpp"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/numerics_matchers.hpp"

namespace principia {
namespace numerics {
namespace quadrature {

using quantities::Angle;
using quantities::Cos;
using quantities::Sin;
using quantities::si::Radian;
using testing_utilities::AlmostEquals;
using testing_utilities::IsNear;
using testing_utilities::RelativeErrorFrom;
using testing_utilities::operator""_⑴;

class QuadratureTest : public ::testing::Test {};

TEST_F(QuadratureTest, Sin) {
auto const f = [](Angle const x) { return Sin(x); };
auto const ʃf = (Cos(2.0 * Radian) - Cos(5.0 * Radian)) * Radian;
EXPECT_THAT(GaussLegendre<1>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(10_⑴)));
EXPECT_THAT(GaussLegendre<2>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(3.3_⑴)));
EXPECT_THAT(GaussLegendre<3>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(4.0e-1_⑴)));
EXPECT_THAT(GaussLegendre<4>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(2.4e-2_⑴)));
EXPECT_THAT(GaussLegendre<5>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(8.6e-4_⑴)));
EXPECT_THAT(GaussLegendre<6>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(2.1e-5_⑴)));
EXPECT_THAT(GaussLegendre<7>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(3.6e-7_⑴)));
EXPECT_THAT(GaussLegendre<8>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(4.7e-9_⑴)));
EXPECT_THAT(GaussLegendre<9>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(4.8e-11_⑴)));
EXPECT_THAT(GaussLegendre<10>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 2498));
EXPECT_THAT(GaussLegendre<11>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 20));
EXPECT_THAT(GaussLegendre<12>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 6));
EXPECT_THAT(GaussLegendre<13>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 1));
EXPECT_THAT(GaussLegendre<14>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 1));
EXPECT_THAT(GaussLegendre<15>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 4));
}

TEST_F(QuadratureTest, Sin10) {
auto const f = [](Angle const x) { return Sin(2 * x); };
auto const ʃf = (Cos(4 * Radian) - Cos(10 * Radian)) / 2 * Radian;
EXPECT_THAT(GaussLegendre<1>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(9.7_⑴)));
EXPECT_THAT(GaussLegendre<2>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(7.6_⑴)));
EXPECT_THAT(GaussLegendre<3>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(7.6_⑴)));
EXPECT_THAT(GaussLegendre<4>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(2.4_⑴)));
EXPECT_THAT(GaussLegendre<5>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(4.2e-1_⑴)));
EXPECT_THAT(GaussLegendre<6>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(4.6e-2_⑴)));
EXPECT_THAT(GaussLegendre<7>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(3.5e-3_⑴)));
EXPECT_THAT(GaussLegendre<8>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(2.0e-4_⑴)));
EXPECT_THAT(GaussLegendre<9>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(8.5e-6_⑴)));
EXPECT_THAT(GaussLegendre<10>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(2.9e-7_⑴)));;
EXPECT_THAT(GaussLegendre<11>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(8.1e-9_⑴)));
EXPECT_THAT(GaussLegendre<12>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(1.9e-10_⑴)));
EXPECT_THAT(GaussLegendre<13>(f, -2.0 * Radian, 5.0 * Radian),
RelativeErrorFrom(ʃf, IsNear(3.7e-12_⑴)));
EXPECT_THAT(GaussLegendre<14>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 422));
EXPECT_THAT(GaussLegendre<15>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 44));
EXPECT_THAT(GaussLegendre<16>(f, -2.0 * Radian, 5.0 * Radian),
AlmostEquals(ʃf, 152));
}

} // namespace quadrature
} // namespace numerics
} // namespace principia