Integrating

Once a flow is built, calling it integrates the underlying ODE. There are two call styles depending on whether you need the full trajectory or just the final state.

---

Call styles

Point integration — final state only

flow(t0, x0, tf)            # returns xf::Vector (StateFlow)
hflow(t0, x0, p0, tf)       # returns (xf, pf) (HamiltonianFlow)

x0 = [1.0, 0.0]
xf = flow(0.0, x0, 1.0)

2-element Vector{Float64}:
 0.36787944127643135
 0.0

p0 = [0.0, 1.0]
xf, pf = hflow(0.0, x0, p0, 1.0)
(xf, pf)

([0.540302305784261, 0.8414709847533872], [-0.8414709847533872, 0.540302305784261])

Trajectory integration — full time history

flow((t0, tf), x0)          # returns VectorFieldTrajectory
hflow((t0, tf), x0, p0)     # returns HamiltonianVectorFieldTrajectory

sol = flow((0.0, 1.0), x0)

VectorFieldTrajectory
├─ result: SciMLIntegrationResult
├─ tspan: (0.0, 1.0)
├─ time points: 250
└─ final state: [0.36787944127643135, 0.0]

hsol = hflow((0.0, 1.0), x0, p0)

HamiltonianVectorFieldTrajectory
├─ result: SciMLIntegrationResult
├─ tspan: (0.0, 1.0)
├─ time points: 250
├─ final state: [0.540302305784261, 0.8414709847533872]
└─ final costate: [-0.8414709847533872, 0.540302305784261]

Output grid and accuracy

A trajectory is returned on a time grid (Integrators.times(sol)) but can be evaluated at any time: sol(t) reads the solver's dense interpolant, accurate to the solver tolerances. When the call sets no output grid, the trajectory is returned on a default grid: if the solver returns fewer than 250 times (an easy integration, hence a coarse polyline when plotted), an AdaptiveGrid(250) is generated, so that the result has 250 distinct times — the grid_size of the solutions of the direct methods, which makes a Solution built by a flow interchangeable with a direct one. With 250 times or more, the result is left as the solver returned it.

length(Integrators.times(flow((0.0, 1.0), x0)))

250

The default grid applies to a top-level trajectory call on a time span (t0, tf) with a dense integration. It is not generated with saveat, a grid of times, grid=, dense=false, a point call f(t0, x0, tf), or a SciMLProblemFlow (which keeps the raw SciML behaviour). To keep the solver steps, pass grid=nothing:

length(Integrators.times(flow((0.0, 1.0), x0; grid=nothing)))

16

To get the trajectory on other or chosen times, use one of:

• grid= at the call, to generate a grid of n distinct times from the span — an integer n (uniform), Configs.UniformGrid(n), or Configs.AdaptiveGrid(n), denser where the plotted curves (state, costate, control) bend:

Integrators.times(flow((0.0, 1.0), x0; grid=11))

11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0

Integrators.times(flow((0.0, 1.0), x0; grid=Configs.AdaptiveGrid(11)))

11-element Vector{Float64}:
 0.0
 0.08193333476902397
 0.16706285327005757
 0.2556324205827278
 0.3479239622845741
 0.4442360641220789
 0.5449232095312161
 0.6503683236433677
 0.7610253531149139
 0.8773785987374856
 1.0

• a grid of times instead of the time span — a vector, a range, or a tuple of at least three times — returned exactly:

sol_grid = flow(range(0.0, 1.0, 11), x0)
Integrators.times(sol_grid)

11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0

• saveat when building the flow (a vector, a range, or a step):

flow_saveat = Flows.Flow(vf; reltol=1e-8, saveat=0.1)
Integrators.times(flow_saveat((0.0, 1.0), x0))

11-element Vector{Float64}:
 0.0
 0.1
 0.2
 0.3
 0.4
 0.5
 0.6
 0.7
 0.8
 0.9
 1.0

Combining two of them (saveat and a grid at the call, or a grid of times and grid=) throws an IncorrectArgument. In every case the grid only shapes the output: the integration (steps, cost) and the dense interpolant are exactly those of flow((t0, tf), x0), so sol(t) between the grid points, the objective and the control of an optimal control flow keep the solver accuracy. Multi-phase flows behave the same way: each phase keeps its own interpolant, and a generated grid (the default one included) counts every switching time once.

AdaptiveGrid(n; uniform=0.1) equidistributes the density ‖y''‖^{1/2} of the plotted curves (each scaled by its range), which minimizes the error of the polyline drawn through the grid, mixed with a uniform share uniform.

Pass dense=false to store only the grid values and save memory: the trajectory is then linear between the saved points, and an optimal control flow warns that its objective is only accurate to the grid spacing.

A point call flow(t0, x0, tf) is not affected by any of these: it ignores saveat and returns the value at tf, whether tf is on the grid or not.

---

Variable parameters

For a NonFixed flow, pass the variable via the variable keyword:

vf_v = Data.VectorField((x, v) -> -v[1] .* x; is_variable=true)

VectorField: autonomous, variable, out-of-place
  natural call: f(x, v)
  uniform call: f(t, x, v)

flow_v = Flows.Flow(vf_v)

Flow
├─ system: VectorFieldSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: NonFixed
│  └─ VectorField: autonomous, variable, out-of-place
│       natural call: f(x, v)
│       uniform call: f(t, x, v)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-8  [default]
   ├─ 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.

xf_v = flow_v(0.0, [1.0, 0.0], 1.0; variable=[2.0])

2-element Vector{Float64}:
 0.13533528348480478
 0.0

The variable argument is required when Traits.is_variable(flow) is true, and rejected for a Fixed flow: both mistakes throw a CTBase.Exceptions.PreconditionError.

---

Variable costate

For a NonFixed HamiltonianFlow, pass variable_costate=true to also integrate the augmented adjoint    (initialized at  ) alongside the state and costate. The point call then returns a triple (xf, pf, pvf) instead of (xf, pf):

h_v = Data.Hamiltonian((x, p, v) -> v[1] * p^2 / 2; is_variable=true)

Hamiltonian: autonomous, variable
  natural call: h(x, p, v)
  uniform call: h(t, x, p, v)

hflow_v = Flows.Flow(h_v)

Flow
├─ system: HamiltonianSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: NonFixed
│  ├─ Hamiltonian: autonomous, variable
│  │    natural call: h(x, p, v)
│  │    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  [default]
   ├─ 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.

julia> xf, pf, pvf = hflow_v(0.0, 1.0, 0.5, 1.0; variable=[2.0], variable_costate=true);

julia> xf
1.9999999999999998

julia> pf
0.5

julia> pvf
-0.12499999999999997

This is only available for point evaluation, and only when variable is provided (NonFixed flows require it).

Free times

The positional tf argument is the evaluation time; it is independent of the variable value, even when a variable component represents a free initial or final time. Passing flow(t0, x0, p0, t1; variable=v) with t1 ≠ v is valid.

When a variable component is a free time, the flow keeps integrating the same naive adjoint — no special-casing — with the augmented costate started at   (the convention variable_costate=true always uses). This zero start is exactly what makes the mitigated free-time transversality conditions valid: written at , they are

where H is obtained from CTFlows.Systems.hamiltonian(flow); a nonzero would shift both sides. The shooting method writes these conditions by hand.

See test/suite/flows/test_variable_costate_free_time.jl for worked shooting residuals built on this mechanism (free , free , and both at once). The Goddard tests (test/suite/integration/test_goddard_ocp.jl) instead close their free final time with the classical   condition, not the adjoint.

---

Hamiltonian / pseudo-Hamiltonian getters

Any Hamiltonian flow exposes its underlying Hamiltonian and — when the flow was built from a control law — its pseudo-Hamiltonian, control law, and their gradients:

Systems.hamiltonian(hflow_v)          # the callable H(t, x, p, v)

Hamiltonian: autonomous, variable
  natural call: h(x, p, v)
  uniform call: h(t, x, p, v)

# functor: (t, x, p, v) -> (∂H/∂x, ∂H/∂p)
Systems.get_hamiltonian_gradient(hflow_v)

CTFlows.Systems.HamiltonianGradient{CTBase.Data.Hamiltonian{Main.var"#11#12", CTBase.Traits.Autonomous, CTBase.Traits.NonFixed}, CTBase.Differentiation.DifferentiationInterface{CTBase.Strategies.CPU, CTBase.Strategies.StrategyOptions{@NamedTuple{ad_backend::CTBase.Options.OptionValue{AutoForwardDiff{nothing, Nothing}}}}}}(Hamiltonian: autonomous, variable
  natural call: h(x, p, v)
  uniform call: h(t, x, p, v), DifferentiationInterface{CPU}(ad_backend=AutoForwardDiff))

Systems.get_variable_gradient(hflow_v)    # functor: (t, x, p, v) -> ∂H/∂v

CTFlows.Systems.HamiltonianVariableGradient{CTBase.Data.Hamiltonian{Main.var"#11#12", CTBase.Traits.Autonomous, CTBase.Traits.NonFixed}, CTBase.Differentiation.DifferentiationInterface{CTBase.Strategies.CPU, CTBase.Strategies.StrategyOptions{@NamedTuple{ad_backend::CTBase.Options.OptionValue{AutoForwardDiff{nothing, Nothing}}}}}}(Hamiltonian: autonomous, variable
  natural call: h(x, p, v)
  uniform call: h(t, x, p, v), DifferentiationInterface{CPU}(ad_backend=AutoForwardDiff))

Systems.pseudo_hamiltonian, Systems.control_law, Systems.get_pseudo_hamiltonian_gradient and Systems.get_pseudo_variable_gradient are available on flows built from a pseudo-Hamiltonian (or an OCP) and a control law — see Control laws. Like the Hamiltonian, the gradient functors also take the short form allowed by the flow's time and variable dependence: (x, p) for an autonomous flow with no variable, (t, x, p) or (x, p, v) otherwise ((x, p, u), ... for the pseudo-Hamiltonian ones).

An accessor that does not apply to a flow throws IncorrectArgument, saying what the flow holds and which accessor to use instead: a state flow has no Hamiltonian (use Systems.vector_field), a flow built from a Hamiltonian vector field has no scalar Hamiltonian (use Systems.hamiltonian_vector_field), and a flow with no control law has no pseudo-Hamiltonian. Systems.control_law also returns the law of a state flow built from a control law, e.g. Flow(ocp, OpenLoop(u)).

---

Configuration objects

The convenience call signatures above internally build configuration objects that bundle the integration parameters. You can also construct them explicitly and pass them to Flows._invoke_flow:

Config type	Usage	Arguments
StateEndPointConfig	state final value	(t0, x0, tf)
StateTrajectoryConfig	state full trajectory	(tspan, x0)
HamiltonianEndPointConfig	state+costate final value	(t0, x0, p0, tf)
HamiltonianTrajectoryConfig	state+costate trajectory	(tspan, x0, p0)
AugmentedHamiltonianEndPointConfig	state+costate+variable-costate final value	(t0, x0, p0, pv0, tf)

cfg = Configs.StateTrajectoryConfig((0.0, 1.0), [1.0, 0.0])

StateTrajectoryConfig
├─ tspan: (0.0, 1.0)
└─ x0: [1.0, 0.0]

Configs.tspan(cfg)

(0.0, 1.0)

Configuration objects separate what to integrate from how to integrate (the flow). This separation is useful when the same config must be passed to several flows.

---

Integrator options

Options are passed as keyword arguments to Flows.Flow(data; opts...) or, when building the integrator explicitly, to Integrators.SciML(; opts...).

The default integrator is SciML backed by OrdinaryDiffEqTsit5 (loaded when using OrdinaryDiffEqTsit5: OrdinaryDiffEqTsit5 appears in your session).

Common options

Option	Default	Description
reltol	1e-8	Relative tolerance
abstol	1e-8	Absolute tolerance
alg	Tsit5()	ODE algorithm (any SciML algorithm)
saveat	—	Output time grid (vector, range or step); does not change the integration; not combinable with a grid at the call — see Output grid and accuracy
dense	:auto	Keep the dense interpolant (true for a trajectory); false saves memory, linear between saved points

# Tighter tolerances
flow_tight = Flows.Flow(vf; reltol=1e-12, abstol=1e-12)

Flow
├─ system: VectorFieldSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  └─ VectorField: autonomous, fixed (no variable), out-of-place
│       natural call: f(x)
│       uniform call: f(t, x, v)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Tsit5  [default]
   ├─ reltol = 1.0e-12  [user]
   ├─ save_everystep = auto  [default]
   ├─ abstol = 1.0e-12  [user]
   ├─ save_start = auto  [default]
   └─ dense = auto  [default]
   Tip: use describe(SciML{CPU}) to see all available options.

Different algorithm (requires the matching OrdinaryDiffEq package to be loaded).

using OrdinaryDiffEqRosenbrock
flow_rodas = Flows.Flow(vf; alg=Rodas4())

Flow
├─ system: VectorFieldSystem
│  ├─ time_dependence: Autonomous
│  ├─ variable_dependence: Fixed
│  └─ VectorField: autonomous, fixed (no variable), out-of-place
│       natural call: f(x)
│       uniform call: f(t, x, v)
└─ integrator: SciML{CPU} (instance, id=:sciml)
   ├─ internalnorm = real_norm  [default]
   ├─ alg = Rodas4  [user]
   ├─ reltol = 1.0e-8  [default]
   ├─ 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.

Unsafe mode

By default, a CTBase.Exceptions.SolverFailure is thrown if the ODE solver returns a non-success retcode. Pass unsafe=true to suppress this check:

xf_unsafe = flow(0.0, [1.0, 0.0], 1.0; unsafe=true)

2-element Vector{Float64}:
 0.36787944127643135
 0.0

Use unsafe=true inside shooting methods or optimisation loops where you want to handle failures gracefully instead of relying on exceptions.

---

SciML integrator internals

The SciML integrator strategy itself — its options, construction, and the CommonSolve.solve method that wraps SciML's solve — is provided by CTSolvers.Integrators. CTFlows contributes only the glue: Integrators.build_problem turns a system and a configuration into a SciML ODEProblem, and Integrators.build_options selects the integrator's cached option bundle for the configuration. Integration is then CommonSolve.solve(prob, integ). Keeping build_problem separate from the solve step lets the same problem definition be re-solved with different parameters efficiently.

integ = Integrators.SciML(; reltol=1e-8)

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.

---

See also

• CTFlows.Configs.StateEndPointConfig, CTFlows.Configs.StateTrajectoryConfig — state configuration objects.

• CTFlows.Configs.HamiltonianEndPointConfig, CTFlows.Configs.HamiltonianTrajectoryConfig — Hamiltonian configuration objects.

• CTFlows.Configs.tspan, CTFlows.Configs.initial_state, CTFlows.Configs.initial_costate — configuration accessors.

• CTFlows.Integrators.build_problem, CTFlows.Integrators.build_options — the CTFlows-side integrator glue.

• CTSolvers.Integrators.SciML, CTSolvers.Integrators.AbstractIntegrator — the integrator strategy (provided by CTSolvers).