diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index ed8007814b..0b91f26b86 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -120,9 +120,12 @@ jobs: shell: pwsh run: | &.\include_what_you_using_all_the_things.ps1 - $files = (git diff --name-only) - foreach ($file in $files) { - $diff = (git diff $file) - echo "::warning file=$file,line=1,title=Include what you using::include_what_you_using_all_the_things.ps1 modifies this file as follows:%0A$([string]::join('%0A', (git diff $file)))" + $files = (git diff --name-only -z) + if ($files) { + $files = $files.Split("`0", [StringSplitOptions]::RemoveEmptyEntries) + foreach ($file in $files) { + $diff = (git diff $file) + echo "::warning file=$file,line=1,title=Include what you using::include_what_you_using_all_the_things.ps1 modifies this file as follows:%0A$([string]::join('%0A', (git diff $file)))" + } + exit $files.Length } - exit $files.Length diff --git a/documentation/bibliography.bib b/documentation/bibliography.bib index fd7a231e95..3dfb1d6b2e 100644 --- a/documentation/bibliography.bib +++ b/documentation/bibliography.bib @@ -44,6 +44,28 @@ @article{AriasCharlotFeisselLestrade1995 volume = {303}, } +@article{BaiJunkins2010, + author = {Bai, Xiaoli and Junkins, John L}, + date = {2010}, + journaltitle = {Advances in the Astronautical Sciences}, + number = {2}, + pages = {1459--1476}, + title = {Solving initial value problems by the Picard-Chebyshev method with NVIDIA GPUs}, + volume = {136}, +} + +@article{BaiJunkins2012, + author = {Bai, Xiaoli and Junkins, John L}, + publisher = {Springer}, + date = {2012}, + doi = {10.1007/s40295-013-0021-6}, + journaltitle = {The Journal of the Astronautical Sciences}, + number = {1}, + pages = {327--351}, + title = {Modified Chebyshev-Picard iteration methods for solution of initial value problems}, + volume = {59}, +} + @article{Bateman1938, author = {Bateman, Harry}, date = {1938-01}, @@ -393,6 +415,18 @@ @article{ChinKidwell2000 volume = {62}, } +@article{ClenshawNorton1963, + author = {Clenshaw, C. W. and Norton, H. J.}, + publisher = {Oxford University Press}, + date = {1963}, + doi = {10.1093/comjnl/6.1.88}, + journaltitle = {The Computer Journal}, + number = {1}, + pages = {88--92}, + title = {The solution of nonlinear ordinary differential equations in Chebyshev series}, + volume = {6}, +} + @article{Cohen1951, author = {Cohen, Jr., A. C.}, date = {1951-06}, @@ -581,6 +615,18 @@ @article{FarrésLaskarBlanesCasasMakazagaMurua2013 volume = {116}, } +@article{FeaginNacozy1983, + author = {Feagin, Terry and Nacozy, Paul}, + publisher = {Springer}, + date = {1983}, + doi = {10.1007/BF01232802}, + journaltitle = {Celestial mechanics}, + number = {2}, + pages = {107--115}, + title = {Matrix formulation of the Picard method for parallel computation}, + volume = {29}, +} + @article{FiengaMacheLaskarGastineau2008, author = {Fienga, A. and Manche, H. and Laskar, J. and Gastineau, M.}, url = {http://dx.doi.org/10.1051/0004-6361:20066607}, @@ -654,6 +700,13 @@ @article{Fornberg1988 volume = {51}, } +@article{Fukushima1997, + author = {Fukushima, Toshio}, + doi = {10.1086/118404}, + journaltitle = {Astronomical Journal v. 113, p. 1909-1914 (1997)}, + title = {Picard iteration method, Chebyshev polynomial approximation, and global numerical integration of dynamical motions}, +} + @article{Fukushima2009a, author = {Fukushima, Toshio}, date = {2009-10-25}, @@ -2242,6 +2295,14 @@ @report{BlancSchgounn1996 title = {AVISO/Altimetry, AVISO User Handbook for Merged TOPEX/\-POSEIDON products, Edition 3.0}, } +@report{ByrdLee1969, + author = {Byrd, P. F. and Lee, K. L.}, + date = {1969}, + number = {19690027545}, + title = {Chebyshev series solution of nonlinear ordinary differential equations-Initial-value problems}, + type = {techreport}, +} + @report{HairerWanner2010, author = {Hairer, Ernst and Wanner, Gerhard}, institution = {Department of Mathematics, University of Geneva, Switzerland}, @@ -2393,6 +2454,14 @@ @thesis{Chandler73 type = {phdthesis}, } +@thesis{Macomber2015, + author = {Macomber, Brent David}, + url = {https://oaktrust.library.tamu.edu/server/api/core/bitstreams/e94bc2a8-7916-4d02-9aef-e7493a07ccb9/content}, + date = {2015}, + title = {Enhancements to Chebyshev-Picard iteration efficiency for generally perturbed orbits and constrained dynamical systems}, + type = {phdthesis}, +} + @unpublished{Andrew2017, author = {“Andrew”}, date = {2017-05-10}, diff --git a/integrators/integrators_body.hpp b/integrators/integrators_body.hpp index 28010e7c34..fa9585ebb2 100644 --- a/integrators/integrators_body.hpp +++ b/integrators/integrators_body.hpp @@ -15,6 +15,7 @@ #include "integrators/methods.hpp" #include "integrators/symmetric_linear_multistep_integrator.hpp" #include "integrators/symplectic_runge_kutta_nyström_integrator.hpp" +#include "integrators/чебышёв_picard_iterator.hpp" // 🧙 For the integrator subclass. // NOLINT #include "quantities/serialization.hpp" // A case branch in a switch on the serialized integrator `kind`. It determines @@ -226,7 +227,7 @@ PRINCIPIA_INTEGRATOR_CASE(FixedStepSizeIntegrator, \ YOSHIDA_1990_ORDER_8E, \ 吉田1990Order8E, \ - sprk_action) + sprk_action) \ namespace principia { namespace integrators { @@ -544,7 +545,8 @@ FixedStepSizeIntegrator::Instance::ReadFromMessage( PRINCIPIA_READ_FSS_INTEGRATOR_INSTANCE_ERK, PRINCIPIA_READ_FSS_INTEGRATOR_INSTANCE_SLMS, PRINCIPIA_READ_FSS_INTEGRATOR_INSTANCE_SPRK, - PRINCIPIA_READ_FSS_INTEGRATOR_INSTANCE_SRKN) + PRINCIPIA_READ_FSS_INTEGRATOR_INSTANCE_SRKN + ) default: LOG(FATAL) << message.DebugString(); } @@ -613,6 +615,7 @@ FixedStepSizeIntegrator::Instance::Instance( std::abort(); \ } + template FixedStepSizeIntegrator const& FixedStepSizeIntegrator::ReadFromMessage( @@ -622,7 +625,8 @@ FixedStepSizeIntegrator::ReadFromMessage( PRINCIPIA_READ_FSS_INTEGRATOR_ERK, PRINCIPIA_READ_FSS_INTEGRATOR_SLMS, PRINCIPIA_READ_FSS_INTEGRATOR_SPRK, - PRINCIPIA_READ_FSS_INTEGRATOR_SRKN) + PRINCIPIA_READ_FSS_INTEGRATOR_SRKN + ) default: LOG(FATAL) << message.kind(); std::abort(); diff --git a/integrators/methods.hpp b/integrators/methods.hpp index de3fe21be4..deae8d7697 100644 --- a/integrators/methods.hpp +++ b/integrators/methods.hpp @@ -184,6 +184,23 @@ struct AsSymplecticRungeKuttaNyström { }; }; +// A ЧебышёвPicard integration method. +// +// M controls the number of nodes at which the function will be evaluated. +// +// Note that this is the highest node _index_ rather than the number of nodes; +// the actual number of nodes is M + 1. +// +// N is the order of the Чебышёв series used to approximate the system state. +template +struct ЧебышёвPicard : not_constructible { + static constexpr std::int64_t M = _M; + static constexpr std::int64_t N = _N; + + static_assert(M >= 1); + static_assert(N >= 1); + static_assert(M >= N); +}; struct AdamsBashforthOrder2 : ExplicitLinearMultistep { static constexpr int order = 2; @@ -1423,6 +1440,7 @@ using internal::Ruth1983; using internal::SymmetricLinearMultistep; using internal::SymplecticPartitionedRungeKutta; using internal::SymplecticRungeKuttaNyström; +using internal::ЧебышёвPicard; using internal::吉田1990Order6A; using internal::吉田1990Order6B; using internal::吉田1990Order6C; diff --git "a/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator.hpp" "b/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator.hpp" new file mode 100644 index 0000000000..4d468648b8 --- /dev/null +++ "b/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator.hpp" @@ -0,0 +1,149 @@ +// The files containing the tree of of child classes of `Integrator` must be +// included in the order of inheritance to avoid circular dependencies. This +// class will end up being reincluded as part of the implementation of its +// parent. +#ifndef PRINCIPIA_INTEGRATORS_INTEGRATORS_HPP_ +#include "integrators/integrators.hpp" +#else +#ifndef PRINCIPIA_INTEGRATORS_ЧЕБЫШЁВ_PICARD_ITERATOR_HPP_ +#define PRINCIPIA_INTEGRATORS_ЧЕБЫШЁВ_PICARD_ITERATOR_HPP_ + +#include +#include + +#include "absl/status/status.h" +#include "base/not_null.hpp" +#include "integrators/ordinary_differential_equations.hpp" +#include "numerics/fixed_arrays.hpp" +#include "quantities/quantities.hpp" +#include "serialization/integrators.pb.h" + +namespace principia { +namespace integrators { +namespace _чебышёв_picard_iterator { +namespace internal { + +using namespace principia::base::_not_null; +using namespace principia::integrators::_integrators; +using namespace principia::integrators::_ordinary_differential_equations; +using namespace principia::numerics::_fixed_arrays; +using namespace principia::quantities::_quantities; + +struct ЧебышёвPicardIterationParams { + // The maximum allowed number of Picard iterations per step. If iteration has + // not stopped (according to the stopping criterion) by the final step, the + // iteration will be considered to have diverged. + std::int64_t max_iterations; + + // If the maximum absolute difference between successive state approximations + // is less than this for two Picard iterations in a row, iteration will be + // considered to have converged. + double stopping_criterion; +}; + +// This class solves ordinary differential equations of the form x′ = f(x, t) +// using Чебышёв-Picard iteration. +// +// Чебышёв-Picard iteration combines Чебышёв interpolation with Picard +// iteration; it was first proposed in [CN63]. It works as follows: +// +// * Start with some initial approximation x⁰ (the constant function +// x⁰(t) = x(t₀) is typically used). +// * Repeatedly construct xⁱ⁺¹ from xⁱ by +// * Approximating f(t, xⁱ) using Чебышёв interpolation. +// * Approximating x₀ + ∫ₜ₀ᵗ f(t, xⁱ(t)) dt by integrating the above +// approximation. This will be xⁱ⁺¹. +// * The sequence xⁱ should converge to (a Чебышёв approximation of) x within +// some interval of t₀. +// +// Due to the properties of Чебышёв polynomials, the iteration may be peformed +// simply using linear algebra operations (in addition to the required +// evaluation of f at Чебышёв nodes) [FN83]. +// +// This code uses the formulation from [Mac15]. + +template +class ЧебышёвPicardIterator : public FixedStepSizeIntegrator { + public: + using ODE = ODE_; + using AppendState = typename Integrator::AppendState; + + class Instance : public FixedStepSizeIntegrator::Instance { + public: + absl::Status Solve(ODE::IndependentVariable const& t_final) override; + + ЧебышёвPicardIterator const& integrator() const override; + + not_null::Instance>> Clone() + const override; + + private: + Instance(InitialValueProblem const& problem, + AppendState const& append_state, + Time const& step, + ЧебышёвPicardIterator const& integrator); + + // The dimension of the ODE. + static constexpr std::int64_t n = + std::tuple_size_v; + + ЧебышёвPicardIterator const& integrator_; + + // Working variables which are stored here so we don't need to reallocate + // them on each Solve call. + + // Stores the nodes rescaled to the current step. + std::vector t_; + + // Controls the boundary condition. + FixedMatrix CₓX₀_; + + // Xⁱ is an (M + 1)×n matrix containing the values of the dependent + // variables at each node. + FixedMatrix Xⁱ_; + FixedMatrix Xⁱ⁺¹_; + + // The computed derivative (at each node, for the current iteration). + FixedMatrix yʹ_; + + friend class ЧебышёвPicardIterator; + }; + + // Constructs a ЧебышёвPicardIterator with the given parameters. + explicit ЧебышёвPicardIterator(ЧебышёвPicardIterationParams const& params); + + ЧебышёвPicardIterationParams const& params() const; + + not_null::Instance>> NewInstance( + InitialValueProblem const& problem, + AppendState const& append_state, + Time const& step) const override; + + void WriteToMessage( + not_null message) const override; + + private: + ЧебышёвPicardIterationParams params_; + + // The nodes used for function evaluation. + // + // These are Чебышёв nodes of the second kind. + FixedVector nodes_; + + // The product of 1.31a and 1.31b from [Mac15]. + FixedMatrix CₓCα_; +}; + +} // namespace internal + +using internal::ЧебышёвPicardIterationParams; +using internal::ЧебышёвPicardIterator; + +} // namespace _чебышёв_picard_iterator +} // namespace integrators +} // namespace principia + +#include "integrators/чебышёв_picard_iterator_body.hpp" + +#endif // PRINCIPIA_INTEGRATORS_ЧЕБЫШЁВ_PICARD_ITERATOR_HPP_ +#endif // PRINCIPIA_INTEGRATORS_INTEGRATORS_HPP_ diff --git "a/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator_body.hpp" "b/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator_body.hpp" new file mode 100644 index 0000000000..23434c7150 --- /dev/null +++ "b/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator_body.hpp" @@ -0,0 +1,332 @@ +#pragma once + +#include "integrators/чебышёв_picard_iterator.hpp" + +#include +#include +#include +#include + +#include "base/for_all_of.hpp" +#include "base/status_utilities.hpp" // 🧙 For RETURN_IF_ERROR. +#include "base/tags.hpp" +#include "geometry/sign.hpp" +#include "numerics/elementary_functions.hpp" +#include "numerics/matrix_views.hpp" +#include "quantities/si.hpp" + +namespace principia { +namespace integrators { + +namespace _чебышёв_picard_iterator { +namespace internal { + +using namespace principia::base::_for_all_of; +using namespace principia::base::_tags; +using namespace principia::geometry::_sign; +using namespace principia::numerics::_elementary_functions; +using namespace principia::numerics::_matrix_views; +using namespace principia::quantities::_si; + +template +ODE_::DependentVariables DependentVariablesFromMatrixRow( + FixedMatrix> const& + matrix, + std::int64_t const row) { + std::int64_t j = 0; + typename ODE_::DependentVariables y; + for_all_of(y).loop([&j, &matrix, row](auto& yⱼ) { + yⱼ = matrix(row, j++) * si::Unit>; + }); + return y; +} + +template +void DependentVariablesToMatrixRow( + typename ODE_::DependentVariables const& y, + std::int64_t const row, + FixedMatrix>& matrix) { + std::int64_t j = 0; + for_all_of(y).loop([row, &matrix, &j](auto const& yⱼ) { + matrix(row, j++) = yⱼ / si::Unit>; + }); +} + +template +void DependentVariableDerivativesToMatrixRow( + typename ODE_::DependentVariableDerivatives const& y, + std::int64_t const row, + FixedMatrix>& matrix) { + std::int64_t j = 0; + for_all_of(y).loop([row, &matrix, &j](auto const& yⱼ) { + matrix(row, j++) = yⱼ / si::Unit>; + }); +} + +// Returns max|aᵢⱼ|. +template +double LInfinityNorm(FixedMatrix const& A) { + double norm = 0.0; + for (std::int64_t i = 0; i < M; ++i) { + for (std::int64_t j = 0; j < N; ++j) { + norm = std::max(norm, std::abs(A(i, j))); + } + } + return norm; +} + +template +absl::Status ЧебышёвPicardIterator::Instance::Solve( + ODE::IndependentVariable const& t_final) { + using DependentVariableDerivatives = + typename ODE::DependentVariableDerivatives; + using State = typename ODE::State; + + auto& append_state = this->append_state_; + auto& current_state = this->current_state_; + auto const& equation = this->equation_; + auto const& step = this->step_; + auto const& params = integrator_.params(); + + // Argument checks. + Sign const integration_direction = Sign(step); + if (integration_direction.is_positive()) { + // Integrating forward. + CHECK_LT(current_state.s.value, t_final); + } else { + // Integrating backward. + CHECK_GT(current_state.s.value, t_final); + } + + while (integration_direction.is_positive() + ? current_state.s.value < t_final + : current_state.s.value > t_final) { + auto const t_initial = current_state.s.value; + + // Rescale the nodes for feeding into the compute_derivative function. + t_.clear(); + for (const double node : integrator_.nodes_) { + t_.push_back(t_initial + (0.5 * node + 0.5) * step); + } + + // Set the boundary condition and store it in CₓX₀_. + std::int64_t j = 0; + for_all_of(current_state.y).loop([this, &j](auto const& yⱼ) { + CₓX₀_(0, j++) = yⱼ.value / si::Unit; + }); + for (std::int64_t i = 1; i <= Method::M; ++i) { + for (std::int64_t j = 0; j < n; ++j) { + CₓX₀_(i, j) = CₓX₀_(0, j); + } + } + + // A good starting guess for X⁰ is uniform current_state.y; as it happens + // that's what we just set CₓX₀_ to. + Xⁱ_ = CₓX₀_; + + double previous_norm = std::numeric_limits::infinity(); + bool converged = false; + for (int64_t iteration = 0; iteration < params.max_iterations; + ++iteration) { + // Evaluate the right hand side of the equation. + for (int64_t i = 0; i < Xⁱ_.rows(); ++i) { + auto const y = + DependentVariablesFromMatrixRow(Xⁱ_, i); + DependentVariableDerivatives yʹᵢ; + RETURN_IF_ERROR(equation.compute_derivative(t_[i], y, yʹᵢ)); + + // Store it in yʹ. + DependentVariableDerivativesToMatrixRow( + yʹᵢ, i, yʹ_); + } + + // Compute new x. + Xⁱ⁺¹_ = integrator_.CₓCα_ * (0.5 * step / Second * yʹ_) + CₓX₀_; + + // Check for convergence by computing the ∞-norm. + const double norm = LInfinityNorm(Xⁱ⁺¹_ - Xⁱ_); + Xⁱ_ = std::move(Xⁱ⁺¹_); + + // We require that ‖Xⁱ⁺¹ - Xⁱ‖ and ‖Xⁱ - Xⁱ⁻¹‖ are _both_ less than + // the given tolerance to account for nonlinearity issues (as suggested in + // [BJ12]). + if (std::max(norm, previous_norm) < params.stopping_criterion) { + converged = true; + break; + } + + previous_norm = norm; + RETURN_IF_STOPPED; + } + + if (converged) { + // We have successfully converged! + for (std::int64_t i = 0; i < Xⁱ_.rows(); ++i) { + append_state( + State(t_[i], DependentVariablesFromMatrixRow(Xⁱ_, i))); + } + + // Set the current state to the final state we appended. + current_state = + State(t_[Xⁱ_.rows() - 1], + DependentVariablesFromMatrixRow(Xⁱ_, Xⁱ_.rows() - 1)); + RETURN_IF_STOPPED; + } else { + // We failed to converge. + return absl::Status(absl::StatusCode::kFailedPrecondition, + "Чебышёв-Picard iteration failed to converge."); + } + } + + return absl::OkStatus(); +} + +template +ЧебышёвPicardIterator const& +ЧебышёвPicardIterator::Instance::integrator() const { + return integrator_; +} + +template +not_null::Instance>> +ЧебышёвPicardIterator::Instance::Clone() const { + return std::unique_ptr(new Instance(*this)); +} + +template +ЧебышёвPicardIterator::Instance::Instance( + InitialValueProblem const& problem, + AppendState const& append_state, + Time const& step, + ЧебышёвPicardIterator const& integrator) + : FixedStepSizeIntegrator::Instance(problem, append_state, step), + integrator_(integrator), + CₓX₀_(uninitialized), + Xⁱ_(uninitialized), + Xⁱ⁺¹_(uninitialized), + yʹ_(uninitialized) { + t_.reserve(integrator.nodes_.size()); +} + +template +ЧебышёвPicardIterator::ЧебышёвPicardIterator( + ЧебышёвPicardIterationParams const& params) + : params_(params), nodes_(uninitialized), CₓCα_(uninitialized) { + // We use the notation from [Mac15], section 1.4.3. + constexpr std::int64_t M = Method::M; + constexpr std::int64_t N = Method::N; + + // Populate nodes. + for (std::int64_t i = 0; i <= M; i++) { + nodes_[i] = -Cos(π / M * i * Radian); + } + + // ᵝT is a (M + 1)×(N + 1) matrix of Чебышёв polynomials evaluated at nodes. + // See [Mac15], equation (1.20). + FixedMatrix ᵝT(uninitialized); + + for (std::int64_t i = 0; i <= M; ++i) { + const auto τᵢ = nodes_[i]; + // The 0-degree polynomial is uniformly 1. + ᵝT(i, 0) = 1; + // The 0-degree polynomial is the identity. + ᵝT(i, 1) = τᵢ; + + // We populate the rest of ᵝT using the recurrence relation. + for (std::int64_t j = 2; j <= N; ++j) { + ᵝT(i, j) = 2 * τᵢ * ᵝT(i, j - 1) - ᵝT(i, j - 2); + } + } + + // ᵝW is a diagonal (N + 1)×(N + 1) matrix with diagonal [½, 1, 1, ..., ½]. + // See [Mac15], equation (1.20). + FixedMatrix ᵝW; + ᵝW(0, 0) = 0.5; + ᵝW(N, N) = 0.5; + for (std::int64_t i = 1; i < N; ++i) { + ᵝW(i, i) = 1; + } + + FixedMatrix Cₓ = ᵝT * ᵝW; + + // R is a diagonal (N + 1)×(N + 1) matrix. + // See [Mac15], equation (1.25). + FixedMatrix R; + R(0, 0) = 1; + R(N, N) = 1.0 / N; + for (std::int64_t i = 1; i < N; ++i) { + R(i, i) = 1.0 / (2 * i); + } + + // S is an (N + 1)×N matrix. + // See equation 1.26 in [Mac15]. + FixedMatrix S; + S(0, 0) = 1; + S(0, 1) = -0.5; + for (std::int64_t k = 2; k < N; ++k) { + S(0, k) = (k % 2 == 1 ? 1 : -1) * (1.0 / (k - 1) - 1.0 / (k + 1)); + } + for (std::int64_t i = 0; i < N; ++i) { + S(i + 1, i) = 1; + } + for (std::int64_t i = 1; i + 2 < N; ++i) { + S(i, i + 1) = -1; + } + + // ᶠTᵀ is ᵝTᵀ with the last row removed. + // See [Mac15], equation (1.22). + FixedMatrix ᶠTᵀ(uninitialized); + for (std::int64_t i = 0; i < N; ++i) { + for (std::int64_t j = 0; j <= M; ++j) { + ᶠTᵀ(i, j) = ᵝT(j, i); + } + } + + // V is is a diagonal (M + 1)×(M + 1) matrix with diagonal [1/M, 2/M, 2/M, + // ..., 1/M]. + FixedMatrix V; + V(0, 0) = V(M, M) = 1.0 / M; + for (std::int64_t i = 1; i < M; ++i) { + V(i, i) = 2.0 / M; + } + + // Cα is R * R * ᶠTᵀ * V (we do not assign it to a variable). + + CₓCα_ = Cₓ * R * S * ᶠTᵀ * V; +} + +template +ЧебышёвPicardIterationParams const& +ЧебышёвPicardIterator::params() const { + return params_; +} + +template +not_null::Instance>> +ЧебышёвPicardIterator::NewInstance( + InitialValueProblem const& problem, + AppendState const& append_state, + Time const& step) const { + // Cannot use `make_not_null_unique` because the constructor of `Instance` is + // private. + return std::unique_ptr( + new Instance(problem, append_state, step, *this)); +} + +template +void ЧебышёвPicardIterator::WriteToMessage( + not_null message) const { + LOG(FATAL) << "Body is neither massive nor massless"; + std::abort(); +} + +} // namespace internal +} // namespace _чебышёв_picard_iterator +} // namespace integrators +} // namespace principia diff --git "a/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator_test.cpp" "b/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator_test.cpp" new file mode 100644 index 0000000000..24a7273944 --- /dev/null +++ "b/integrators/\321\207\320\265\320\261\321\213\321\210\321\221\320\262_picard_iterator_test.cpp" @@ -0,0 +1,442 @@ +#include "integrators/чебышёв_picard_iterator.hpp" + +#include +#include +#include + +#include "geometry/instant.hpp" +#include "gmock/gmock.h" +#include "gtest/gtest.h" +#include "integrators/methods.hpp" +#include "integrators/ordinary_differential_equations.hpp" +#include "numerics/elementary_functions.hpp" +#include "quantities/quantities.hpp" +#include "quantities/si.hpp" +#include "testing_utilities/almost_equals.hpp" +#include "testing_utilities/matchers.hpp" + +namespace principia { +namespace integrators { + +using ::testing::DoubleNear; +using ::testing::Test; +using ::testing::Types; +using ::testing::Values; + +using namespace principia::geometry::_instant; +using namespace principia::integrators::_methods; +using namespace principia::integrators::_ordinary_differential_equations; +using namespace principia::integrators::_чебышёв_picard_iterator; +using namespace principia::numerics::_elementary_functions; +using namespace principia::quantities::_quantities; +using namespace principia::quantities::_si; +using namespace principia::testing_utilities::_almost_equals; +using namespace principia::testing_utilities::_matchers; + +using ODE = ExplicitFirstOrderOrdinaryDifferentialEquation; + +namespace { + +// An initial value problem with a known solution. +struct SolvedInitialValueProblem { + ODE::IndependentVariable t₀() const { + return problem.initial_state.s.value; + } + + InitialValueProblem problem; + std::function + solution; +}; + +// The first order ODE y′ = y with y₀ = 1. +// +// The solution is y = eᵗ. +SolvedInitialValueProblem LinearProblem() { + ODE linear_ode; + linear_ode.compute_derivative = + [](Instant const& t, + ODE::DependentVariables const& dependent_variables, + ODE::DependentVariableDerivatives& dependent_variable_derivatives) { + auto const& [y] = dependent_variables; + auto& [yʹ] = dependent_variable_derivatives; + yʹ = y / Second; + return absl::OkStatus(); + }; + + InitialValueProblem problem; + problem.equation = linear_ode; + problem.initial_state = {Instant(), {1 * Metre}}; + + return SolvedInitialValueProblem{ + .problem = problem, + .solution = [](Instant const& t) -> std::tuple { + return std::exp((t - Instant()) / Second) * Metre; + }}; +} + +// The first order ODE y′ = π/8 (1 + y²) with y(-1) = 0. +// +// The solution is y = tan(π/8 (t + 1)). +SolvedInitialValueProblem TangentProblem() { + ODE ode; + ode.compute_derivative = + [](Instant const& t, + ODE::DependentVariables const& dependent_variables, + ODE::DependentVariableDerivatives& dependent_variable_derivatives) { + auto const& [y] = dependent_variables; + auto& [yʹ] = dependent_variable_derivatives; + yʹ = π / 8 * (1 + y * y / (Metre * Metre)) * Metre / Second; + return absl::OkStatus(); + }; + + InitialValueProblem problem; + problem.equation = ode; + problem.initial_state = {Instant() - 1 * Second, {0 * Metre}}; + + return SolvedInitialValueProblem{ + .problem = problem, + .solution = [](Instant const& t) -> std::tuple { + return Tan(π / 8 * ((t - Instant()) / Second + 1) * Radian) * Metre; + }}; +} + +// The first order ODE y′ = cos(t + εy). +// +// The solution ([Fuk97], equations 31 and 32) is +// y = -γt + 2/ε tan⁻¹((β + σ cos φ)/(1 + σ sin φ)) +// where +// σ = α(sin φ + β cos φ) +// φ = ½ (1 - γε) t +// α = 2ε/(1 - ε + √(1 - ε²)) +// β = 1/(1 + α) tan(εy₀/2) +// γ = ε/(1 + √(1 - ε²)) +SolvedInitialValueProblem PerturbedSinusoidProblem(double const ε, + double const y₀) { + ODE ode; + ode.compute_derivative = + [ε](Instant const& t, + ODE::DependentVariables const& dependent_variables, + ODE::DependentVariableDerivatives& dependent_variable_derivatives) { + auto const& [y] = dependent_variables; + auto& [yʹ] = dependent_variable_derivatives; + yʹ = Cos(((t - Instant()) / Second + ε * y / Metre) * Radian) * Metre / + Second; + return absl::OkStatus(); + }; + + InitialValueProblem problem; + problem.equation = ode; + problem.initial_state = {Instant(), {y₀ * Metre}}; + + const double one_plus_sqrt_one_minus_ε² = 1 + Sqrt(1 - ε * ε); + const double α = 2 * ε / (one_plus_sqrt_one_minus_ε² - ε); + const double β = 1 / (1 + α) * Tan(ε * y₀ / 2 * Radian); + const double γ = ε / one_plus_sqrt_one_minus_ε²; + + return SolvedInitialValueProblem{ + .problem = problem, + .solution = [ε, α, β, γ](Instant const& t) -> std::tuple { + const Angle φ = 0.5 * (1 - γ * ε) * (t - Instant()) / Second * Radian; + const double σ = α * (Sin(φ) + β * Cos(φ)); + + return (-γ * (t - Instant()) / Second + + 2 / ε * ArcTan((β + σ * Cos(φ)) / (1 + σ * Sin(φ))) / Radian) * + Metre; + }}; +} + +TEST(ЧебышёвPicardIteratorTest, MultipleSteps) { + // Set up the initial value problem. + SolvedInitialValueProblem const& problem = LinearProblem(); + Time const step = 1 * Second; + + std::vector states; + auto const append_state = [&states](ODE::State const& state) { + states.push_back(state); + }; + + // Build the integrator and solve the problem. + ЧебышёвPicardIterator<ЧебышёвPicard<64, 64>, ODE> const integrator( + ЧебышёвPicardIterationParams{ + .max_iterations = 64, + .stopping_criterion = 1e-16, + }); + + auto const instance = + integrator.NewInstance(problem.problem, append_state, step); + auto const t_final = problem.t₀() + 2.5 * step; + + EXPECT_OK(instance->Solve(t_final)); + + // Verify we have reached the desired final time. + EXPECT_GE(states.back().s.value, t_final); + + // Verify the results are close to the known solution. + for (const auto& state : states) { + auto t = state.s.value; + auto y = std::get<0>(state.y).value; + EXPECT_THAT(y, AlmostEquals(std::get<0>(problem.solution(t)), 0, 5)) + << "t=" << t; + } +} + +TEST(ЧебышёвPicardIteratorTest, Backwards) { + // Set up the initial value problem. + SolvedInitialValueProblem const& problem = LinearProblem(); + Time const step = -1 * Second; // Backwards! + + std::vector states; + auto const append_state = [&states](ODE::State const& state) { + states.push_back(state); + }; + + // Build the integrator and solve the problem. + ЧебышёвPicardIterator<ЧебышёвPicard<64, 64>, ODE> const integrator( + ЧебышёвPicardIterationParams{ + .max_iterations = 64, + .stopping_criterion = 2e-16, + }); + + auto const instance = + integrator.NewInstance(problem.problem, append_state, step); + auto const t_final = problem.t₀() + 2.5 * step; + + EXPECT_OK(instance->Solve(t_final)); + + // Verify we have reached the desired final time. + EXPECT_LE(states.back().s.value, t_final); // ≤ because we are backwards. + + // Verify the results are close to the known solution. + for (const auto& state : states) { + auto t = state.s.value; + auto y = std::get<0>(state.y).value; + EXPECT_THAT(y, AlmostEquals(std::get<0>(problem.solution(t)), 0, 12)) + << "t=" << t; + } +} + +TEST(ЧебышёвPicardIteratorTest, Divergence) { + // Set up the initial value problem. + SolvedInitialValueProblem const& problem = LinearProblem(); + Time const step = 60 * Second; // Way too big; iteration won't converge. + + std::vector states; + auto const append_state = [&states](ODE::State const& state) { + states.push_back(state); + }; + + // Build the integrator and solve the problem. + ЧебышёвPicardIterator<ЧебышёвPicard<64, 64>, ODE> const integrator( + ЧебышёвPicardIterationParams{ + .max_iterations = 64, + .stopping_criterion = 0.1, // Differences never even get this low! + }); + + auto const instance = + integrator.NewInstance(problem.problem, append_state, step); + auto const t_final = problem.t₀() + step; + + EXPECT_THAT(instance->Solve(t_final), + StatusIs(absl::StatusCode::kFailedPrecondition)); +} + +template +concept ЧебышёвPicardIteratorTestParameters = requires { + // The ODE (with solution) under test. + { T::problem() } -> std::convertible_to; + + // The ЧебышёвPicard method. + typename T::Method; + + // The step size used for this test. + { T::step } -> std::convertible_to; + + // The stopping criterion used for this test. + { T::stopping_criterion } -> std::convertible_to; + + // The distance the integrated solution is allowed to vary from the analytic + // solution. + { T::tolerance } -> std::convertible_to; +}; + +template<ЧебышёвPicardIteratorTestParameters T> +class ЧебышёвPicardIteratorParameterizedTest : public Test {}; + +TYPED_TEST_SUITE_P(ЧебышёвPicardIteratorParameterizedTest); + +TYPED_TEST_P(ЧебышёвPicardIteratorParameterizedTest, Convergence) { + // Set up the initial value problem. + SolvedInitialValueProblem const problem = TypeParam::problem(); + Time const step = TypeParam::step; + + std::vector states; + auto const append_state = [&states](ODE::State const& state) { + states.push_back(state); + }; + + // Build the integrator and solve the problem. + ЧебышёвPicardIterator const integrator( + ЧебышёвPicardIterationParams{ + .max_iterations = 64, + .stopping_criterion = TypeParam::stopping_criterion, + }); + + auto const instance = + integrator.NewInstance(problem.problem, append_state, step); + EXPECT_OK(instance->Solve(problem.t₀() + step)); + + // Verify the results are close to the known solution. + for (const auto& state : states) { + auto t = state.s.value; + auto y = std::get<0>(state.y).value; + EXPECT_THAT(y / Metre, + DoubleNear(std::get<0>(problem.solution(t)) / Metre, + TypeParam::tolerance)) + << "t=" << t; + } +} + +REGISTER_TYPED_TEST_SUITE_P(ЧебышёвPicardIteratorParameterizedTest, + Convergence); + +// Although in theory the iteration for y′ = y should converge for +// intervals +// of length < 40 for sufficiently high N (see [BJ12]), in +// practice the convergence degrades far earlier. +struct Linear1Second : not_constructible { + static SolvedInitialValueProblem problem() { + return LinearProblem(); + } + + using Method = ЧебышёвPicard<64>; + static constexpr Time step = 1 * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 9e-16; +}; +struct Linear2Seconds : not_constructible { + static SolvedInitialValueProblem problem() { + return LinearProblem(); + } + + using Method = ЧебышёвPicard<64>; + static constexpr Time step = 2 * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 4.5e-15; +}; +struct Linear4Seconds : not_constructible { + static SolvedInitialValueProblem problem() { + return LinearProblem(); + } + + using Method = ЧебышёвPicard<64>; + static constexpr Time step = 4 * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 8e-14; +}; +struct Linear8Seconds : not_constructible { + static SolvedInitialValueProblem problem() { + return LinearProblem(); + } + + using Method = ЧебышёвPicard<64>; + static constexpr Time step = 8 * Second; + static constexpr double stopping_criterion = 1e-12; + static constexpr double tolerance = 4e-10; +}; +struct Linear16Seconds : not_constructible { + static SolvedInitialValueProblem problem() { + return LinearProblem(); + } + + using Method = ЧебышёвPicard<64>; + static constexpr Time step = 16 * Second; + static constexpr double stopping_criterion = 1e-7; + // Absolute error is high due to the + // exponential growth. + static constexpr double tolerance = 4e-3; +}; +struct LinearMGreaterThanN : not_constructible { + static SolvedInitialValueProblem problem() { + return LinearProblem(); + } + + using Method = ЧебышёвPicard<128, 64>; + static constexpr Time step = 1 * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 2e-15; +}; + +using Linear = Types; +INSTANTIATE_TYPED_TEST_SUITE_P(Linear, + ЧебышёвPicardIteratorParameterizedTest, + Linear); + +// This problem appears in [BL69]. +struct TangentA : not_constructible { + static SolvedInitialValueProblem problem() { + return TangentProblem(); + } + + // These are the parameters from [BL69]. + using Method = ЧебышёвPicard<16>; + static constexpr Time step = 2 * Second; + static constexpr double stopping_criterion = 0.5e-10; + static constexpr double tolerance = 8.7e-10; +}; +struct TangentB : not_constructible { + static SolvedInitialValueProblem problem() { + return TangentProblem(); + } + + // We can achieve better accuracy with higher N and + // a more stringent stopping criterion. + using Method = ЧебышёвPicard<32>; + static constexpr Time step = 2 * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 4.5e-16; +}; + +using Tangent = Types; +INSTANTIATE_TYPED_TEST_SUITE_P(Tangent, + ЧебышёвPicardIteratorParameterizedTest, + Tangent); + +// From [BJ10]. Figures 4, 5, 7, 8 are slow and omitted for now. +struct PerturbedSinusoidFigure3 : not_constructible { + static SolvedInitialValueProblem problem() { + return PerturbedSinusoidProblem(/*ε=*/0.01, + /*y₀=*/1); + } + + using Method = ЧебышёвPicard<500>; + static constexpr Time step = 64 * π * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 1e-9; +}; +struct PerturbedSinusoidFigure6 : not_constructible { + static SolvedInitialValueProblem problem() { + return PerturbedSinusoidProblem(/*ε=*/0.001, + /*y₀=*/1); + } + + using Method = ЧебышёвPicard<400>; + static constexpr Time step = 64 * π * Second; + static constexpr double stopping_criterion = 1e-16; + static constexpr double tolerance = 2.5e-11; +}; + +using PerturbedSinusoid = + Types; +INSTANTIATE_TYPED_TEST_SUITE_P(PerturbedSinusoid, + ЧебышёвPicardIteratorParameterizedTest, + PerturbedSinusoid); + +} // namespace + +} // namespace integrators +} // namespace principia diff --git a/serialization/integrators.proto b/serialization/integrators.proto index a809ab13ff..1b16253852 100644 --- a/serialization/integrators.proto +++ b/serialization/integrators.proto @@ -29,6 +29,7 @@ message FixedStepSizeIntegrator { extend Integrator { optional FixedStepSizeIntegrator extension = 3001; } + // Next ID: 44 enum Kind { ADAMS_BASHFORTH_ORDER_2 = 39; ADAMS_BASHFORTH_ORDER_3 = 40; @@ -168,3 +169,11 @@ message State { repeated DoublePrecision velocity = 2; required DoublePrecision time = 3; } + +// This corresponds to the ChebyshevPicardIterationParams struct. +message ChebyshevPicardIterationParams { + optional int32 m = 1; + optional int32 n = 2; + optional int32 max_iterations = 3; + optional double stopping_criterion = 4; +}