CTFlows.jl

Appearance

Flow(Hamiltonian) compatibility ​

This page is a living compatibility reference for the Hamiltonian flow built from a scalar Data.Hamiltonian function: which state/costate types and call styles it accepts, each shown with a minimal, executable example. Every ✓ / ⚠ / ✗ in the table below is demonstrated by a code block on this page and is re-run on every documentation build, so the page cannot drift from the code.

Scope: CPU, default AD backend (AutoForwardDiff). The table is generated from probe/cpu/probe_cpu.jl (run it locally with julia --project=probe/cpu probe/cpu/probe_cpu.jl). GPU compatibility (method=:gpu, default backend AutoMooncake) is a separate effort — see GPU flows and probe/gpu.

Unlike Flow(VectorField) and Flow(HamiltonianVectorField), a scalar Hamiltonian has no in-place variant — there is no mutability trait on Data.Hamiltonian at all. The flow computes  ,    by automatic differentiation instead of from a user-supplied vector field, so the axis that varies on this page is not OOP/IP but the input type accepted by the AD backend.

All examples use the harmonic oscillator   (autonomous, fixed) — the same dynamics as Flow(HamiltonianVectorField), so the analytic solution is identical:    at  .

Last probed: 2026-07-23

This page's table shape — which containers and axes are tested — reflects probe/cpu/probe_cpu.jl as of this date. Every ✓/✗ cell below is still re-executed on every documentation build regardless (see Compatibility overview) — only the scope of what's tested can go stale, not the results shown.

The flow is built once — no ad_backend keyword needed, since AutoForwardDiff is already the CPU default:

julia
h     = Data.Hamiltonian((x, p) -> 0.5 * (sum(abs2, x) + sum(abs2, p)))
hflow = Flows.Flow(h; reltol=1e-8)
Flow
├─ system: HamiltonianSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  ├─ Hamiltonian: autonomous, fixed (no variable)
│  │    natural call: h(x, p)
│  │    uniform call: h(t, x, p, v)
│  └─ backend: DifferentiationInterface{CPU}(ad_backend=AutoForwardDiff)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-8  [default]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.

An explicit backend, or a GPU one, can be selected the same way: Flows.Flow(h; ad_backend=ADTypes.AutoForwardDiff()) — see GPU flows for method=:gpu and AutoMooncake/AutoZygote.

Compatibility table ​

The two call styles are the point call hflow(t0, x0, p0, tf) → final (xf, pf), and the trajectory call hflow((t0, tf), x0, p0) → Trajectories.HamiltonianVectorFieldTrajectory. There is no OOP/IP split.

State/costate typepointtraj
Scalar Real✓✓
Vector Real✓✓
MVector Real✓✓
SVector Real✓✓
Matrix Real (batch)✓✓
MMatrix Real✓✓
SMatrix Real✓✓
Complex (any container)✗ (a)✗ (a)

Legend ​

  • ✓ — works; the result matches the analytic    and is executed on this page.

  • ✗ — raises an error; see the corresponding note.

(a) Complex is not supported by the default backend

AutoForwardDiff cannot differentiate through a Complex input: constructing a dual number over ComplexF64 raises ArgumentError: Cannot create a dual over scalar type ComplexF64. This is a ForwardDiff.jl limitation, not specific to CTFlows.

If your Hamiltonian naturally involves complex state (e.g. from a Schrödinger-type equation), it can typically be rewritten as a real system by splitting the state into real/imaginary (or conjugate) components — this is not a limitation of CTFlows, just a modeling choice that avoids differentiating through Complex directly.

The table above is about the flow's state/costate input. A separate question — how to differentiate the flow itself with respect to (x0, p0), e.g. for a shooting method — is covered next.

Real states ​

Scalar ​

julia
julia> hflow(0.0, 1.0, 0.0, pi/2)   # (xf, pf) ≈ (0.0, -1.0)
(-1.5541996200011742e-10, -0.9999999997009266)
julia
julia> sol = hflow((0.0, pi/2), 1.0, 0.0);

julia> (Trajectories.state(sol)(pi/2), Trajectories.costate(sol)(pi/2))   # both scalars
(-1.5541963451763603e-10, -0.9999999997009266)

Vector ​

julia
julia> hflow(0.0, [1.0, 0.0], [0.0, 1.0], pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

MVector and SVector ​

julia
julia> hflow(0.0, MVector{2}(1.0, 0.0), MVector{2}(0.0, 1.0), pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

julia> hflow(0.0, SA[1.0, 0.0], SA[0.0, 1.0], pi/2)
([-1.5541996200011742e-10, 0.9999999997009266], [-0.9999999997009266, -1.5541996200011742e-10])

Matrix (batch) ​

H is separable and quadratic (sum(abs2, x)), so it reduces a matrix state/costate to a scalar the same way it does a vector, and the induced dynamics are independent per column:

julia
julia> hflow(0.0, [1.0 2.0; 3.0 4.0], [0.0 0.0; 1.0 1.0], pi/2)
([-9.327931901967881e-11 -1.8655863803935762e-10; 0.9999999995245173 0.9999999994312379], [-0.9999999998043553 -1.9999999996087106; -2.9999999995063455 -3.9999999993107])

julia> hflow(0.0, MMatrix{2,2}(1.0, 3.0, 2.0, 4.0), MMatrix{2,2}(0.0, 1.0, 0.0, 1.0), pi/2)   # mutable
([-9.327931901967881e-11 -1.8655863803935762e-10; 0.9999999995245173 0.9999999994312379], [-0.9999999998043553 -1.9999999996087106; -2.9999999995063455 -3.9999999993107])

julia> hflow(0.0, SMatrix{2,2}(1.0, 3.0, 2.0, 4.0), SMatrix{2,2}(0.0, 1.0, 0.0, 1.0), pi/2)   # immutable
([-9.327931901967881e-11 -1.8655863803935762e-10; 0.9999999995245173 0.9999999994312379], [-0.9999999998043553 -1.9999999996087106; -2.9999999995063455 -3.9999999993107])

Complex states — not supported ​

julia
julia> try
           hflow(0.0, 1.0 + 2.0im, 0.0 + 0.0im, pi/2)
       catch e
           showerror(stdout, e)
       end
ArgumentError: Cannot create a dual over scalar type ComplexF64. If the type behaves as a scalar, define ForwardDiff.can_dual(::Type{ComplexF64}) = true.

See note (a) above — for a system that is naturally complex, rewrite it as a real system on the real/imaginary (or conjugate) components instead.

Automatic differentiation: sensitivities of the flow ​

Since Flow(Hamiltonian) is itself AutoForwardDiff-backed internally, differentiating it with respect to (x0, p0) — e.g. for a shooting method — means nested AD. Do this by differentiating the flow call, never by constructing a ForwardDiff.Dual by hand and passing it as x0/p0: a hand-built, untagged Dual collides with the flow's own internal tag (DualMismatchError), and a hand-supplied custom tag does not fix it either — an arbitrary marker type lacks ForwardDiff.Tag's internal ordering machinery. Wrapping the call in an outer ForwardDiff.jacobian / ForwardDiff.gradient sidesteps this entirely, because ForwardDiff then generates its own properly ordered tag for the perturbation — this is exactly how a NonlinearSolve-style shooting method differentiates a flow that is itself ForwardDiff-backed.

This computes the exact Jacobian of with respect to — for the linear harmonic oscillator it is the constant matrix :

julia
julia> shoot(z) = collect(hflow(0.0, z[1], z[2], pi/2))
shoot (generic function with 1 method)

julia> ForwardDiff.jacobian(shoot, [1.0, 0.0])
2×2 Matrix{Float64}:
 -1.5542e-10   1.0
 -1.0         -1.5542e-10

The same closure-wrapping pattern works with ForwardDiff.gradient for a scalar shooting residual. No internal backend change (e.g. switching to AutoZygote/AutoMooncake to dodge the nested Dual) is needed for this to work.

See also ​