• At the moment we support DtO for both, of course.
• Reversible solvers only exist for Stratonovich so that's a moot point.
• The "ODE" OtD corresponds to Stratonovich.

The only one remaining is Ito OtD.

At the moment it's fixed to an LU-based solver, but it'd be nice to make this an option: in particular to make using an SVD-based solver possible, which is helpful for poorly-conditioned systems.

In particular Brownian motions, and possibly lazy diffeq solves.

Evaluating the control over [t_0, t_0.5] instead of over [t_0, t_1] raises the weak order on additive niose.

The latent ODE example seems to fail if using IController(rtol=1e-5, atol=1e-5, unvmap_dt=True). Possibly indicative of something going wrong with the adaptivity?

• Optimise-then-discretise
• Reversible solvers
• Interpolating adjoints
• Checkpointing?
• Quadrature?

This is relevant for CDEs in particular, but in principle one could also want to solve an ODE or SDE continuously as part of some larger program.

Related to #5?

• My thesis
• Relevant papers to each example

(c.f. https://diffeq.sciml.ai/stable/extras/timestepping/)

• PIController
• PIDController
• SDE-suitable controllers based on half-step-sizes
• Gustaffson acceleration / predictive controllers

At the moment all of our SDE solvers are Levy-area-free. It should be relatively straightforward to add support for different kinds of Levy area, by extending the evaluate interface.

At least as an option, this might be a nice thing to have.
Would need to think about the correct way of backpropagating through it + making sure XLA doesn't compile it away.

C.f. also Section VIII.5 of Hairer-Lubich-Wanner, which is the only reference I know of that actually discusses this.

When jit-compiling on the GPU: errors are currently only printed to stdout/stderr rather than properly raising. See google/jax#9457, including the bare-bones of a possible workaround. In practice it'd be better to format the traceback appropriately, using the traceback module; see also https://github.com/google/jax/blob/7a6986c4c8fd8469bae36306efb0417b0a2f6d8c/jax/_src/traceback_util.py for inspiration.

• Semi-implicit Euler
• Reversible Heun
• Leapfrog/midpoint

• In particular for vmap'ing.
• DETest

At the moment all global interpolation routines require that the sequence of times be strictly increasing. In practice it's often helpful to admit weakly increasing times as well, e.g. to pad variable-length time series.

In particular this is in contrast to dense interpolation (generated by diffeqsolve(saveat=SaveAt(dense=True))) which will sometimes use weakly increasing times. (Although this is a detail that is hidden from the user.)

Operations that need addressing:

• linear_interpolation
• rectilinear_interpolation
• backward_hermite_coeffs

(The corresponding classes should already handle this.)

The following produces very small mathematics in the documentation:

??? note

    $math$

It's a bit of a wart, used only for passing direction.

Likewise for AbstractStepSizeController.wrap_solver.

Available at: https://github.com/google/tree-math

Still need to see if this is implements every required operation. (In particular around implicit solvers.)

https://arxiv.org/abs/2009.09457

One thing that needs thinking about is how this compares to the use of jump processes as a driving control, which morally speaking do something very similar.

• Discrete terminating events
• Discrete non-terminating events
• Continuous terminating events
• Continuous non-terminating events

The continuous events can be implemented by using the most recent dense_info and performing a nonlinear solve to locate the event location.

SDE:

• Milstein
• Euler-Heun
• SRKs (Additive)
• Talay (Ito commutative)

Implict:

• BDF
• Rosenbrock methods
• IMEX methods, e.g. http://dx.doi.org/10.1175/mwr-d-13-00132.1

DAE:

• Lots of overlap with Rosenbrock methods, get them in at the same time.
• Standard DAE->ODE conversion (+optional nonlinear projections)
• Set made_jump=True if dconstraint/dz=0
• Mass matrices: work with these or work with constraint functions? Former is a little more general; maybe less efficient in some cases? Also think about traced mass matrices, possibly space-varying (e.g. quasilinear problems), possibly rank-varying.

Symplectic:

• Higher-order symplectic methods
• A lot of the lower-order methods (Stormer-Verlet, Leapfrog, Kick-drift-kick, Semi-implicit Euler) are all essentially the same thing. It would be good to write out something making this explicit, plus maybe add some special variants if you really want a particular variant.

Other:

• Higher-order ERK methods (Verner etc.)
• See also a few of the comments made about DifferentialEquations.jl in this (old) blog post: https://www.stochasticlifestyle.com/comparison-differential-equation-solver-suites-matlab-r-julia-python-c-fortran/

We now have a couple of simple ODE benchmarks, but could definitely afford to add a lot more.

Currently we're using local quadratics. It should maybe be local semi-ovals?

Hi,
first of all, let me say that this looks like an amazing project.
I am looking forward to playing around with this :).

In a concrete problem I am dealing with, I have a forced system where the external force is piecewise constant. The external force changes at specific time points (t1, ..., tn), causing a discontinuity of the time derivative.
I would like to use adaptive step-size solvers for increased accuracy, but naively applying adaptive step-size solvers will "waste" a lot of steps to find the point of change.

Would including the change points in SaveAt avoid this problem?
Or is there some other recommended way to handle this?

We have a lot of BibTeX references specified throughout the documentation. It would be nice to add:

• DOI references to the BibTeX itself

@article{blahblah,
    ...
    doi={...}
}

• arXiv/DOI hyperlinks adjacent to the BibTeX

[arXiv](...)
@article{...
}

• Brownian Tree
• Piecewise linear
• Karhunen--Loeve

At the moment it is just the same as a jax.jvp of solution.evaluate; in particular it is the derivative of the numerical spline approximation; not an approximation to the true derivative. The latter might be more helpful? Not sure.

• Linear interpolation
• Rectilinear interpolation
• Hermite cubic splines with backward differences

At the moment:

• the ODE order is specified
• an SDE order is specified -- in principle for whatever the most general type of noise that solver is expecting. (e.g. general noise for Euler, commutative noise for Milstein)

This isn't very flexible. For example what about solving an additive-noise problem with Heun? In this case Heun gets strong order 1, rather than the 0.5 it is specified with. At the moment orders are (only) used for adaptive stepping, so practically speaking this can be handled by passing the appropriate local_order to diffrax.PIDController, but this isn't well-advertised.

Generally speaking determining this automatically seems to be a huge can of worms. The strong convergence order for every kind of solver for every kind of noise simply aren't known. (Which is the reason we punt the problem into user-land instead.) Moreover we'd probably need to do some pretty involved introspection of the information we're passed to determine what kind of order to expect, that may well simply be wrong in edge cases.

At the moment there are several stop_gradients in there, primarily to avoid instability observed in the backward pass. There's been a few bugs fixed in IController since then, so this may no longer be the case?

It may still be nice to keep as an option, since non-backpropagation through rejected steps is an increase in efficiency with minimal change to the gradient.

During the early phases of the development of this library, we pretty much copied the interpolation routines wholesale from torchdiffeq without much thought. In particular this means a 4th order polynomial interpolation scheme is used for several solvers -- which Dopri5 uses with the appropriate care, but most other algorithms do not.

Where no better choice is available then these other algorithms should switch to use 3rd order Hermite interpolation, which is really the standard choice when a custom dense interpolation scheme hasn't been developed.

Hypersolvers are learnt solvers, as per https://arxiv.org/abs/2007.09601.
See also their implementation in https://github.com/DiffEqML/torchdyn

This is a nice-to-have stretch goal: a "lazy" Solution when using saveat.dense that only computes the solution at query time, rather than in advance.

That is, specify step locations and jump locations for the solver.

• Initial automated stepsize selection is giving different values to torchdiffeq (and needs using on the latent ODE example)
• _step is getting recompiled 61 steps into a CNF?

Very cool work!

It would be great to have an example on how to solve an ODE that is directly forced by some signal x(t), e.g. a forced mass-spring-damper

m y'' + r y' + k y = x(t).

Do I understand correctly, that the controlled ODEs are "forced" by the derivative of x(t)?

• docstrings etc.
• The full extent of things we can vmap -- vmap all the things!
• Make sure to include solver.step as an explicit part of the public API, for use independent of diffeqsolve.

At the moment these custom operations aren't yet part of the public API, but they should make it possible to replace unvmap_* and nondifferentiable_output, which are currently implemented as custom primitives.

Leave a comment here: google/jax#2628 once Diffrax goes public.