We introduce the package ManifoldDiffEq.jl to solve differential equations on Riemannian manifolds. The package implements algorithms using interfaces of ManifoldsBase.jl and OrdinaryDiffEq.jl, combining strengths of both of these worlds.
Currently it offers solvers for two kinds of differential equations. The first one is of the form
\frac{\mathrm{d}}{\mathrm{dt}}u = f(u,p,t)
on a time interval t \in [a,b], where p are some parameters and f(u,p,t) is a tangent vector to the point u(t) on a Riemannian manifold \mathcal M.
The second one is similar but the tangent vector is indirectly specified by an action of a Lie group on \mathcal M.
Algorithms in this package can be used with all manifolds from Manifolds.jl, and utilize the exponential map (or its approximations) and vector transports.
For algorithms based on Lie groups, also the group decorator for Riemannian manifolds from Manifolds.jl is used, and any kind of Lie group action can be easily applied.
Additionally, the solution is interpolated using manifold interpolation methods.
Currently the following algorithms are available
- A Forward Euler solver both for manifolds as well as one using the Lie group structure
- Crouch-Grossmann (CG) 2 and 3 for the frozen coefficient approach
- Runge-Kutta Munthe-Kaas (RKMK) 4
Kind regards
Mateusz (@mateuszbaran) and Ronny (@kellertuer)
17 Likes
Cool package! It looks like it only implements a subset of DiffEq functionality though (less solvers and no timestep adaptation?), is that by design? Can you explain what it does better than just using a DiffEq callback to project back to the manifold after each step?
Well for not it implements (only) 5 solvers since we just started. We plan to do more Lie group integrators for sure, but some of them are quite involved. The same holds for intrinsic manifold integrators.
That is also the main difference. For Lie groups the differential equation can be tackled in the Lie algebra and you never leave the algebra or the group, so you do not need a projection (or an embedding you work in such that you can/have to project).
Similarly for example the forward Euler on a manifold uses a retraction (or if available you can use the exponential map). One possible retraction for embedded manifolds that have projections available, is the projection for sure.
But for SPD matrices for example /which are an open set in the set of matrices) there is no projection.
Also the projection, while it is sometimes available and then it is a retraction, might not be the best choice (computationally or in the exactness/precision you aim for).
Using the intrinsic methods without embed/project might – depending on the manifold at hand – also reduce the memory you need since you can stay manifold_dimension(M) sized instead of having to use representations in the embedding. That of course is not always the case (The Sphere is easier if you use unit vectors for example), but it might be an advantage, too.
Time steps and such are planned as well, for sure – and they might even be easily available from DiffEq, I have not et checked. When possible, we for sure want to use tools from DiffEq. But to handle the Manifold, we also always use the tools from Manifolds.jl.
Edit: Oh for Lie groups there is even another challenge. Since tangent vectors (and hence the differential equation) are phrased in the Lie algebra (or the tangent space of the identity) before projecting you would have to (parallel or vector) transport the Lie algebra elements to the right place, which is an operation our approach can avoid there.
4 Likes
Yes, we’re definitely going to add more integrators, and include methods that support timestep adaptation.
The main things I don’t know how to do in plain DiffEq is support for arbitrary exponential maps and vector transports for Crouch-Grossman methods, arbitrary Lie group actions in RKMK methods, and on-manifold interpolation of the solution.
Projection is definitely one way to solve ODEs on manifolds. Maybe these methods would be enough for me but for CG and RKMK methods looked much more suitable. Do multi-stage time-reversible projection methods actually work better? The number of projections that needs to happen would definitely be challenging to track in proofs of order of convergence.
4 Likes