Manifold Projection

The following callbacks are designed to provide post-step modifications to preserve geometric behaviors in the solution.

DiffEqCallbacks.ManifoldProjectionType
ManifoldProjection(g;nlsolve=NLSOLVEJL_SETUP(),save=true)

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

  • g: The residual function for the manifold. This is an inplace function of form g(resid, u) or g(resid, u, p, t) which writes to the residual resid the difference from the manifold components. Here, it is assumed that resid is of the same shape as u.

Keyword Arguments

  • nlsolve: A nonlinear solver as defined in the nlsolve format
  • save: Whether to do the standard saving (applied after the callback)
  • autonomous: Whether g is an autonomous function of the form g(resid, u).
  • nlopts: Optional arguments to nonlinear solver which can be any of the NLsolve keywords.

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 guerenteed to match the order of the integrator even with the ManifoldProjection.

References

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

Example

Here we solve the harmonic oscillator:

using OrdinaryDiffEq, DiffEqCallbacks, Plots

u0 = ones(2)
function f(du,u,p,t)
  du[1] = u[2]
  du[2] = -u[1]
end
prob = ODEProblem(f,u0,(0.0,100.0))
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
timespan: (0.0, 100.0)
u0: 2-element Vector{Float64}:
 1.0
 1.0

However, this problem is supposed to conserve energy, and thus we define our manifold to conserve the sum of squares:

function g(resid,u,p,t)
  resid[1] = u[2]^2 + u[1]^2 - 2
  resid[2] = 0
end
g (generic function with 1 method)

To build the callback, we just call

cb = ManifoldProjection(g)
DiscreteCallback{DiffEqCallbacks.var"#7#8", ManifoldProjection{false, DiffEqCallbacks.NonAutonomousFunction{typeof(Main.g), false}, DiffEqCallbacks.NLSOLVEJL_SETUP{0, true}, Dict{Symbol, Any}}, typeof(DiffEqCallbacks.Manifold_initialize), typeof(SciMLBase.FINALIZE_DEFAULT)}(DiffEqCallbacks.var"#7#8"(), ManifoldProjection{false, DiffEqCallbacks.NonAutonomousFunction{typeof(Main.g), false}, DiffEqCallbacks.NLSOLVEJL_SETUP{0, true}, Dict{Symbol, Any}}(DiffEqCallbacks.NonAutonomousFunction{typeof(Main.g), false}(Main.g, 0, 0), DiffEqCallbacks.NonAutonomousFunction{typeof(Main.g), false}(Main.g, 0, 0), DiffEqCallbacks.NLSOLVEJL_SETUP{0, true}(), Dict{Symbol, Any}()), DiffEqCallbacks.Manifold_initialize, SciMLBase.FINALIZE_DEFAULT, Bool[0, 1])

Using this callback, the Runge-Kutta method Vern7 conserves energy. Note that the standard saving occurs after the step and before the callback, and thus we set save_everystep=false to turn off all standard saving and let the callback save after the projection is applied.

sol = solve(prob,Vern7(),save_everystep=false,callback=cb)
@show sol[end][1]^2 + sol[end][2]^2 ≈ 2
true
using Plots
plot(sol,idxs=(1,2))