Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Numerical Integration of Ordinary Differential Equations and Collocation

Numerical optimal-control methods depend fundamentally on numerical methods for ordinary differential equations. Before introducing direct transcription and collocation, it is essential to understand how differential equations are approximated, how time-marching schemes work, and why discretization changes the mathematical problem.

This section develops:

Initial-Value Problems

Consider the ordinary differential equation

x˙(t)=f(t,x(t)),\dot{\boldsymbol{x}}(t) = \boldsymbol{f}(t,\boldsymbol{x}(t)),
(1)

with initial condition

x(t0)=x0.\boldsymbol{x}(t_0) = \boldsymbol{x}_0.
(2)

The exact solution is a continuous function

x:[t0,tf]→Rn.\boldsymbol{x}:[t_0,t_f]\rightarrow\mathbb{R}^n.
(3)

A numerical method replaces this continuous solution by approximations at a finite set of points.

Time Discretization

Introduce a sequence of grid points:

t0<t1<t2<⋯<tN=tf.t_0<t_1<t_2<\cdots<t_N=t_f.
(4)

Define the step size

hk=tk+1−tk.h_k = t_{k+1}-t_k.
(5)

The grid need not be uniform. If all steps are equal,

hk=h,h_k=h,
(6)

the mesh is uniform.

Let

xk≈x(tk).\boldsymbol{x}_k \approx \boldsymbol{x}(t_k).
(7)

A time-marching method computes

x0→x1→x2→⋯→xN.\boldsymbol{x}_0 \rightarrow \boldsymbol{x}_1 \rightarrow \boldsymbol{x}_2 \rightarrow \cdots \rightarrow \boldsymbol{x}_N.
(8)

Finite-Difference Approximation of the Derivative

A forward finite-difference approximation is

x˙(tk)≈xk+1−xkhk.\dot{\boldsymbol{x}}(t_k) \approx \frac{\boldsymbol{x}_{k+1}-\boldsymbol{x}_k}{h_k}.
(9)

Taylor expansion gives

x(tk+1)=x(tk)+hkx˙(tk)+hk22x¨(ξk),\boldsymbol{x}(t_{k+1}) = \boldsymbol{x}(t_k) + h_k\dot{\boldsymbol{x}}(t_k) + \frac{h_k^2}{2}\ddot{\boldsymbol{x}}(\xi_k),
(10)

for some

ξk∈(tk,tk+1).\xi_k\in(t_k,t_{k+1}).
(11)

Therefore,

x(tk+1)−x(tk)hk=x˙(tk)+O(hk).\frac{\boldsymbol{x}(t_{k+1})-\boldsymbol{x}(t_k)}{h_k} = \dot{\boldsymbol{x}}(t_k) + \mathcal{O}(h_k).
(12)

The derivative approximation is first-order accurate.

Forward Euler Method

Substitute the forward difference into the differential equation and evaluate the vector field at the beginning of the interval:

xk+1−xkhk=f(tk,xk).\frac{\boldsymbol{x}_{k+1}-\boldsymbol{x}_k}{h_k} = \boldsymbol{f}(t_k,\boldsymbol{x}_k).
(13)

Rearranging gives

xk+1=xk+hkf(tk,xk).\boxed{ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k + h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k). }
(14)

This is the forward Euler method.

Geometric interpretation

The method uses the tangent direction at the beginning of the interval and extrapolates along a straight line.

Forward Euler extrapolates from the current state along the initial tangent, producing a local discretization error relative to the exact trajectory.

Figure 1:Forward Euler extrapolates from the current state along the initial tangent, producing a local discretization error relative to the exact trajectory.

Advantages

Disadvantages

Local and Global Error

The local truncation error is the error made in one step assuming the starting value is exact.

For forward Euler,

τk+1=O(hk2).\tau_{k+1} = \mathcal{O}(h_k^2).
(15)

The global error after many steps is

ek=x(tk)−xk.\boldsymbol{e}_k = \boldsymbol{x}(t_k)-\boldsymbol{x}_k.
(16)

For a sufficiently smooth and stable problem,

ek=O(h).\boldsymbol{e}_k = \mathcal{O}(h).
(17)

Thus, forward Euler is globally first-order accurate.

Truncation Error and Roundoff Error

Reducing the step size decreases truncation error, but it increases the number of arithmetic operations.

Two competing effects occur:

  1. Large step size: truncation error dominates.

  2. Extremely small step size: roundoff and cancellation errors may dominate.

If

hkf(tk,xk)h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k)
(18)

is smaller than machine resolution relative to xk\boldsymbol{x}_k, the update may numerically satisfy

xk+1=xk,\boldsymbol{x}_{k+1} = \boldsymbol{x}_k,
(19)

even when the true solution is changing.

Backward Euler Method

Evaluate the vector field at the end of the interval:

xk+1−xkhk=f(tk+1,xk+1).\frac{\boldsymbol{x}_{k+1}-\boldsymbol{x}_k}{h_k} = \boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}).
(20)

Therefore,

xk+1=xk+hkf(tk+1,xk+1).\boxed{ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k + h_k\boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}). }
(21)

This is the backward Euler method.

The method is implicit because xk+1\boldsymbol{x}_{k+1} appears on both sides.

Solving the Implicit Equation

Define the residual

R(xk+1)=xk+1−xk−hkf(tk+1,xk+1).\boldsymbol{R}(\boldsymbol{x}_{k+1}) = \boldsymbol{x}_{k+1} - \boldsymbol{x}_k - h_k\boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}).
(22)

At every step, solve

R(xk+1)=0.\boldsymbol{R}(\boldsymbol{x}_{k+1})=\boldsymbol{0}.
(23)

A fixed-point iteration is

xk+1(j)=xk+hkf(tk+1,xk+1(j−1)).\boldsymbol{x}_{k+1}^{(j)} = \boldsymbol{x}_k + h_k\boldsymbol{f} \left( t_{k+1}, \boldsymbol{x}_{k+1}^{(j-1)} \right).
(24)

A simple initial guess is

xk+1(0)=xk,\boldsymbol{x}_{k+1}^{(0)} = \boldsymbol{x}_k,
(25)

or a forward Euler prediction.

Newton’s method may instead be used:

xk+1(j+1)=xk+1(j)−[∂R∂xk+1]−1R(xk+1(j)).\boldsymbol{x}_{k+1}^{(j+1)} = \boldsymbol{x}_{k+1}^{(j)} - \left[ \frac{\partial\boldsymbol{R}}{\partial\boldsymbol{x}_{k+1}} \right]^{-1} \boldsymbol{R} \left( \boldsymbol{x}_{k+1}^{(j)} \right).
(26)

The Jacobian is

∂R∂xk+1=I−hk∂f∂x(tk+1,xk+1).\frac{\partial\boldsymbol{R}}{\partial\boldsymbol{x}_{k+1}} = \boldsymbol{I} - h_k \frac{\partial\boldsymbol{f}}{\partial\boldsymbol{x}} \left( t_{k+1},\boldsymbol{x}_{k+1} \right).
(27)

Explicit and Implicit Schemes

A method is explicit if

xk+1\boldsymbol{x}_{k+1}
(28)

can be computed directly from known quantities.

A method is implicit if the new state satisfies a nonlinear algebraic equation.

Explicit methods are typically cheaper per step.

Implicit methods are typically more stable and are often preferred for stiff problems.

The Generalized θ\theta-Method

A unified one-step family is

xk+1−xkhk=θf(tk,xk)+(1−θ)f(tk+1,xk+1).\frac{\boldsymbol{x}_{k+1}-\boldsymbol{x}_k}{h_k} = \theta\boldsymbol{f}(t_k,\boldsymbol{x}_k) + (1-\theta) \boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}).
(29)

Equivalently,

xk+1=xk+hk[θfk+(1−θ)fk+1].\boxed{ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k + h_k \left[ \theta\boldsymbol{f}_k + (1-\theta)\boldsymbol{f}_{k+1} \right]. }
(30)

Special cases are:

θ=1⇒forward Euler,θ=0⇒backward Euler,θ=12⇒trapezoidal or Crank–Nicolson method.\begin{aligned} \theta=1 &\quad\Rightarrow\quad \text{forward Euler}, \\ \theta=0 &\quad\Rightarrow\quad \text{backward Euler}, \\ \theta=\frac12 &\quad\Rightarrow\quad \text{trapezoidal or Crank–Nicolson method}. \end{aligned}
(31)

The trapezoidal scheme is

xk+1=xk+hk2[f(tk,xk)+f(tk+1,xk+1)].\boxed{ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \frac{h_k}{2} \left[ \boldsymbol{f}(t_k,\boldsymbol{x}_k) + \boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}) \right]. }
(32)

Unlike Euler methods, it is second-order accurate.

Stability Analysis Using the Test Equation

Consider

x˙=λx,\dot{x} = \lambda x,
(33)

with

Re⁡(λ)<0.\operatorname{Re}(\lambda)<0.
(34)

The exact solution decays:

x(t)=x0eλt.x(t) = x_0e^{\lambda t}.
(35)

A numerical method produces

xk+1=R(z)xk,z=hλ,x_{k+1} = R(z)x_k, \qquad z=h\lambda,
(36)

where R(z)R(z) is the stability function.

The method is stable for a given zz if

∣R(z)∣≤1.|R(z)|\leq1.
(37)

Forward Euler

RFE(z)=1+z.R_{\mathrm{FE}}(z) = 1+z.
(38)

Stability requires

∣1+z∣≤1.|1+z|\leq1.
(39)

For real negative λ\lambda,

−2≤hλ≤0.-2\leq h\lambda\leq0.
(40)

Therefore, the step size is restricted.

Backward Euler

RBE(z)=11−z.R_{\mathrm{BE}}(z) = \frac{1}{1-z}.
(41)

Backward Euler is stable for the entire left half-plane:

Re⁡(z)≤0.\operatorname{Re}(z)\leq0.
(42)

It is therefore A-stable.

Trapezoidal rule

RTR(z)=1+z/21−z/2.R_{\mathrm{TR}}(z) = \frac{1+z/2}{1-z/2}.
(43)

It is also A-stable, but unlike backward Euler, it does not strongly damp very fast modes.

Accuracy Versus Stability

Accuracy and stability are distinct concepts.

A method may be:

Backward Euler and forward Euler are both first-order accurate, but backward Euler has much stronger stability properties.

Multiple-Stage Methods

Instead of estimating the slope at only one point, multiple-stage methods estimate slopes at several locations within a step.

Runge–Kutta methods are one-step, multiple-stage methods.

They do not require values from earlier time steps beyond xk\boldsymbol{x}_k, but they use several vector-field evaluations during the interval

..
(44)

Classical Fourth-Order Runge–Kutta Method

Let

h=hk.h=h_k.
(45)

Define:

k1=hf(tk,xk),k2=hf(tk+h2,xk+k12),k3=hf(tk+h2,xk+k22),k4=hf(tk+h,xk+k3).\begin{aligned} \boldsymbol{k}_1 &= h\boldsymbol{f}(t_k,\boldsymbol{x}_k), \\ \boldsymbol{k}_2 &= h\boldsymbol{f} \left( t_k+\frac{h}{2}, \boldsymbol{x}_k+\frac{\boldsymbol{k}_1}{2} \right), \\ \boldsymbol{k}_3 &= h\boldsymbol{f} \left( t_k+\frac{h}{2}, \boldsymbol{x}_k+\frac{\boldsymbol{k}_2}{2} \right), \\ \boldsymbol{k}_4 &= h\boldsymbol{f} \left( t_k+h, \boldsymbol{x}_k+\boldsymbol{k}_3 \right). \end{aligned}
(46)

The update is

xk+1=xk+16(k1+2k2+2k3+k4).\boxed{ \boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \frac{1}{6} \left( \boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4 \right). }
(47)

The method requires four vector-field evaluations per step.

Its global error is

O(h4).\mathcal{O}(h^4).
(48)

Interpretation of RK4

RK4 computes:

The weighted average

16(k1+2k2+2k3+k4)\frac{1}{6} \left( \boldsymbol{k}_1+2\boldsymbol{k}_2+2\boldsymbol{k}_3+\boldsymbol{k}_4 \right)
(49)

approximates the integrated effect of the vector field across the interval.

Computational Cost

A higher-order method typically reduces error for a given step size but requires more work per step.

The relevant cost is not simply the number of stages.

One should compare:

total computational cost=cost per step×number of steps required.\text{total computational cost} = \text{cost per step} \times \text{number of steps required}.
(50)

A fourth-order method may be cheaper overall because it permits much larger steps for the same accuracy.

The Discrete Problem Is Not the Continuous Problem

The continuous problem is

x˙=f(t,x).\dot{\boldsymbol{x}} = \boldsymbol{f}(t,\boldsymbol{x}).
(51)

After discretization, the numerical method defines a different problem:

xk+1=Ψ(tk,xk,hk).\boldsymbol{x}_{k+1} = \boldsymbol{\Psi} \left( t_k,\boldsymbol{x}_k,h_k \right).
(52)

The step size hkh_k is a parameter that does not appear in the original differential equation.

For each selected mesh, the discrete problem generally has a different solution.

Therefore, one should say:

The computation provides a numerical solution of a discrete approximation to the differential equation.

It does not provide the exact continuous solution unless a special exactness property happens to hold.

Consistency, Stability, and Convergence

Three concepts are fundamental.

Consistency

A method is consistent if the local truncation error vanishes as

h→0.h\rightarrow0.
(53)

Stability

A method is stable if perturbations do not grow uncontrollably under the discrete propagation.

Convergence

A method is convergent if

xk→x(tk)\boldsymbol{x}_k \rightarrow \boldsymbol{x}(t_k)
(54)

as the mesh is refined.

For standard linear multistep settings, consistency plus zero-stability implies convergence. For one-step methods, analogous stability and consistency results apply.

Sequential Form of Forward Euler

The forward Euler equations are:

x1=x0+h0f0,x2=x1+h1f1,⋮xk+1=xk+hkfk.\begin{aligned} \boldsymbol{x}_1 &= \boldsymbol{x}_0+h_0\boldsymbol{f}_0, \\ \boldsymbol{x}_2 &= \boldsymbol{x}_1+h_1\boldsymbol{f}_1, \\ &\vdots \\ \boldsymbol{x}_{k+1} &= \boldsymbol{x}_k+h_k\boldsymbol{f}_k. \end{aligned}
(55)

In time marching, each equation is used after the preceding state has been computed.

Residual Form

Each equation can instead be written as a residual:

rk=xk+1−xk−hkf(tk,xk)=0.\boldsymbol{r}_k = \boldsymbol{x}_{k+1} - \boldsymbol{x}_k - h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k) = \boldsymbol{0}.
(56)

Stack the unknown states:

y=[x0x1⋮xN].\boldsymbol{y} = \begin{bmatrix} \boldsymbol{x}_0\\ \boldsymbol{x}_1\\ \vdots\\ \boldsymbol{x}_N \end{bmatrix}.
(57)

Stack all residuals:

G(y,h)=[r0r1⋮rN−1]=0.\boxed{ \boldsymbol{G}(\boldsymbol{y},\boldsymbol{h}) = \begin{bmatrix} \boldsymbol{r}_0\\ \boldsymbol{r}_1\\ \vdots\\ \boldsymbol{r}_{N-1} \end{bmatrix} = \boldsymbol{0}. }
(58)

This is a system of nonlinear algebraic equations.

Simultaneous Solution

Instead of computing the states sequentially, solve all state values simultaneously:

G(y,h)=0.\boldsymbol{G}(\boldsymbol{y},\boldsymbol{h})=\boldsymbol{0}.
(59)

Possible root-finding methods include:

For forward Euler, sequential and simultaneous solutions are mathematically equivalent when solved exactly.

However, the simultaneous viewpoint is the conceptual bridge to direct transcription.

Backward Euler in Residual Form

Backward Euler gives

rk=xk+1−xk−hkf(tk+1,xk+1)=0.\boldsymbol{r}_k = \boldsymbol{x}_{k+1} - \boldsymbol{x}_k - h_k\boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}) = \boldsymbol{0}.
(60)

Stacking all intervals yields

GBE(y,h)=0.\boldsymbol{G}_{\mathrm{BE}}(\boldsymbol{y},\boldsymbol{h}) = \boldsymbol{0}.
(61)

The entire trajectory can be solved as one algebraic system.

This avoids performing a separate nonlinear iteration at each time step, although the global nonlinear solver still performs iterations.

Introduction to Collocation

Collocation enforces equality between:

at selected points.

For forward Euler:

xk+1−xkhk−f(tk,xk)=0.\boxed{ \frac{\boldsymbol{x}_{k+1}-\boldsymbol{x}_k}{h_k} - \boldsymbol{f}(t_k,\boldsymbol{x}_k) = \boldsymbol{0}. }
(62)

For backward Euler:

xk+1−xkhk−f(tk+1,xk+1)=0.\boxed{ \frac{\boldsymbol{x}_{k+1}-\boldsymbol{x}_k}{h_k} - \boldsymbol{f}(t_{k+1},\boldsymbol{x}_{k+1}) = \boldsymbol{0}. }
(63)

These equations force the derivative approximation and vector field to coincide at the selected collocation point.

Defect Constraints

The residual

Δk=xk+1−xk−hkf(tk,xk)\boldsymbol{\Delta}_k = \boldsymbol{x}_{k+1} - \boldsymbol{x}_k - h_k\boldsymbol{f}(t_k,\boldsymbol{x}_k)
(64)

is often called a defect.

The discrete dynamics are enforced by

Δk=0.\boxed{ \boldsymbol{\Delta}_k=\boldsymbol{0}. }
(65)

In direct optimal control, state values are optimization variables and defects become equality constraints in a nonlinear program.

Time Marching Versus Collocation

Time marching

Collocation

Connection to Direct Optimal Control

In a direct transcription method, one may optimize over:

{x0,…,xN,u0,…,uN}.\left\{ \boldsymbol{x}_0,\ldots,\boldsymbol{x}_N, \boldsymbol{u}_0,\ldots,\boldsymbol{u}_N \right\}.
(66)

The differential equations become algebraic defect constraints:

Δk(xk,xk+1,uk,uk+1)=0.\boldsymbol{\Delta}_k(\boldsymbol{x}_k,\boldsymbol{x}_{k+1},\boldsymbol{u}_k,\boldsymbol{u}_{k+1})=\boldsymbol{0}.
(67)

The objective is discretized:

J≈JN.J \approx J_N.
(68)

Path constraints are enforced at nodes or collocation points:

c(xk,uk,tk)≤0.\boldsymbol{c}(\boldsymbol{x}_k,\boldsymbol{u}_k,t_k)\leq\boldsymbol{0}.
(69)

Thus, the continuous optimal-control problem becomes a sparse nonlinear programming problem.

Common Errors

  1. Calling forward Euler accurate merely because the step is small.

  2. Ignoring roundoff when reducing the step size.

  3. Confusing local and global truncation error.

  4. Assuming implicit methods are automatically higher order.

  5. Assuming stability and accuracy are the same.

  6. Comparing methods only by the number of function evaluations.

  7. Saying a numerical result is the exact solution of the differential equation.

  8. Forgetting that the step size changes the discrete problem.

  9. Treating collocation as unrelated to finite differences.

  10. Ignoring defect convergence under mesh refinement.

Verification Guidelines

A numerical ODE solution should be checked by:

  1. reducing the step size;

  2. comparing multiple integration methods;

  3. monitoring local error estimates;

  4. checking invariants or conserved quantities;

  5. comparing with analytic solutions when available;

  6. evaluating residuals;

  7. inspecting stability restrictions;

  8. testing variable-step methods;

  9. checking sensitivity to tolerances; and

  10. documenting convergence rates.

Summary

The main conclusions are:

  1. Numerical integration approximates a continuous ODE on a discrete mesh.

  2. Forward Euler is explicit, simple, and first-order accurate.

  3. Backward Euler is implicit, first-order accurate, and much more stable.

  4. The θ\theta-method unifies forward Euler, backward Euler, and trapezoidal integration.

  5. RK4 uses four stages and provides fourth-order global accuracy.

  6. Truncation, roundoff, stability, and cost must all be balanced.

  7. A discretized equation is a different mathematical problem from the original ODE.

  8. Time-marching equations can be rewritten as simultaneous algebraic residuals.

  9. Enforcing derivative approximations and vector fields at common points leads to collocation.

  10. Collocation is a central foundation of direct numerical optimal-control methods.

Connection. A time-marching formula can be imposed simultaneously rather than sequentially; Euler collocation makes this algebraic viewpoint explicit.

Chapter 4: Direct Numerical Methods and Collocation
Direct Numerical Methods and Direct Shooting
Chapter 4: Direct Numerical Methods and Collocation
Direct Collocation with Euler Discretization