Чебышёв-Picard Iteration by rnlahaye · Pull Request #4409 · mockingbirdnest/Principia

rnlahaye · 2026-01-16T06:17:51Z

This change adds an Integrator subclass which implements Чебышёв-Picard iteration.

Чебышёв-Picard iteration (a combination of Чебышёв approximation with Picard iteration) was first proposed by Clenshaw and Norton in 1963. It has since been studied by Shaver, Feagin and Narcozy, Fukushima, and recently by Junkins et al (see also this NASA report from 1969 which provides an implementation in FORTRAN).

Given a first order ordinary differential equation $x'=f(x,t)$ with initial condition $x(t_0)=x_0$, it can be used to solve the ODE on an interval $[t_0, t_f]$:

Start with some initial approximation $x^0$ (typically the constant $x^0(t) = x_0$ is used).
For each $x^i$:
- Compute a Чебышёв approximation $g$ to $f(x^i, t)$ over $[t_0, t_f]$ (this requires evaluating $f(x^i, t)$ at Чебышёв nodes).
- Compute $x^{i+1}(t) = x_0 + \int_{t_0}^tg(s)\ ds$. Note that this integral can be evaluated analytically and yields another Чебышёв series.
For sufficiently small intervals (and well-behaved $f$), the sequence $x^i$ converges uniformly.

See the linked papers for more details.

There is evidence that, together with some of the optimizations below, it has better performance than… some other integrators which Principia does not use.

Regarding the implementation

The linear algebra code uses the formulation (and notation) from Macomber's thesis.

I made the integrator a subclass of FixedStepSizeIntegrator because it has a size of something like a step which is fixed. I implemented WriteToMessage but not ReadFromMessage as I wasn't sure of the best way to return a const reference while avoiding a memory leak.

Regarding future work

The following is yet to be implemented:

Improvements to the linear algebra code (by adding lazy operations a la Eigen).
Support for second order ODEs.
Adaptive step size/Чебышёв order. Macomber's thesis hints at how this can be done; there is also a paper on this which I don't have access to.
Terminal convergence approximation (discussed in Macomber's thesis and the subsequent published paper). A local Taylor series gravity approximation can be used in later iterations, reducing the cost of the gravity function evaluation.
Constraint restoration (discussed in Macomber's thesis). While Чебышёв-Picard iteration naturally preserves the Hamiltonian, it is possible to use conserved quantities to improve the convergence rate and enlarge the convergence domain.

…tolerance.

pleroy · 2026-01-16T10:32:06Z

Thanks for the contribution. It's going to take a while to review because (1) it requires a functioning brain and (2) we are in the middle of a move.

In the meantime, can you fix the lint/compilation errors?

rnlahaye · 2026-01-24T19:01:48Z

I'm not sure how to fix the check-iwyu failure; it seems to be complaining that the file doesn't exist.

pleroy · 2026-01-24T19:22:00Z

@eggrobin: Can you take a look at the check-iwyu failure? It seems to happen as part of the Powershell script that reports "something changed". I have a hunch that this is not Unicode-aware: it may be the first time we have an IWYU failure in a file with a non-ASCII name.

@rnlahaye: I'll try to do a first pass over the PR tomorrow, but it won't be particularly deep. I'd like to read the papers you linked, but that will take some time.

pleroy

A very shallow first review (I didn't look at the test) to get things going. Apologies for the million nitpicks.

pleroy · 2026-01-26T17:34:56Z

documentation/bibliography.bib

@@ -1,4 +1,4 @@
-
+


What changed here?

I'm not sure; I checked in a hex editor and the first few bytes of the files were identical.

pleroy · 2026-01-26T17:40:48Z

documentation/bibliography.bib

+  author       = {Clenshaw, C. W. and Norton, H. J.},
+  publisher    = {Oxford University Press},
+  date         = {1963},
+  doi          = {10.2/3},


The DOI should be 10.1093/comjnl/6.1.88 I think.

pleroy · 2026-01-26T17:42:08Z

documentation/bibliography.bib

+  pages        = {107--115},
+  title        = {Matrix formulation of the Picard method for parallel computation},
+  volume       = {29},
+}


Let's add the DOI 10.1007/BF01232802.

pleroy · 2026-01-26T17:42:43Z

documentation/bibliography.bib

+@article{Fukushima1997,
+  author       = {Fukushima, Toshio},
+  date         = {1997},
+  journaltitle = {Astronomical Journal v. 113, p. 1909-1914 (1997)},


Remove volume, pages, date from this.

pleroy · 2026-01-26T17:43:59Z

documentation/bibliography.bib

+  pages        = {1909--1914},
+  title        = {Picard iteration method, Chebyshev polynomial approximation, and global numerical integration of dynamical motions},
+  volume       = {113},
+}


Add the DOI 10.1086/118404.

pleroy · 2026-01-26T18:35:07Z

integrators/чебышёв_picard_iterator_body.hpp

+                DependentVariablesFromMatrixRow<ODE>(Xⁱ_, Xⁱ_.rows() - 1));
+    } else {
+      // We failed to converge.
+      return absl::Status(absl::StatusCode::kFailedPrecondition,


Would it make sense to use this status? There is a magic long distance coupling where the status that gets returned here determines the explanation given to the user in the UI, so we must make sure that we pick the right status.

The comment of that status is "The integration may be retried with the same arguments and progress will happen", which is not the case here—if the iteration doesn't converge on some interval, it will repeatedly fail to converge unless the interval is made shorter.

pleroy · 2026-01-26T18:37:07Z

integrators/чебышёв_picard_iterator_body.hpp

+    }
+  }
+
+  return absl::OkStatus();


Note that our integrators are generally peppered with RETURN_IF_STOPPED so that we don't continue integrating if the user changed their mind. If you plan on integrating this in Principia, these statements will be needed.

Added some RETURN_IF_STOPPED.

pleroy · 2026-01-26T18:38:35Z

integrators/чебышёв_picard_iterator_body.hpp

+    Time const& step, ЧебышёвPicardIterator const& integrator)
+    : FixedStepSizeIntegrator<ODE>::Instance(problem, append_state, step),
+      integrator_(integrator),
+      t_(),


I'd omit this field, we generally don't list default initializations.

pleroy · 2026-01-26T18:41:37Z

integrators/чебышёв_picard_iterator_body.hpp

+      ᵝT(i, j) = 2 * τᵢ * ᵝT(i, j - 1) - ᵝT(i, j - 2);
+
+      // Make sure the zeroes are actually zero.
+      if (std::abs(ᵝT(i, j)) < 1e-14) ᵝT(i, j) = 0;


{} around the statement that follows the if.

More importantly, we are not going to like the magic number 1e-14 here. I don't understand the algorithm, so I'm not sure what to recommend, but these kinds of tests are disasters waiting to happen.

I removed the line; it wasn't actually necessary.

pleroy · 2026-01-26T18:45:06Z

serialization/integrators.proto

  optional CompositionMethod composition_method = 2;  // Only for SRKN.
+
+  // Only for CHEBYSHEV_PICARD.
+  optional ChebyshevPicardIterationParams chebyshev_picard_params = 3;


Since some of these params clearly belong to the instance, I'm not sure if this field belongs here.

I've removed the proto params (and the rest of the serialization code) for now; I'll add them back in a future PR.

(I don't want to add proto enums until I know what $M$ and $N$ are good)

… of N.

eggrobin · 2026-02-05T23:43:34Z

integrators/чебышёв_picard_iterator_body.hpp

+  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<std::remove_reference_t<decltype(yⱼ)>>;


It is always a bit sad when we have to drop the units like that. It happens, in flight plan optimization or in Newhall. Long ago all the integrators used to do that too; this is fine for now.

I wonder though, what would the units of those matrices be, and which matrix operations are used on them? The rows seem to be dependent variables multiplied by some scalar? ~~Are there any matrix multiplications or matrix-vector multiplications happening?~~ (EDIT: The answer to the last question is yes, e.g., ᵝT * ᵝW. But ᵝW is diagonal and dimensionless.)

I think all the matrices defined in the ЧебышёвPicardIterator constructor are dimensionless (as they don't depend on the ODE at all, only on $M$ and $N$).

The update operation is $X^{i+1} = C_xC_\alpha\cdot\frac12\cdot\mathrm{step}\cdot y'+C_xX_0$, where:

$y'$ is an $(M+1)\times n$ matrix where each row is a derivative at some time $t$ and each column has dimension of the corresponding element of DependentVariableDerivatives.

$C_xX_0$ is a $(M+1)\times n$ matrix where each row is a copy of the initial state $y(t_0)$ and each column has dimension of the corresponding element of DependentVariables.

So each column of $\mathrm{step}\cdot y'$ should have dimension as DependentVariableDifferences, multiplying this by the dimensionless $\frac12 C_xC_\alpha$ should not change this.

If there was a FixedMatrix variant which accepted different dimensions for each column this would probably work; in that case I wouldn't need to divide by Second in the update step.

…_using

eggrobin · 2026-02-06T00:15:45Z

I couldn’t think of a good way to test a fix to the GitHub Actions script, so I pushed changes to the check-iwyu check in here. It looks like it now manages to work well enough to tell you what changes in the log, though somehow the messages don’t show up in the diff (but I don’t know if that’s a bug on my side, GitHub having trouble with non-ASCII paths, or the fancy new GitHub UI not showing things where I expect them).

I’ll cherry-pick those changes into an other PR.

eggrobin · 2026-02-06T02:19:19Z

documentation/bibliography.bib

  volume       = {116},
 }

+@article{FeaginNarcozy1983,


It looks like the author is called Nacozy, no r?

Whoops! Fixed.

integrators/чебышёв_picard_iterator_body.hpp

eggrobin · 2026-02-06T03:44:02Z

integrators/чебышёв_picard_iterator_body.hpp

+
+      // Set the current state to the final state we appended.
+      current_state =
+          State(t_[Xⁱ_.rows() - 1],


It looks like everything is dropped down to single precision (well, that’s going to be double, but single in the sense of not being DoublePrecision) at each step.

This is typically disastrous for long-term integration with the kinds of methods we currently have; for many methods compensated summation of increments computed at each step is enough, for some (the multistep methods, see #1184) full double-double additions are needed. I don’t know how this one fares, but at least it looks like the independent variable is being straightforwardly incremented, so it will get its full dose of roundoff accumulation.

This is probably not immediately actionable though. The right way to figure this out would be to add this method to https://github.com/mockingbirdnest/Principia/blob/master/mathematica/integrator_plots.cpp; the before/after of #1184 was striking:

I am curious to see what the plots show. I'll add this method to integrator_plots.cpp after I implement the second-order version.

eggrobin · 2026-02-06T03:56:13Z

integrators/чебышёв_picard_iterator_test.cpp

@@ -0,0 +1,461 @@
+#include "integrators/чебышёв_picard_iterator.hpp"


I can’t immediately find a discussion of the convergence order in [Mac15]. This is usually something we test about our fixed-step integrators, though since this is driven by a stopping criterion, you would have to count actual evaluations rather than steps, like for the variable stepsize integrators. I don’t think we bothered doing that in the tests of the variable stepsize integrators.

The order (assuming that there is some fixed convergence order) would still show up as the slope of the error-work graph, see, e.g., #4194 (comment).

I haven't found any information on the convergence order either. I have found remarks that the iteration converges geometrically, but that isn't really relevant.

If we see a fixed convergence order in the error-work graph we can always add a test.

Co-authored-by: Robin Leroy <egg.robin.leroy@gmail.com>

rnlahaye added 24 commits January 12, 2026 14:19

MCPI: it's alive!

dfecc88

ChebyshevPicardIterationParams.

b95ba19

ChebyshevPicardIteratorTest TEST_P.

d3585e7

More test cases for ChebyshevPicardIteratorTest.

31bb19c

Move the solved IVP into the test parameters.

6143f45

Add the tangent problem from the NASA paper.

7103c72

Add the perturbed sinusoid problem.

e908836

MultipleSteps test.

5e1cf93

ChebyshevPicardIterator: backwards stepping.

f03fd4f

ChebyshevPicardIteratorTest divergence test.

5f9bb7b

Implement ChebyshevPicardIterator WriteToMessage.

2bc7369

Integrator #ifndef trick for ChebyshevPicardIterator.

c4429b7

Go back to absolute tolerances in CPI test.

72b52df

Use -1 as starting t for tangent test now that we are using absolute …

cb02b0e

…tolerance.

Rename variables to better match the source.

37e426b

Function for MaxNorm.

7895c44

Don't store Cₓ; it isn't really necessary.

43eef1b

Rename x to Xⁱ.

3c22631

Try not to allocate in the loops.

a73c0d2

Rename chebyshev to чебышёв in filenames.

0192266

Rename Chebyshev to Чебышёв inside the files.

b957136

Documentation and citations.

ac29708

Remove the norm printing and the decomposition include.

2849aa2

Move Fukushima1997 to the correct position.

6c5925a

rnlahaye added 3 commits January 16, 2026 07:21

Lint.

afe326d

Handle CHEBYSHEV_PICARD case in ReadFromMessage.

8c9c620

More lint.

1965e8d

pleroy reviewed Jan 26, 2026

View reviewed changes

rnlahaye added 9 commits January 26, 2026 19:31

Updates to bibliography as requested.

a346203

Address some comments.

1a66062

Address more comments.

6bbbaf2

L∞Norm -> LInfinity norm (fix Windows build).

16bcf39

Add ЧебышёвPicard method.

56a2a2e

Change more i++ to ++i (also j, k).

903079c

Parameterize ЧебышёвPicardIteratorParameterizedTest by method instead…

a11e5b5

… of N.

Fix computation of CₓCα_, and add a test for M > N.

0f235d7

Use FixedMatrix instead of UnboundedMatrix in ЧебышёвPicardIterator.

fed4c5f

eggrobin reviewed Feb 5, 2026

View reviewed changes

eggrobin added 2 commits February 6, 2026 01:00

Don’t escape filenames when listing files changed by include_what_you…

ec46317

…_using

[StringSplitOptions]::RemoveEmptyEntries

17be53f

eggrobin mentioned this pull request Feb 6, 2026

Don’t escape filenames when listing files changed by include_what_you_using #4419

Open

Try to make it work when it works

bfb2881

eggrobin reviewed Feb 6, 2026

View reviewed changes

rnlahaye and others added 4 commits February 5, 2026 21:48

IWYU.

d320f0e

Spell Nacozy's name correctly.

0139a50

Remove ЧебышёвPicard serialization for now.

abf6275

Use ‖ instead of ||.

697bf9a

Co-authored-by: Robin Leroy <egg.robin.leroy@gmail.com>

		@@ -0,0 +1,461 @@
		#include "integrators/чебышёв_picard_iterator.hpp"

Conversation

rnlahaye commented Jan 16, 2026

Regarding the implementation

Regarding future work

Uh oh!

pleroy commented Jan 16, 2026

Uh oh!

rnlahaye commented Jan 24, 2026

Uh oh!

pleroy commented Jan 24, 2026

Uh oh!

pleroy left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

eggrobin Feb 5, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

eggrobin commented Feb 6, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

eggrobin Feb 6, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

eggrobin Feb 5, 2026 •

edited

Loading

eggrobin Feb 6, 2026 •

edited

Loading