Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
44 commits
Select commit Hold shift + click to select a range
dfecc88
MCPI: it's alive!
rnlahaye Jan 12, 2026
b95ba19
ChebyshevPicardIterationParams.
rnlahaye Jan 13, 2026
d3585e7
ChebyshevPicardIteratorTest TEST_P.
rnlahaye Jan 13, 2026
31bb19c
More test cases for ChebyshevPicardIteratorTest.
rnlahaye Jan 13, 2026
6143f45
Move the solved IVP into the test parameters.
rnlahaye Jan 13, 2026
7103c72
Add the tangent problem from the NASA paper.
rnlahaye Jan 13, 2026
e908836
Add the perturbed sinusoid problem.
rnlahaye Jan 13, 2026
5e1cf93
MultipleSteps test.
rnlahaye Jan 13, 2026
f03fd4f
ChebyshevPicardIterator: backwards stepping.
rnlahaye Jan 13, 2026
5f9bb7b
ChebyshevPicardIteratorTest divergence test.
rnlahaye Jan 13, 2026
2bc7369
Implement ChebyshevPicardIterator WriteToMessage.
rnlahaye Jan 14, 2026
c4429b7
Integrator #ifndef trick for ChebyshevPicardIterator.
rnlahaye Jan 14, 2026
72b52df
Go back to absolute tolerances in CPI test.
rnlahaye Jan 14, 2026
cb02b0e
Use -1 as starting t for tangent test now that we are using absolute …
rnlahaye Jan 14, 2026
37e426b
Rename variables to better match the source.
rnlahaye Jan 14, 2026
7895c44
Function for MaxNorm.
rnlahaye Jan 14, 2026
43eef1b
Don't store Cₓ; it isn't really necessary.
rnlahaye Jan 14, 2026
3c22631
Rename x to Xⁱ.
rnlahaye Jan 14, 2026
a73c0d2
Try not to allocate in the loops.
rnlahaye Jan 14, 2026
0192266
Rename chebyshev to чебышёв in filenames.
rnlahaye Jan 15, 2026
b957136
Rename Chebyshev to Чебышёв inside the files.
rnlahaye Jan 15, 2026
ac29708
Documentation and citations.
rnlahaye Jan 15, 2026
2849aa2
Remove the norm printing and the decomposition include.
rnlahaye Jan 15, 2026
6c5925a
Move Fukushima1997 to the correct position.
rnlahaye Jan 16, 2026
afe326d
Lint.
rnlahaye Jan 16, 2026
8c9c620
Handle CHEBYSHEV_PICARD case in ReadFromMessage.
rnlahaye Jan 16, 2026
1965e8d
More lint.
rnlahaye Jan 16, 2026
a346203
Updates to bibliography as requested.
rnlahaye Jan 27, 2026
1a66062
Address some comments.
rnlahaye Jan 27, 2026
6bbbaf2
Address more comments.
rnlahaye Jan 27, 2026
16bcf39
L∞Norm -> LInfinity norm (fix Windows build).
rnlahaye Jan 27, 2026
56a2a2e
Add ЧебышёвPicard method.
rnlahaye Feb 5, 2026
903079c
Change more i++ to ++i (also j, k).
rnlahaye Feb 5, 2026
a11e5b5
Parameterize ЧебышёвPicardIteratorParameterizedTest by method instead…
rnlahaye Feb 5, 2026
0f235d7
Fix computation of CₓCα_, and add a test for M > N.
rnlahaye Feb 5, 2026
fed4c5f
Use FixedMatrix instead of UnboundedMatrix in ЧебышёвPicardIterator.
rnlahaye Feb 5, 2026
ec46317
Don’t escape filenames when listing files changed by include_what_you…
eggrobin Feb 5, 2026
17be53f
[StringSplitOptions]::RemoveEmptyEntries
eggrobin Feb 6, 2026
bfb2881
Try to make it work when it works
eggrobin Feb 6, 2026
d320f0e
IWYU.
rnlahaye Feb 6, 2026
0139a50
Spell Nacozy's name correctly.
rnlahaye Feb 6, 2026
abf6275
Remove ЧебышёвPicard serialization for now.
rnlahaye Feb 6, 2026
697bf9a
Use ‖ instead of ||.
rnlahaye Feb 6, 2026
6b30b63
Merge branch 'master' into mcpi
rnlahaye Feb 6, 2026
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 8 additions & 5 deletions .github/workflows/lint.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -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
69 changes: 69 additions & 0 deletions documentation/bibliography.bib
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand Down Expand Up @@ -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},
Expand Down Expand Up @@ -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},
Expand Down Expand Up @@ -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},
Expand Down Expand Up @@ -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},
Expand Down Expand Up @@ -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},
}
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While the published paper contains the second-order vector/matrix formulation, it does not contain the first-order one (which is what I'm using here).

I added a URL to the citation.


@unpublished{Andrew2017,
author = {“Andrew”},
date = {2017-05-10},
Expand Down
10 changes: 7 additions & 3 deletions integrators/integrators_body.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -226,7 +227,7 @@
PRINCIPIA_INTEGRATOR_CASE(FixedStepSizeIntegrator, \
YOSHIDA_1990_ORDER_8E, \
吉田1990Order8E, \
sprk_action)
sprk_action) \

namespace principia {
namespace integrators {
Expand Down Expand Up @@ -544,7 +545,8 @@ FixedStepSizeIntegrator<ODE_>::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();
}
Expand Down Expand Up @@ -613,6 +615,7 @@ FixedStepSizeIntegrator<ODE_>::Instance::Instance(
std::abort(); \
}


template<typename ODE_>
FixedStepSizeIntegrator<ODE_> const&
FixedStepSizeIntegrator<ODE_>::ReadFromMessage(
Expand All @@ -622,7 +625,8 @@ FixedStepSizeIntegrator<ODE_>::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();
Expand Down
18 changes: 18 additions & 0 deletions integrators/methods.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<std::int64_t _M, std::int64_t _N = _M>
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;
Expand Down Expand Up @@ -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;
Expand Down
149 changes: 149 additions & 0 deletions integrators/чебышёв_picard_iterator.hpp
Original file line number Diff line number Diff line change
@@ -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 <memory>
#include <vector>

#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;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let's have an include for this.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The include is near the top, line 6 (like the other integrators).

using namespace principia::integrators::_ordinary_differential_equations;
using namespace principia::numerics::_fixed_arrays;
using namespace principia::quantities::_quantities;

struct ЧебышёвPicardIterationParams {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So we generally distinguish the method (see methods.hpp) which is used to parameterize the integrator, from the parameters which are used when creating instances.

For instance for the ephemeris, the method is QuinlanTremaine1990Order12 and the instances are created with a step of 10 minutes.

It's clear that max_iterations and stopping_criterion are instance parameters, but I wonder if M and N would make sense as methods? WDYT?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The working matrices are fixed once M and N are chosen so it is convenient to make them part of the method and store the matrices in the class (of course, this could be done in other ways as well).

I see that all the methods in methods.hpp have fixed values, which makes sense for them but maybe less so here (where any M ≥ N is valid). One possibility is to have

struct ЧебышёвPicard {
  std::int64_t M;
  std::int64_t N;
};

(in methods.hpp)

alternatively one could have

template<std::int64_t M, std::int64_t N>
struct ЧебышёвPicard: not_constructible {};

(or am I misunderstanding?)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think template<std::int64_t M, std::int64_t N>struct ЧебышёвPicard would be most consistent with the way the others work, assuming that various combinations of M and N are interesting (for AdamsBashforth we just went with 5 separate structs so that we don’t need to implement the computation of the coefficients, but that would probably get silly here).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, added template<std::int64_t M, std::int64_t N>struct ЧебышёвPicard.

This enabled me to switch from UnboundedMatrix to FixedMatrix, which resulted in me finding and fixing a bug in the M > N case.

// 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<typename Method, typename ODE_>
class ЧебышёвPicardIterator : public FixedStepSizeIntegrator<ODE_> {
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am a bit unhappy that this is called ...Iterator and not ...Integrator. I realize that this technique is called "iteration", but really this class is an integrator in the sense of the Principia abstractions. Furthermore, "iterator" makes me think of STL iterators or other objects that are used to iterate over some collection.

(If we decide to change the name, let's do it last.)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Happy to change the name after everything else is resolved.

public:
using ODE = ODE_;
using AppendState = typename Integrator<ODE>::AppendState;

class Instance : public FixedStepSizeIntegrator<ODE>::Instance {
public:
absl::Status Solve(ODE::IndependentVariable const& t_final) override;

ЧебышёвPicardIterator const& integrator() const override;

not_null<std::unique_ptr<typename Integrator<ODE>::Instance>> Clone()
const override;

private:
Instance(InitialValueProblem<ODE> 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<typename ODE::DependentVariables>;

Чебышёв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<typename ODE::IndependentVariable> t_;

// Controls the boundary condition.
FixedMatrix<double, Method::M + 1, n> CₓX₀_;

// Xⁱ is an (M + 1)×n matrix containing the values of the dependent
// variables at each node.
FixedMatrix<double, Method::M + 1, n> Xⁱ_;
FixedMatrix<double, Method::M + 1, n> Xⁱ⁺¹_;

// The computed derivative (at each node, for the current iteration).
FixedMatrix<double, Method::M + 1, n> yʹ_;

friend class ЧебышёвPicardIterator;
};

// Constructs a ЧебышёвPicardIterator with the given parameters.
explicit ЧебышёвPicardIterator(ЧебышёвPicardIterationParams const& params);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As mentioned above, some of the things in params should be passed to the instance, I think.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How does this work with NewInstance? As defined in FixedStepSizeIntegrator, it doesn't have any spare arguments to pass them (or am I misunderstanding something?)


ЧебышёвPicardIterationParams const& params() const;

not_null<std::unique_ptr<typename Integrator<ODE>::Instance>> NewInstance(
InitialValueProblem<ODE> const& problem,
AppendState const& append_state,
Time const& step) const override;

void WriteToMessage(
not_null<serialization::FixedStepSizeIntegrator*> message) const override;

private:
ЧебышёвPicardIterationParams params_;

// The nodes used for function evaluation.
//
// These are Чебышёв nodes of the second kind.
FixedVector<double, Method::M + 1> nodes_;

// The product of 1.31a and 1.31b from [Mac15].
FixedMatrix<double, Method::M + 1, Method::M + 1> 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_
Loading