Discretisation schemes

continuo solves the perfect-foresight transition as a two-point boundary value problem by global collocation: the trajectory is discretised on a grid of N intervals, and every node value is stacked into one large nonlinear system solved by Newton (see Linear solvers for the linear core). The scheme is how each interval is discretised — it sets the accuracy and stability of the transition for a given grid.

The default is Crank–Nicolson (the implicit midpoint rule): second-order and A-stable, a robust general-purpose choice. Three higher-order collocation families are also available, selectable with their order; on a smooth problem they reach a target accuracy on a far coarser grid.

This is independent of the linear backend (Linear solvers) and the nonlinear steady-state algorithm (Steady-state solvers): a run mixes them freely. Every scheme here produces the same one-step (block-triangular) coupling, so the linear backends and their cross-segment warm-start are unaffected by the choice.

Choosing a scheme

Three entry points, later overriding earlier (CLI > directive > default):

The simulate directive

Pin the scheme (and order) with the model (see Commands):

simulate(T = 120, N = 300, scheme = radau, order = 5);

Python API

Pass scheme= / order= to continuo.Model.simul():

model.simul(scheme="gauss", order=4)

Command line

Override at the call site:

$ continuo model.mod --scheme radau --order 5

order selects the collocation order within a family; omit it for the family default. crank_nicolson is fixed second-order and takes no order.

The families

All four are collocation methods — they differ only in where the s stage nodes sit inside each interval; the order then follows from the node count. The bound on accuracy is the global order p: halving the step cuts the error by about 2^p.

scheme	order	order ∈	stability	use it for
crank_nicolson	2	— (fixed)	A-stable, symmetric	the robust default
gauss	2s	2, 4, 6	A-stable, symmetric	smooth / conservative (orbits)
radau	2s−1	1, 3, 5	L-stable, stiffly accurate	stiff / sharp transitions
lobatto_iiia	2s−2	2, 4, 6	A-stable, endpoint nodes	BVP-native collocation

gauss at order 2 is exactly Crank–Nicolson (the one-stage Gauss method is the implicit midpoint), so the families extend the default rather than replace it. radau is L-stable and stiffly accurate — the right choice when the dynamics are stiff or the path has a sharp transition that A-stable schemes resolve only on a fine grid.

How it works

Each interval carries s internal stage unknowns (the stage derivatives, plus the algebraic stage values for index-1 variables) appended to the stacked vector after the node block. The model residual is composed with the stage relations and differentiated by CasADi automatic differentiation, so the exact Jacobian is preserved and the per-interval coupling stays one-step. Stage unknowns are seeded from the node guess and dropped from the returned path. Crank–Nicolson keeps its dedicated one-stage form.

Accuracy

On a smooth problem the error of a scheme of order p falls like h^p. The Goodwin example makes this concrete on a closed orbit: against a fine reference, the max error over the trajectory is

N	crank_nicolson (2)	gauss (4)	radau (5)
150	2.7e-02	2.3e-05	5.3e-07
300	6.9e-03	1.4e-06	1.6e-08
600	1.7e-03	9.0e-08	5.2e-10
1200	4.3e-04	5.6e-09	1.6e-11

— the higher-order schemes reach at N = 150 an accuracy Crank–Nicolson does not match at N = 1200. Higher order pays off when the solution is smooth; near a kink (an occasionally-binding constraint) or a discontinuous shock the gain is limited by the non-smoothness, and refining N (or Crank–Nicolson) can be the steadier choice.

Note

The schemes above run on a uniform grid by default, but the grid need not be uniform: non-uniform and adaptively refined meshes — placing nodes where the solution varies fastest — build on this same collocation machinery. See Time grids and adaptive refinement.