Manifold Projection

ManifoldProjection(
    manifold; nlsolve = missing, save = true, autonomous = nothing,
    manifold_jacobian = nothing, autodiff = nothing,
    resid_prototype = nothing, kwargs...) -> DiscreteCallback

In many cases, you may want to declare a manifold on which a solution lives. Mathematically, a manifold M is defined by a function g as the set of points where g(u) = 0. An embedded manifold can be a lower dimensional object which constrains the solution. For example, g(u) = E(u) - C where E is the energy of the system in state u, meaning that the energy must be constant (energy preservation). Thus by defining the manifold the solution should live on, you can retain desired properties of the solution.

ManifoldProjection projects the solution of the differential equation to the chosen manifold g, conserving a property while conserving the order. It is a consequence of convergence proofs both in the deterministic and stochastic cases that post-step projection to manifolds keep the same convergence rate, thus any algorithm can be easily extended to conserve properties. If the solution is supposed to live on a specific manifold or conserve such property, this guarantees the conservation law without modifying the convergence properties.

Arguments

• manifold: residual function defining the manifold. For an in-place problem, define manifold(resid, u, p) or manifold(resid, u, p, t). For an out-of-place problem, define manifold(u, p) or manifold(u, p, t). The residual is zero on the manifold.

Keywords

• nlsolve = missing: use the built-in single-factorization projection algorithm. Pass a nonlinear solver in the NonlinearSolve.jl format to select another algorithm. Pass nothing to use the NonlinearSolve polyalgorithm.
• save::Bool = true: save immediately after the projection.
• autonomous = nothing: whether manifold omits the time argument. Pass Val(true) or Val(false) to avoid runtime branching; nothing infers the form during initialization.
• resid_prototype = nothing: prototype defining the residual shape for an in-place problem. Without one, the residual is assumed to have the same shape as u.
• autodiff = nothing: DifferentiationInterface automatic-differentiation backend used when manifold_jacobian is not supplied.
• manifold_jacobian = nothing: analytic Jacobian of manifold with respect to u. Use the same calling form as manifold, with the Jacobian output as the first argument for an in-place problem.
• kwargs...: additional keywords passed to the built-in projection algorithm or to NonlinearSolve.jl if nlsolve is not missing.

Returns

• DiscreteCallback: a callback that projects the state after each accepted step. If the nonlinear projection does not converge, the callback terminates the integrator with the projection solver's unsuccessful return code.

Throws

• ErrorException: during callback initialization if both manifold_jacobian and autodiff are nothing.

Saveat Warning

Note that the ManifoldProjection callback modifies the endpoints of the integration intervals and thus breaks assumptions of internal interpolations. Because of this, the values for given by saveat will not be order-matching. However, the interpolation error can be proportional to the change by the projection, so if the projection is making small changes then one is still safe. However, if there are large changes from each projection, you should consider only saving at stopping/projection times. To do this, set tstops to the same values as saveat. There is a performance hit by doing so because now the integrator is forced to stop at every saving point, but this is guaranteed to match the order of the integrator even with the ManifoldProjection.

References

[1] Ernst Hairer, Christian Lubich, Gerhard Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Berlin ; New York :Springer, 2002.

Examples

using ADTypes, DiffEqCallbacks, OrdinaryDiffEq

function unit_circle(resid, u, p, t)
    resid[1] = sum(abs2, u) - 1
end

function rotation!(du, u, p, t)
    du[1] = -u[2]
    du[2] = u[1]
end

prob = ODEProblem(rotation!, [1.0, 0.0], (0.0, 10.0))
cb = ManifoldProjection(unit_circle; resid_prototype = [0.0], autodiff = AutoFiniteDiff())

sol = solve(prob, Tsit5(); callback = cb)