In system dynamics and control engineering, most mathematical models are expressed as initial-value problems (IVPs) for ordinary differential equations (ODEs) of the form

\[ y'(t) = f(t, y(t)), \quad t \in [t_0, T], \quad y(t_0) = y_0, \]

where \( f : \mathbb{R} \times \mathbb{R}^m \to \mathbb{R}^m \) is sufficiently smooth. Analytic solutions are rarely available, so we approximate the continuous-time trajectory \( y(t) \) on a discrete grid

\[ t_n = t_0 + n h, \quad n = 0,1,\dots,N, \quad h = \frac{T - t_0}{N}, \]

by a discrete sequence \( \{y_n\}_{n=0}^N \). A numerical scheme (or time-stepping method) is a recursion that maps \( y_n \) to \( y_{n+1} \). Understanding how errors arise and propagate in this recursion is essential for:

• Choosing appropriate step sizes \( h \) for simulation accuracy.
• Ensuring that numerical trajectories respect the stability of the underlying continuous system.
• Comparing integration algorithms (e.g., explicit vs. implicit) in control-oriented simulations.

flowchart TD
  A["Continuous-time ODE model y' = f(t,y)"] --> B["Choose grid t_n = t0 + n*h"]
  B --> C["Define discrete scheme y_{n+1} = Phi_h(t_n, y_n)"]
  C --> D["One-step defect (local truncation error)"]
  D --> E["Accumulated discrepancy (global error)"]
  E --> F["Analyze consistency, stability, convergence"]
        

In later chapters, these concepts will underpin numerical simulation of mechanical, electrical, and multi-domain systems in tools such as Python/SciPy, MATLAB/Simulink, and control libraries.

A broad and important class of numerical integrators for IVPs are one-step methods. They have the general form

\[ y_{n+1} = \Phi_h(t_n, y_n), \quad n = 0,1,\dots,N-1, \]

where \( \Phi_h \) is called the update or step operator. It may depend on \( f \), on \( h \), and possibly on evaluations of \( f \) at intermediate points between \( t_n \) and \( t_{n+1} \).

Important examples that we will use as running test cases:

• Explicit (forward) Euler:

  \[ y_{n+1} = y_n + h f(t_n, y_n). \]

• Implicit (backward) Euler:

  \[ y_{n+1} = y_n + h f(t_{n+1}, y_{n+1}), \]

  which requires solving an algebraic equation for \( y_{n+1} \) at each step.
• Trapezoidal rule (Crank–Nicolson):

  \[ y_{n+1} = y_n + \frac{h}{2}\bigl(f(t_n, y_n) + f(t_{n+1}, y_{n+1})\bigr). \]

These schemes approximate the action of the exact flow map \( \varphi_h \) defined by

\[ \varphi_h(t_n, y(t_n)) = y(t_n + h), \]

and we will measure their quality and robustness by comparing \( \Phi_h \) to \( \varphi_h \) in various ways.

The local truncation error (LTE) measures the defect of a numerical scheme in a single step, assuming that the starting value of the step is exact. Let \( y(t) \) be the exact solution of the IVP. For a one-step method \( y_{n+1} = \Phi_h(t_n,y_n) \) we define

\[ \delta_{n+1}(h) = y(t_{n+1}) - \Phi_h\bigl(t_n, y(t_n)\bigr). \]

This is the per-step defect introduced by replacing the exact flow \( \varphi_h \) with the numerical step \( \Phi_h \), when starting on the exact trajectory.

A one-step method is said to be consistent of order \( p \ge 1 \) if there exists a constant \( C > 0 \) independent of \( h \) such that

\[ \max_{0 \le n \le N-1} |\delta_{n+1}(h)| \le C h^{p+1} \quad \text{as } h \to 0. \]

The exponent \( p+1 \) characterizes the rate at which the one-step defect vanishes as the step size is refined.

3.1 Example: LTE of Forward Euler

For forward Euler

\[ y_{n+1} = y_n + h f(t_n,y_n), \]

and the exact solution, Taylor expansion around \( t_n \) gives

\[ y(t_{n+1}) = y(t_n) + h y'(t_n) + \frac{h^2}{2} y''(\xi_n), \quad \xi_n \in (t_n, t_{n+1}). \]

Since \( y'(t_n) = f(t_n, y(t_n)) \), we can rewrite this as

\[ y(t_{n+1}) = y(t_n) + h f(t_n, y(t_n)) + \frac{h^2}{2} y''(\xi_n). \]

The local truncation error of Euler is therefore

\[ \delta_{n+1}(h) = y(t_{n+1}) - \bigl(y(t_n) + h f(t_n,y(t_n))\bigr) = \frac{h^2}{2} y''(\xi_n) = \mathcal{O}(h^2), \]

provided \( y''(t) \) remains bounded. Hence forward Euler is consistent of order \( p = 1 \).

While LTE measures a defect per step, the global error at grid point \( t_n \) is

\[ e_n = y(t_n) - y_n. \]

A numerical method is called convergent of order \( p \) if there exists \( C > 0 \) independent of \( h \) such that

\[ \max_{0 \le n \le N} |e_n| \le C h^p \quad \text{as } h \to 0. \]

Thus, the global approximation error decays like \( h^p \) when the grid is refined.

4.1 Error Recursion

For a general one-step method with LTE \( \delta_{n+1}(h) \), we can write

\[ \begin{aligned} e_{n+1} &= y(t_{n+1}) - y_{n+1} \\ &= y(t_{n+1}) - \Phi_h(t_n, y_n) \\ &= \bigl(y(t_{n+1}) - \Phi_h(t_n, y(t_n))\bigr) + \bigl(\Phi_h(t_n, y(t_n)) - \Phi_h(t_n, y_n)\bigr) \\ &= \delta_{n+1}(h) + \bigl[\Phi_h(t_n, y(t_n)) - \Phi_h(t_n, y_n)\bigr]. \end{aligned} \]

If \( \Phi_h \) is Lipschitz in its second argument with constant \( L_h \), i.e.,

\[ |\Phi_h(t, v) - \Phi_h(t, w)| \le L_h |v - w|, \]

then

\[ |e_{n+1}| \le |\delta_{n+1}(h)| + L_h |e_n|. \]

The factors \( L_h \) describe how perturbations are amplified step by step and are thus directly related to stability. Roughly:

• Consistency controls the size of the local errors \( \delta_{n+1} \).
• Stability controls how these errors are amplified or damped as they propagate.

4.2 Consistency + Stability ⇒ Convergence (Idea)

A fundamental result of numerical analysis states that, under mild conditions, for one-step (and more generally linear multistep) methods:

Theorem (informal). If a method is consistent of order \( p \) and stable (in a precise sense described below), then the global error satisfies

\[ \max_{0 \le n \le N} |e_n| \le C h^p \]

for sufficiently small \( h \), i.e., the method is convergent of order \( p \).

The proof relies on inequalities of Grönwall type applied to the error recursion and will be revisited (in a simplified form) when we discuss specific algorithms in later numerical chapters.

In control engineering, continuous-time stability of a linear system \( x'(t) = A x(t) \) is characterized by the locations of the eigenvalues of \( A \) in the complex plane (e.g., all having negative real parts). For numerical schemes, we must additionally ensure that the discretized recursion does not introduce artificial instability.

5.1 Test Equation and Stability Function

Stability analysis is typically carried out on the scalar test equation

\[ y'(t) = \lambda y(t), \quad \lambda \in \mathbb{C}, \]

whose exact solution is \( y(t) = e^{\lambda (t - t_0)} y_0 \). If \( \operatorname{Re}(\lambda) < 0 \), then \( |y(t)| \) decays to zero as \( t \to \infty \).

Applying a one-step method to this test equation yields a recursion of the form

\[ y_{n+1} = R(z) \, y_n, \quad z = h \lambda, \]

where \( R(z) \) is the stability function. The set

\[ \mathcal{S} = \{ z \in \mathbb{C} : |R(z)| \le 1 \} \]

is called the region of absolute stability of the method. A discretization is said to be stable for a given \( \lambda \) and step \( h \) if \( z = h\lambda \in \mathcal{S} \), ensuring that \( |y_n| \) does not grow without bound as \( n \to \infty \).

5.2 Examples: Euler and Trapezoidal Schemes

• Forward Euler. Applied to \( y' = \lambda y \):

  \[ y_{n+1} = y_n + h \lambda y_n = (1 + z) y_n, \quad z = h\lambda, \]

  so

  \[ R_{\text{FE}}(z) = 1 + z, \quad \mathcal{S}_{\text{FE}} = \{ z \in \mathbb{C} : |1 + z| \le 1 \}, \]

  the closed disk of radius \( 1 \) centered at \( -1 \).
• Backward Euler. Applied to \( y' = \lambda y \):

  \[ y_{n+1} = y_n + h \lambda y_{n+1} \quad\Rightarrow\quad (1 - z) y_{n+1} = y_n, \]

  so

  \[ R_{\text{BE}}(z) = \frac{1}{1 - z}. \]

• Trapezoidal rule. Applying \( y_{n+1} = y_n + \frac{h}{2}(\lambda y_n + \lambda y_{n+1}) \) gives

  \[ \left(1 - \frac{z}{2}\right) y_{n+1} = \left(1 + \frac{z}{2}\right) y_n, \quad R_{\text{TR}}(z) = \frac{1 + z/2}{1 - z/2}. \]

5.3 A-Stability

A method is called A-stable if its stability region contains the entire left half-plane:

\[ \{ z \in \mathbb{C} : \operatorname{Re}(z) \le 0 \} \subseteq \mathcal{S}. \]

For backward Euler,

\[ |R_{\text{BE}}(z)| = \frac{1}{|1 - z|}. \]

If \( \operatorname{Re}(z) \le 0 \), then

\[ |1 - z|^2 = (1 - \operatorname{Re}(z))^2 + (\operatorname{Im}(z))^2 \ge 1, \]

so \( |1 - z| \ge 1 \) and therefore \( |R_{\text{BE}}(z)| \le 1 \). Hence backward Euler is A-stable. In contrast, forward Euler is stable only in the disk \( \mathcal{S}_{\text{FE}} \), so large negative \( \lambda \) require very small \( h \).

flowchart TD
  T["Test equation: y' = lambda y"] --> M["Apply scheme: get y_{n+1} = R(z) y_n"]
  M --> Z["Define z = lambda*h"]
  Z --> C["Solve condition |R(z)| <= 1"]
  C --> H["Choose h so all eigenvalues lambda*h lie in stability region"]
        

Consider a linear time-invariant (LTI) system

\[ x'(t) = A x(t), \quad A \in \mathbb{R}^{m \times m}. \]

For such systems, continuous-time stability is governed by eigenvalues \( \lambda_j \) of \( A \): if \( \operatorname{Re}(\lambda_j) < 0 \) for all \( j \), the origin is asymptotically stable.

Applying a one-step method with stability function \( R(z) \) to the linear system yields a discrete-time recursion of the form

\[ x_{n+1} = \Psi_h x_n, \]

where \( \Psi_h \) is a matrix that depends on \( A \) and \( h \). For forward Euler,

\[ x_{n+1} = \bigl(I + h A\bigr) x_n. \]

If \( A \) is diagonalizable, each eigenmode evolves as

\[ x_{n+1}^{(j)} = R(h\lambda_j)\,x_n^{(j)}. \]

Discrete-time stability requires \( |R(h\lambda_j)| \le 1 \) for all eigenvalues \( \lambda_j \). Thus, in control-oriented simulation, step-size selection is constrained by the requirement that \( h \lambda_j \in \mathcal{S} \) for all modes of interest.

We implement the forward Euler method for the scalar test equation \( y'(t) = \lambda y(t) \) with \( \lambda = -2 \), \( y(0) = 1 \), on \( [0,1] \). The exact solution is \( y(t) = e^{-2t} \). We compute both:

• the global error \( |y(t_n) - y_n| \),
• the local defect (LTE) by applying the scheme to the exact values.

7.1 Python (NumPy, SciPy Ecosystem)

In Python, scientific computing for system dynamics typically uses numpy for arrays and scipy.integrate (e.g., solve_ivp) for higher-order ODE solvers. Here we implement Euler explicitly and use the analytic solution as reference.

import numpy as np

def forward_euler(f, t0, y0, h, N):
    t = np.linspace(t0, t0 + N * h, N + 1)
    y = np.zeros_like(t)
    y[0] = y0
    for n in range(N):
        y[n + 1] = y[n] + h * f(t[n], y[n])
    return t, y

lam = -2.0

def f(t, y):
    return lam * y

t0 = 0.0
y0 = 1.0
h = 0.1
N = int((1.0 - t0) / h)

t, y_num = forward_euler(f, t0, y0, h, N)
y_exact = np.exp(lam * t)

global_error = np.abs(y_exact - y_num)

# Local truncation error: apply one step of Euler to the exact solution
local_error = np.zeros_like(t)
local_error[1:] = np.abs(
    y_exact[1:] - (y_exact[:-1] + h * f(t[:-1], y_exact[:-1]))
)

print("max global error:", np.max(global_error))
print("max local error:", np.max(local_error))
      

7.2 C++ (from Scratch, with Note on odeint)

In C++, one may either implement integrators manually (as below) or use libraries such as boost::numeric::odeint for higher-order methods and adaptive step sizing.

#include <iostream>
#include <vector>
#include <cmath>

double f(double t, double y) {
    const double lambda = -2.0;
    return lambda * y;
}

int main() {
    double t0 = 0.0;
    double y0 = 1.0;
    double h  = 0.1;
    int N = static_cast<int>((1.0 - t0) / h);

    std::vector<double> t(N + 1), y(N + 1), y_exact(N + 1);
    t[0] = t0;
    y[0] = y0;
    y_exact[0] = std::exp(-2.0 * t0);

    for (int n = 0; n < N; ++n) {
        t[n + 1] = t0 + (n + 1) * h;
        y[n + 1] = y[n] + h * f(t[n], y[n]);
        y_exact[n + 1] = std::exp(-2.0 * t[n + 1]);
    }

    double max_global_error = 0.0;
    for (int n = 0; n <= N; ++n) {
        double e = std::fabs(y_exact[n] - y[n]);
        if (e > max_global_error) {
            max_global_error = e;
        }
    }

    std::cout << "max global error = " << max_global_error << '\n';
    return 0;
}
      

7.3 Java (Manual Integrator, Apache Commons Math)

In Java, the Apache Commons Math library provides ODE solvers (e.g., Runge–Kutta). Here we show a basic Euler implementation.

public class ForwardEulerExample {
    static double lambda = -2.0;

    static double f(double t, double y) {
        return lambda * y;
    }

    public static void main(String[] args) {
        double t0 = 0.0;
        double y0 = 1.0;
        double h  = 0.1;
        int N = (int) ((1.0 - t0) / h);

        double[] t = new double[N + 1];
        double[] y = new double[N + 1];
        double[] yExact = new double[N + 1];

        t[0] = t0;
        y[0] = y0;
        yExact[0] = Math.exp(lambda * t0);

        for (int n = 0; n < N; ++n) {
            t[n + 1] = t0 + (n + 1) * h;
            y[n + 1] = y[n] + h * f(t[n], y[n]);
            yExact[n + 1] = Math.exp(lambda * t[n + 1]);
        }

        double maxGlobalError = 0.0;
        for (int n = 0; n <= N; ++n) {
            double e = Math.abs(yExact[n] - y[n]);
            if (e > maxGlobalError) {
                maxGlobalError = e;
            }
        }

        System.out.println("max global error = " + maxGlobalError);
    }
}
      

7.4 MATLAB / Simulink

MATLAB offers built-in ODE solvers such as ode45 (Runge–Kutta) that automatically control local error. Below is a simple forward Euler implementation.

lambda = -2.0;
f = @(t, y) lambda * y;

t0 = 0.0;
y0 = 1.0;
h  = 0.1;
T  = 1.0;

N = round((T - t0) / h);
t = linspace(t0, T, N + 1);
y = zeros(1, N + 1);
y(1) = y0;

for n = 1:N
    y(n + 1) = y(n) + h * f(t(n), y(n));
end

y_exact = exp(lambda * t);
global_error = abs(y_exact - y);

fprintf('max global error = %g\n', max(global_error));

% Reference integration with an embedded Runge-Kutta solver:
% [t_ref, y_ref] = ode45(@(t, y) lambda * y, [t0 T], y0);
      

In Simulink, forward Euler can be implemented using:

• A Gain block with gain \( h \lambda \),
• A Sum block implementing y + h*lambda*y,
• A Unit Delay or Memory block to store \( y_n \) across steps.

7.5 Wolfram Mathematica

Mathematica offers NDSolve for high-order, adaptive integration. The code below shows a manual Euler implementation for comparison.

lambda = -2.0;
f[t_, y_] := lambda y;

t0 = 0.0;
y0 = 1.0;
h  = 0.1;
tFinal = 1.0;
nSteps = Round[(tFinal - t0)/h];

tList = Table[t0 + n*h, {n, 0, nSteps}];
yList = ConstantArray[0.0, nSteps + 1];
yList[[1]] = y0;

Do[
  yList[[n + 1]] = yList[[n]] + h*f[tList[[n]], yList[[n]]],
  {n, 1, nSteps}
];

yExact = Exp[lambda*tList];
globalError = Abs[yExact - yList];

Max[globalError]
      

For production simulations, NDSolve (Mathematica), ode45 (MATLAB), and solve_ivp (SciPy) use higher-order, adaptive schemes precisely to control local and global error.

Problem 1 (LTE of Forward Euler for \( y' = y \)). Consider the IVP \( y'(t) = y(t), \; y(0) = 1 \). Apply the forward Euler method with step size \( h \). Compute the local truncation error \( \delta_{n+1}(h) \) and show that it is of order \( \mathcal{O}(h^2) \).

Solution. The exact solution is \( y(t) = e^t \). With grid \( t_n = n h \), forward Euler gives

\[ y_{n+1} = y_n + h y_n. \]

The local truncation error is

\[ \begin{aligned} \delta_{n+1}(h) &= y(t_{n+1}) - \bigl(y(t_n) + h y(t_n)\bigr) \\ &= e^{t_{n+1}} - e^{t_n}(1 + h). \end{aligned} \]

Since \( t_{n+1} = t_n + h \), we have

\[ e^{t_{n+1}} = e^{t_n + h} = e^{t_n} e^{h} = e^{t_n}\Bigl(1 + h + \tfrac{h^2}{2} + \mathcal{O}(h^3)\Bigr). \]

Therefore

\[ \begin{aligned} \delta_{n+1}(h) &= e^{t_n}\Bigl(1 + h + \tfrac{h^2}{2} + \mathcal{O}(h^3)\Bigr) - e^{t_n}(1 + h) \\ &= e^{t_n}\Bigl(\tfrac{h^2}{2} + \mathcal{O}(h^3)\Bigr) = \mathcal{O}(h^2). \end{aligned} \]

Thus forward Euler has local truncation error of order \( h^2 \).

Problem 2 (Global Error Order for \( y' = \lambda y \)). Consider \( y'(t) = \lambda y(t), \; y(0) = y_0 \) with constant \( \lambda \in \mathbb{C} \). Show that forward Euler has global error of order \( \mathcal{O}(h) \) at fixed final time \( T \).

Solution. For forward Euler,

\[ y_{n+1} = (1 + h\lambda) y_n, \quad y_0 \text{ given}. \]

Solving the recursion yields

\[ y_n = (1 + h\lambda)^n y_0. \]

With \( t_n = n h \), the exact solution is \( y(t_n) = e^{\lambda t_n} y_0 \). The global error is

\[ e_n = y(t_n) - y_n = \bigl(e^{\lambda t_n} - (1 + h\lambda)^n\bigr) y_0. \]

For small \( h \), using the expansion \( \log(1 + h\lambda) = h\lambda + \mathcal{O}(h^2) \) we have

\[ \begin{aligned} (1 + h\lambda)^n &= \exp\bigl(n \log(1 + h\lambda)\bigr) \\ &= \exp\bigl(n h\lambda + n \mathcal{O}(h^2)\bigr) = \exp\bigl(\lambda t_n + \mathcal{O}(n h^2)\bigr). \end{aligned} \]

For fixed \( T \), we have \( n h = t_n \le T \), so \( n h^2 = t_n h = \mathcal{O}(h) \). Thus,

\[ (1 + h\lambda)^n = e^{\lambda t_n}\bigl(1 + \mathcal{O}(h)\bigr), \]

and therefore

\[ e_n = e^{\lambda t_n} y_0 \cdot \mathcal{O}(h) = \mathcal{O}(h). \]

Hence the global error of forward Euler is of order \( h \).

Problem 3 (Absolute Stability Region of Forward Euler). Derive the region of absolute stability \( \mathcal{S}_{\text{FE}} \) for the forward Euler method.

Solution. For the test equation \( y' = \lambda y \), forward Euler gives

\[ y_{n+1} = (1 + z) y_n, \quad z = h\lambda. \]

Stability requires \( |1 + z| \le 1 \). Hence the region of absolute stability is

\[ \mathcal{S}_{\text{FE}} = \{ z \in \mathbb{C} : |1 + z| \le 1 \}, \]

the closed disk of radius \( 1 \) centered at \( -1 \) in the complex plane. Note that only a bounded portion of the left half-plane belongs to this region.

Problem 4 (Step-Size Restriction from Stability). Suppose we discretize \( y' = \lambda y \) with real \( \lambda < 0 \) using forward Euler. Derive a condition on \( h \) for stability, and determine the maximal stable step size for \( \lambda = -20 \).

Solution. For real \( \lambda < 0 \), we have \( z = h\lambda \in \mathbb{R} \). Stability requires

\[ |1 + z| < 1. \]

For real \( z \), this inequality is equivalent to

\[ -1 < 1 + z < 1. \]

Subtracting \( 1 \) yields

\[ -2 < z < 0. \]

With \( z = h\lambda \) and \( \lambda < 0 \), this becomes

\[ -2 < h\lambda < 0. \]

For \( \lambda = -20 \) we have \( h\lambda = -20h \), so

\[ -2 < -20h < 0. \]

Dividing by \( -20 \) reverses the inequalities:

\[ 0 < h < \frac{1}{10}. \]

Therefore, for \( \lambda = -20 \), forward Euler is stable only for \( h < 0.1 \).

• Time-stepping schemes approximate the exact flow map of an ODE on a discrete grid \( t_n = t_0 + n h \).
• The local truncation error measures the one-step defect assuming an exact starting value and is used to define the order of consistency.
• The global error quantifies the accumulated deviation from the exact trajectory; under stability, its order matches the consistency order.
• Numerical stability is analyzed via the test equation \( y' = \lambda y \) and the stability function \( R(z) \), whose magnitude determines the region of absolute stability \( \mathcal{S} \).
• Forward Euler has first-order convergence and a bounded stability region \( \mathcal{S}_{\text{FE}} \), while backward Euler and the trapezoidal rule are A-stable, making them attractive for stiff or fast-decaying modes.
• Implementations in Python, C++, Java, MATLAB/Simulink, and Mathematica all realize the same underlying recursion; their libraries provide higher-order and adaptive methods that use these concepts to control error and ensure stable simulations of dynamic systems.