We consider the continuous-time LTI state equation (from earlier lessons in this chapter): \( \dot{\mathbf{x}}(t) \): the state derivative, \( \mathbf{x}(t)\in\mathbb{R}^n \): the state, \( \mathbf{u}(t)\in\mathbb{R}^m \): the input.

\[ \dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t), \quad \mathbf{x}(t_0)=\mathbf{x}_0 \]

The exact solution (Lessons 1–2) is \( \mathbf{x}(t) \): the state trajectory.

\[ \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 + \int_{t_0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\mathbf{u}(\tau)\,d\tau. \]

In practice, we often need numerical simulation because (i) \( \mathbf{u}(t) \) may be complex or only available as samples, (ii) we may need fast repeated simulations, and (iii) we may want a discrete-time recursion for implementation and verification.

Write the dynamics as a general ODE: \( \dot{\mathbf{x}}(t)=\mathbf{f}(t,\mathbf{x}(t)) \): where here \( \mathbf{f}(t,\mathbf{x})=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}(t) \). Choose a step size \( h > 0 \) and grid \( t_k=t_0+k h \). A one-step method has the form

\[ \mathbf{x}_{k+1}=\mathbf{x}_k + h\,\boldsymbol{\Phi}(t_k,\mathbf{x}_k;h), \quad \mathbf{x}_0=\mathbf{x}(t_0). \]

2.1 Forward Euler

Forward Euler uses \( \boldsymbol{\Phi}(t_k,\mathbf{x}_k;h)=\mathbf{f}(t_k,\mathbf{x}_k) \), hence

\[ \mathbf{x}_{k+1}=\mathbf{x}_k + h\,\mathbf{f}(t_k,\mathbf{x}_k). \]

2.2 Local truncation error of Euler (rigorous Taylor proof)

Define the (one-step) local truncation error by applying the numerical formula to the true solution: \( \boldsymbol{\tau}_{k+1} \):

\[ \boldsymbol{\tau}_{k+1} = \mathbf{x}(t_{k+1}) - \Big(\mathbf{x}(t_k) + h\,\mathbf{f}(t_k,\mathbf{x}(t_k))\Big). \]

Since \( \mathbf{x}(t) \) is differentiable (guaranteed for LTI with continuous input), Taylor’s theorem with remainder gives, for some \( \xi_k\in(t_k,t_{k+1}) \):

\[ \mathbf{x}(t_{k+1}) = \mathbf{x}(t_k) + h\,\dot{\mathbf{x}}(t_k) + \frac{h^2}{2}\,\ddot{\mathbf{x}}(\xi_k). \]

Using \( \dot{\mathbf{x}}(t_k)=\mathbf{f}(t_k,\mathbf{x}(t_k)) \) yields

\[ \boldsymbol{\tau}_{k+1} = \frac{h^2}{2}\,\ddot{\mathbf{x}}(\xi_k), \quad \Rightarrow \quad \|\boldsymbol{\tau}_{k+1}\| \le \frac{h^2}{2}\max_{\tau\in[t_k,t_{k+1}]}\|\ddot{\mathbf{x}}(\tau)\|. \]

Hence Euler has \( \boldsymbol{\tau}_{k+1}=\mathcal{O}(h^2) \) local truncation error.

2.3 Global error order (Grönwall-style bound)

Let the global error be \( \mathbf{e}_k=\mathbf{x}(t_k)-\mathbf{x}_k \). Assume \( \mathbf{f} \) is Lipschitz in \( \mathbf{x} \) with constant \( L \) on the trajectory region:

\[ \|\mathbf{f}(t,\mathbf{x})-\mathbf{f}(t,\mathbf{y})\| \le L\|\mathbf{x}-\mathbf{y}\|. \]

Standard analysis (Euler stability recursion) gives

\[ \|\mathbf{e}_{k+1}\| \le (1+hL)\|\mathbf{e}_k\| + C h^2, \]

where \( C \) bounds \( \|\ddot{\mathbf{x}}(t)\| \). Iterating the inequality and using \( (1+hL)^k \le e^{L(t_k-t_0)} \) yields

\[ \|\mathbf{e}_k\| \le \frac{C}{L}\left(e^{L(t_k-t_0)}-1\right)h \quad \Rightarrow \quad \mathbf{e}_k=\mathcal{O}(h). \]

Therefore, forward Euler is first-order accurate globally.

2.4 Classical Runge–Kutta 4 (RK4)

RK4 uses four stages:

\[ \begin{aligned} \mathbf{k}_1 &= \mathbf{f}(t_k,\mathbf{x}_k),\\ \mathbf{k}_2 &= \mathbf{f}(t_k+\tfrac{h}{2},\mathbf{x}_k+\tfrac{h}{2}\mathbf{k}_1),\\ \mathbf{k}_3 &= \mathbf{f}(t_k+\tfrac{h}{2},\mathbf{x}_k+\tfrac{h}{2}\mathbf{k}_2),\\ \mathbf{k}_4 &= \mathbf{f}(t_k+h,\mathbf{x}_k+h\mathbf{k}_3),\\ \mathbf{x}_{k+1} &= \mathbf{x}_k + \tfrac{h}{6}\left(\mathbf{k}_1+2\mathbf{k}_2+2\mathbf{k}_3+\mathbf{k}_4\right). \end{aligned} \]

A (nontrivial) Taylor-series matching argument shows RK4 has local truncation error \( \mathcal{O}(h^5) \) and global error \( \mathcal{O}(h^4) \), making it substantially more accurate than Euler at the same step size.

For LTI systems, a particularly important simulation mode assumes the input is held constant over each interval (a standard digital implementation assumption): \( \mathbf{u}(t)=\mathbf{u}_k \) for \( t\in[t_k,t_{k+1}) \).

3.1 Exact discrete-time recursion

Starting from the exact solution and substituting \( t=t_{k+1}=t_k+h \) with ZOH input:

\[ \mathbf{x}(t_{k+1}) = e^{\mathbf{A}h}\mathbf{x}(t_k) + \int_{0}^{h} e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \; \mathbf{u}_k. \]

Define \( \mathbf{A}_d=e^{\mathbf{A}h} \) and \( \mathbf{B}_d=\int_{0}^{h} e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \). Then the exact ZOH update is

\[ \mathbf{x}_{k+1} = \mathbf{A}_d \mathbf{x}_k + \mathbf{B}_d \mathbf{u}_k. \]

3.2 Closed-form formula when \( \mathbf{A} \) is invertible

If \( \mathbf{A} \) is nonsingular, we can derive a compact expression for \( \mathbf{B}_d \). Since \( \frac{d}{d\tau}e^{\mathbf{A}\tau}=\mathbf{A}e^{\mathbf{A}\tau} \), we have

\[ \mathbf{A}\int_{0}^{h} e^{\mathbf{A}\tau}\,d\tau = \int_{0}^{h} \mathbf{A}e^{\mathbf{A}\tau}\,d\tau = \int_{0}^{h} \frac{d}{d\tau}\left(e^{\mathbf{A}\tau}\right)\,d\tau = e^{\mathbf{A}h}-\mathbf{I}. \]

Multiplying by \( \mathbf{A}^{-1} \) gives

\[ \int_{0}^{h} e^{\mathbf{A}\tau}\,d\tau = \mathbf{A}^{-1}(e^{\mathbf{A}h}-\mathbf{I}), \quad \Rightarrow \quad \mathbf{B}_d = \mathbf{A}^{-1}(\mathbf{A}_d-\mathbf{I})\mathbf{B}. \]

3.3 Van Loan augmentation (robust even when \( \mathbf{A} \) is singular)

A numerically robust method computes \( (\mathbf{A}_d,\mathbf{B}_d) \) using one matrix exponential of an augmented matrix:

\[ \exp\!\left( \begin{bmatrix} \mathbf{A} & \mathbf{B}\\ \mathbf{0} & \mathbf{0} \end{bmatrix} h \right) = \begin{bmatrix} \mathbf{A}_d & \mathbf{B}_d\\ \mathbf{0} & \mathbf{I} \end{bmatrix}. \]

Proof sketch (block ODE argument): define \( \mathbf{Z}(t)=\exp(\mathbf{M}t) \) where \( \mathbf{M}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\\\mathbf{0}&\mathbf{0}\end{bmatrix} \). Then \( \dot{\mathbf{Z}}(t)=\mathbf{M}\mathbf{Z}(t) \) with \( \mathbf{Z}(0)=\mathbf{I} \). Partition \( \mathbf{Z}(t)=\begin{bmatrix}\mathbf{Z}_{11}(t)&\mathbf{Z}_{12}(t)\\\mathbf{0}&\mathbf{I}\end{bmatrix} \). The block dynamics satisfy \( \dot{\mathbf{Z}}_{11}=\mathbf{A}\mathbf{Z}_{11} \Rightarrow \mathbf{Z}_{11}(h)=e^{\mathbf{A}h} \), and \( \dot{\mathbf{Z}}_{12}=\mathbf{A}\mathbf{Z}_{12}+\mathbf{B} \Rightarrow \mathbf{Z}_{12}(h)=\int_0^h e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \).

Numerical instability can occur even when the continuous-time LTI system is stable. The standard analysis begins with the scalar test equation \( \dot{x}=\lambda x \).

4.1 Forward Euler stability region (exact derivation)

Euler yields \( x_{k+1}=(1+h\lambda)x_k \). The numerical trajectory decays only if

\[ |1+h\lambda| < 1. \]

This is a disk in the complex plane centered at \( -1 \) with radius \( 1 \). For real negative \( \lambda \), the condition becomes

\[ -2 < h\lambda < 0 \quad \Rightarrow \quad 0 < h < \frac{2}{|\lambda|}. \]

4.2 Extension to LTI systems

For \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) and Euler:

\[ \mathbf{x}_{k+1} = (\mathbf{I}+h\mathbf{A})\mathbf{x}_k. \]

A conservative stability requirement is \( \rho(\mathbf{I}+h\mathbf{A}) < 1 \), where \( \rho(\cdot) \) is spectral radius. In practice, choose \( h \) relative to the fastest time constants (largest magnitude negative real parts of eigenvalues). RK4 typically allows larger stable steps than Euler for the same system, but still has a bounded stability region (unlike implicit methods, which are beyond this chapter’s scope).

The following flow captures the essential steps for numerically simulating \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u}(t) \) either by direct ODE stepping (Euler/RK) or by exact ZOH discretization.

flowchart TD
  S["Start: choose A, B, input u(t), x0, t0, tf"] --> G["Choose step h and time grid t[k]"]
  G --> M["Choose method: Euler / RK4 / Adaptive RK / Exact ZOH"]
  M -->|Euler/RK| L1["Loop k=0..N-1: compute f(t[k], x[k])"]
  L1 --> U1["Update x[k+1] using chosen formula"]
  M -->|Exact ZOH| D1["Compute Ad=expm(A*h) and Bd via augmented expm"]
  D1 --> L2["Loop k=0..N-1: sample u[k]=u(t[k])"]
  L2 --> U2["Update x[k+1] = Ad*x[k] + Bd*u[k]"]
  U1 --> Y["(optional) compute y(t)=C*x(t)+D*u(t)"]
  U2 --> Y
  Y --> O["Store/plot trajectories; compute errors if reference is available"]
  O --> E["End"]
        

The code below simulates a stable 2-state LTI example with a smooth input: \( \mathbf{u}(t)=\sin(2t) \). Each implementation includes: (i) fixed-step RK4 from scratch, and (ii) exact ZOH discretization via Van Loan augmentation.

6.1 Python (NumPy/SciPy; optional python-control)

File: Chapter7_Lesson5.py

"""
Chapter7_Lesson5.py
Modern Control — Chapter 7, Lesson 5: Numerical Simulation of State Equations

This script demonstrates:
1) Continuous-time simulation of x_dot = A x + B u(t) using:
   - fixed-step RK4 (implemented from scratch)
   - SciPy's adaptive RK45 (solve_ivp)
2) Exact discrete-time simulation under Zero-Order Hold (ZOH) using:
   x_{k+1} = A_d x_k + B_d u_k
   where A_d = expm(A h), B_d = integral_0^h expm(A tau) B d tau
   computed via Van Loan's augmented-matrix exponential.

Dependencies:
  numpy, scipy

Optional:
  control (python-control) for state-space helper functions.
"""

import numpy as np
from scipy.integrate import solve_ivp
from scipy.linalg import expm

def u_of_t(t: float) -> float:
    # Smooth input; for ZOH discretization we sample u_k = u(t_k)
    return float(np.sin(2.0 * t))

def f(t: float, x: np.ndarray, A: np.ndarray, B: np.ndarray) -> np.ndarray:
    return A @ x + B.flatten() * u_of_t(t)

def rk4_step(t: float, x: np.ndarray, h: float, A: np.ndarray, B: np.ndarray) -> np.ndarray:
    k1 = f(t, x, A, B)
    k2 = f(t + 0.5*h, x + 0.5*h*k1, A, B)
    k3 = f(t + 0.5*h, x + 0.5*h*k2, A, B)
    k4 = f(t + h, x + h*k3, A, B)
    return x + (h/6.0) * (k1 + 2*k2 + 2*k3 + k4)

def simulate_rk4(A: np.ndarray, B: np.ndarray, x0: np.ndarray, t0: float, tf: float, h: float):
    n_steps = int(np.floor((tf - t0)/h))
    ts = np.zeros(n_steps + 1)
    xs = np.zeros((n_steps + 1, x0.size))
    ts[0] = t0
    xs[0] = x0
    t = t0
    x = x0.copy()
    for k in range(n_steps):
        x = rk4_step(t, x, h, A, B)
        t = t + h
        ts[k+1] = t
        xs[k+1] = x
    return ts, xs

def van_loan_discretization(A: np.ndarray, B: np.ndarray, h: float):
    """
    Computes (A_d, B_d) for ZOH using Van Loan's method:

      expm( [A  B; 0  0] h ) = [A_d  B_d; 0  I]

    For a continuous-time system x_dot = A x + B u with u constant over [k h, (k+1) h).
    """
    n = A.shape[0]
    m = B.shape[1]
    M = np.zeros((n + m, n + m))
    M[:n, :n] = A
    M[:n, n:] = B
    # bottom-right already zeros
    E = expm(M * h)
    Ad = E[:n, :n]
    Bd = E[:n, n:]
    return Ad, Bd

def simulate_zoh_exact(A: np.ndarray, B: np.ndarray, x0: np.ndarray, t0: float, tf: float, h: float):
    Ad, Bd = van_loan_discretization(A, B, h)
    n_steps = int(np.floor((tf - t0)/h))
    ts = np.zeros(n_steps + 1)
    xs = np.zeros((n_steps + 1, x0.size))
    ts[0] = t0
    xs[0] = x0
    x = x0.copy()
    for k in range(n_steps):
        t = t0 + k*h
        uk = np.array([[u_of_t(t)]])  # ZOH sample
        x = Ad @ x + (Bd @ uk).flatten()
        ts[k+1] = t + h
        xs[k+1] = x
    return ts, xs, Ad, Bd

def simulate_scipy_rk45(A: np.ndarray, B: np.ndarray, x0: np.ndarray, t0: float, tf: float, h_out: float):
    sol = solve_ivp(lambda t, x: f(t, x, A, B), (t0, tf), x0, method="RK45", rtol=1e-9, atol=1e-12)
    # sample on a uniform grid for comparison
    ts = np.arange(t0, tf + 1e-12, h_out)
    xs = np.vstack([np.interp(ts, sol.t, sol.y[i, :]) for i in range(sol.y.shape[0])]).T
    return ts, xs

def main():
    # Example: stable 2-state system
    A = np.array([[0.0, 1.0],
                  [-2.0, -3.0]])
    B = np.array([[0.0],
                  [1.0]])
    x0 = np.array([1.0, 0.0])
    t0, tf = 0.0, 10.0
    h = 0.01

    ts_rk4, xs_rk4 = simulate_rk4(A, B, x0, t0, tf, h)
    ts_zoh, xs_zoh, Ad, Bd = simulate_zoh_exact(A, B, x0, t0, tf, h)
    ts_ref, xs_ref = simulate_scipy_rk45(A, B, x0, t0, tf, h_out=h)

    # Compute errors against reference at grid points
    err_rk4 = np.linalg.norm(xs_rk4 - xs_ref, axis=1)
    err_zoh = np.linalg.norm(xs_zoh - xs_ref, axis=1)

    print("A_d (ZOH exact):\n", Ad)
    print("B_d (ZOH exact):\n", Bd)
    print("Final state RK4: ", xs_rk4[-1])
    print("Final state ZOH: ", xs_zoh[-1])
    print("Final state ref: ", xs_ref[-1])
    print("Max error RK4 vs ref: ", float(np.max(err_rk4)))
    print("Max error ZOH vs ref: ", float(np.max(err_zoh)))

    # Optional: python-control demo if installed
    try:
        import control  # type: ignore
        sysc = control.ss(A, B, np.eye(2), np.zeros((2, 1)))
        t = ts_rk4
        u = np.sin(2.0 * t)
        tout, yout, xout = control.forced_response(sysc, T=t, U=u, X0=x0, return_x=True)
        print("python-control forced_response final state:", xout[:, -1])
    except Exception as e:
        print("python-control not used (optional). Reason:", str(e))

if __name__ == "__main__":
    main()
      

6.2 C++ (Eigen; RK4 + Van Loan with matrix exponential)

File: Chapter7_Lesson5.cpp

/*
Chapter7_Lesson5.cpp
Modern Control — Chapter 7, Lesson 5: Numerical Simulation of State Equations

Demonstrates:
  - Fixed-step RK4 simulation of x_dot = A x + B u(t)
  - Exact ZOH discretization via Van Loan method using matrix exponential

Dependencies:
  - Eigen (including unsupported MatrixFunctions module)

Build example (Linux/macOS):
  g++ -O2 -std=c++17 Chapter7_Lesson5.cpp -I /path/to/eigen -o sim

Note: Eigen's matrix exponential is in unsupported/Eigen/MatrixFunctions.
*/

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>

static double u_of_t(double t) {
    return std::sin(2.0 * t);
}

static Eigen::VectorXd f(double t, const Eigen::VectorXd& x,
                         const Eigen::MatrixXd& A, const Eigen::MatrixXd& B) {
    return A * x + B * Eigen::VectorXd::Constant(1, u_of_t(t));
}

static Eigen::VectorXd rk4_step(double t, const Eigen::VectorXd& x, double h,
                                const Eigen::MatrixXd& A, const Eigen::MatrixXd& B) {
    Eigen::VectorXd k1 = f(t, x, A, B);
    Eigen::VectorXd k2 = f(t + 0.5*h, x + 0.5*h*k1, A, B);
    Eigen::VectorXd k3 = f(t + 0.5*h, x + 0.5*h*k2, A, B);
    Eigen::VectorXd k4 = f(t + h, x + h*k3, A, B);
    return x + (h/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
}

static void van_loan_discretization(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, double h,
                                    Eigen::MatrixXd& Ad, Eigen::MatrixXd& Bd) {
    // expm([A B; 0 0] h) = [Ad Bd; 0 I]
    const int n = static_cast<int>(A.rows());
    const int m = static_cast<int>(B.cols());
    Eigen::MatrixXd M = Eigen::MatrixXd::Zero(n + m, n + m);
    M.block(0, 0, n, n) = A;
    M.block(0, n, n, m) = B;

    Eigen::MatrixXd E = (M * h).exp(); // matrix exponential
    Ad = E.block(0, 0, n, n);
    Bd = E.block(0, n, n, m);
}

int main() {
    Eigen::Matrix2d A;
    A << 0.0, 1.0,
        -2.0, -3.0;

    Eigen::Vector2d x0;
    x0 << 1.0, 0.0;

    Eigen::Matrix<double, 2, 1> B;
    B << 0.0,
         1.0;

    double t0 = 0.0, tf = 10.0, h = 0.01;
    int N = static_cast<int>(std::floor((tf - t0)/h));

    // RK4 simulation
    std::vector<double> ts(N + 1);
    std::vector<Eigen::Vector2d> xs(N + 1);
    ts[0] = t0;
    xs[0] = x0;

    double t = t0;
    Eigen::Vector2d x = x0;
    for (int k = 0; k < N; ++k) {
        x = rk4_step(t, x, h, A, B);
        t += h;
        ts[k+1] = t;
        xs[k+1] = x;
    }

    // Exact ZOH discretization
    Eigen::Matrix2d Ad;
    Eigen::Matrix<double, 2, 1> Bd;
    van_loan_discretization(A, B, h, Ad, Bd);

    Eigen::Vector2d xz = x0;
    for (int k = 0; k < N; ++k) {
        double tk = t0 + k*h;
        double uk = u_of_t(tk); // ZOH sample
        xz = Ad * xz + Bd * uk;
    }

    std::cout << "A_d (ZOH exact):\n" << Ad << "\n\n";
    std::cout << "B_d (ZOH exact):\n" << Bd << "\n\n";
    std::cout << "Final state RK4: " << xs.back().transpose() << "\n";
    std::cout << "Final state ZOH: " << xz.transpose() << "\n";
    return 0;
}
      

6.3 Java (Apache Commons Math; RK4 + compact matrix exponential)

File: Chapter7_Lesson5.java

/*
Chapter7_Lesson5.java
Modern Control — Chapter 7, Lesson 5: Numerical Simulation of State Equations

Demonstrates:
  - Fixed-step RK4 simulation of x_dot = A x + B u(t)
  - Exact ZOH discretization via Van Loan method, requiring matrix exponential.

This file includes a compact scaling-and-squaring + Padé(13) matrix exponential
implementation (sufficient for small/medium matrices in teaching contexts).

Dependencies:
  - Apache Commons Math 3.x (linear algebra)

Compile (example):
  javac -cp commons-math3-3.6.1.jar Chapter7_Lesson5.java
Run:
  java -cp .:commons-math3-3.6.1.jar Chapter7_Lesson5
*/

import org.apache.commons.math3.linear.*;
import java.util.Arrays;

public class Chapter7_Lesson5 {

    static double uOfT(double t) {
        return Math.sin(2.0 * t);
    }

    static RealVector f(double t, RealVector x, RealMatrix A, RealMatrix B) {
        // x_dot = A x + B u(t) ; here u is scalar
        return A.operate(x).add(B.getColumnVector(0).mapMultiply(uOfT(t)));
    }

    static RealVector rk4Step(double t, RealVector x, double h, RealMatrix A, RealMatrix B) {
        RealVector k1 = f(t, x, A, B);
        RealVector k2 = f(t + 0.5*h, x.add(k1.mapMultiply(0.5*h)), A, B);
        RealVector k3 = f(t + 0.5*h, x.add(k2.mapMultiply(0.5*h)), A, B);
        RealVector k4 = f(t + h, x.add(k3.mapMultiply(h)), A, B);
        return x.add(k1.add(k2.mapMultiply(2.0)).add(k3.mapMultiply(2.0)).add(k4).mapMultiply(h/6.0));
    }

    // ---------- Matrix exponential: scaling/squaring + Padé(13) ----------
    static RealMatrix expm(RealMatrix A) {
        int n = A.getRowDimension();
        double norm1 = oneNorm(A);
        double theta13 = 5.371920351148152; // Padé(13) threshold
        int s = Math.max(0, (int)Math.ceil(Math.log(norm1/theta13)/Math.log(2.0)));

        RealMatrix As = A.scalarMultiply(1.0 / Math.pow(2.0, s));
        RealMatrix E = pade13(As);

        for (int i = 0; i < s; i++) {
            E = E.multiply(E);
        }
        return E;
    }

    static double oneNorm(RealMatrix A) {
        int n = A.getColumnDimension();
        double max = 0.0;
        for (int j = 0; j < n; j++) {
            double colSum = 0.0;
            for (int i = 0; i < A.getRowDimension(); i++) {
                colSum += Math.abs(A.getEntry(i, j));
            }
            if (colSum > max) max = colSum;
        }
        return max;
    }

    static RealMatrix pade13(RealMatrix A) {
        double[] b = new double[] {
            64764752532480000.0,
            32382376266240000.0,
            7771770303897600.0,
            1187353796428800.0,
            129060195264000.0,
            10559470521600.0,
            670442572800.0,
            33522128640.0,
            1323241920.0,
            40840800.0,
            960960.0,
            16380.0,
            182.0,
            1.0
        };

        int n = A.getRowDimension();
        RealMatrix I = MatrixUtils.createRealIdentityMatrix(n);

        RealMatrix A2 = A.multiply(A);
        RealMatrix A4 = A2.multiply(A2);
        RealMatrix A6 = A4.multiply(A2);
        RealMatrix A8  = A4.multiply(A4);
        RealMatrix A10 = A6.multiply(A4);
        RealMatrix A12 = A6.multiply(A6);
        RealMatrix A3  = A2.multiply(A);
        RealMatrix A5  = A4.multiply(A);
        RealMatrix A7  = A6.multiply(A);
        RealMatrix A9  = A8.multiply(A);
        RealMatrix A11 = A10.multiply(A);
        RealMatrix A13 = A12.multiply(A);

        RealMatrix U = A13.scalarMultiply(b[13]).add(A11.scalarMultiply(b[11])).add(A9.scalarMultiply(b[9]))
                .add(A7.scalarMultiply(b[7])).add(A5.scalarMultiply(b[5])).add(A3.scalarMultiply(b[3])).add(A.scalarMultiply(b[1]));

        RealMatrix V = A12.scalarMultiply(b[12]).add(A10.scalarMultiply(b[10])).add(A8.scalarMultiply(b[8]))
                .add(A6.scalarMultiply(b[6])).add(A4.scalarMultiply(b[4])).add(A2.scalarMultiply(b[2])).add(I.scalarMultiply(b[0]));

        RealMatrix P = V.add(U);
        RealMatrix Q = V.subtract(U);

        DecompositionSolver solver = new LUDecomposition(Q).getSolver();
        return solver.solve(P);
    }

    static class Discretization {
        RealMatrix Ad;
        RealMatrix Bd;
        Discretization(RealMatrix Ad, RealMatrix Bd) { this.Ad = Ad; this.Bd = Bd; }
    }

    static Discretization vanLoanDiscretization(RealMatrix A, RealMatrix B, double h) {
        int n = A.getRowDimension();
        int m = B.getColumnDimension();
        RealMatrix M = MatrixUtils.createRealMatrix(n + m, n + m);
        M.setSubMatrix(A.getData(), 0, 0);
        M.setSubMatrix(B.getData(), 0, n);

        RealMatrix E = expm(M.scalarMultiply(h));
        RealMatrix Ad = E.getSubMatrix(0, n-1, 0, n-1);
        RealMatrix Bd = E.getSubMatrix(0, n-1, n, n+m-1);
        return new Discretization(Ad, Bd);
    }

    public static void main(String[] args) {
        RealMatrix A = MatrixUtils.createRealMatrix(new double[][]{
                {0.0, 1.0},
                {-2.0, -3.0}
        });
        RealMatrix B = MatrixUtils.createRealMatrix(new double[][]{
                {0.0},
                {1.0}
        });
        RealVector x0 = MatrixUtils.createRealVector(new double[]{1.0, 0.0});

        double t0 = 0.0, tf = 10.0, h = 0.01;
        int N = (int)Math.floor((tf - t0)/h);

        RealVector x = x0.copy();
        double t = t0;
        for (int k = 0; k < N; k++) {
            x = rk4Step(t, x, h, A, B);
            t += h;
        }

        Discretization disc = vanLoanDiscretization(A, B, h);
        RealVector xz = x0.copy();
        for (int k = 0; k < N; k++) {
            double tk = t0 + k*h;
            double uk = uOfT(tk);
            xz = disc.Ad.operate(xz).add(disc.Bd.getColumnVector(0).mapMultiply(uk));
        }

        System.out.println("Final state RK4 : " + Arrays.toString(x.toArray()));
        System.out.println("Final state ZOH : " + Arrays.toString(xz.toArray()));
    }
}
      

6.4 MATLAB/Simulink (ode45, RK4, expm, programmatic Simulink)

File: Chapter7_Lesson5.m

% Chapter7_Lesson5.m
% Modern Control — Chapter 7, Lesson 5: Numerical Simulation of State Equations
%
% Demonstrates:
%  1) Continuous-time simulation using ode45 (adaptive RK)
%  2) Fixed-step RK4 (implemented explicitly)
%  3) Exact ZOH discretization via Van Loan augmented exponential.
%
% Run:
%   Chapter7_Lesson5

clear; clc;

A = [0 1; -2 -3];
B = [0; 1];
x0 = [1; 0];

u = @(t) sin(2*t);
f = @(t,x) A*x + B*u(t);

t0 = 0; tf = 10; h = 0.01;
tgrid = (t0:h:tf)';

%% 1) ode45 reference
opts = odeset('RelTol',1e-9,'AbsTol',1e-12);
[ts_ref, xs_ref] = ode45(f, [t0 tf], x0, opts);
x_ref = interp1(ts_ref, xs_ref, tgrid);

%% 2) Fixed-step RK4
N = numel(tgrid);
x_rk4 = zeros(N, numel(x0));
x_rk4(1,:) = x0.';
x = x0;
for k = 1:N-1
    t = tgrid(k);
    k1 = f(t, x);
    k2 = f(t + 0.5*h, x + 0.5*h*k1);
    k3 = f(t + 0.5*h, x + 0.5*h*k2);
    k4 = f(t + h, x + h*k3);
    x = x + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
    x_rk4(k+1,:) = x.';
end

%% 3) Exact ZOH discretization (Van Loan)
[Ad, Bd] = van_loan_discretization(A, B, h);

x_zoh = zeros(N, numel(x0));
x_zoh(1,:) = x0.';
x = x0;
for k = 1:N-1
    t = tgrid(k);
    uk = u(t); % ZOH sample
    x = Ad*x + Bd*uk;
    x_zoh(k+1,:) = x.';
end

%% Error comparison
err_rk4 = vecnorm(x_rk4 - x_ref, 2, 2);
err_zoh = vecnorm(x_zoh - x_ref, 2, 2);

fprintf('Max error RK4 vs ode45: %.3e\n', max(err_rk4));
fprintf('Max error ZOH vs ode45: %.3e\n', max(err_zoh));

%% Optional: Simulink programmatic model (no images)
try
    mdl = 'Chapter7_Lesson5_Simulink';
    if bdIsLoaded(mdl); close_system(mdl, 0); end
    new_system(mdl); open_system(mdl);

    add_block('simulink/Sources/Sine Wave', [mdl '/u']);
    set_param([mdl '/u'], 'Amplitude', '1', 'Frequency', '2'); % sin(2 t)

    add_block('simulink/Continuous/State-Space', [mdl '/SS']);
    set_param([mdl '/SS'], 'A', mat2str(A), 'B', mat2str(B), 'C', mat2str(eye(2)), 'D', 'zeros(2,1)');

    add_block('simulink/Sinks/To Workspace', [mdl '/x_out']);
    set_param([mdl '/x_out'], 'VariableName', 'x_sim', 'SaveFormat', 'Array');

    add_line(mdl, 'u/1', 'SS/1');
    add_line(mdl, 'SS/1', 'x_out/1');

    set_param(mdl, 'StopTime', num2str(tf));
    sim(mdl);
    disp('Simulink model created and simulated: Chapter7_Lesson5_Simulink');
catch ME
    disp('Simulink step skipped (Simulink not available or error):');
    disp(ME.message);
end

%% ---- Local function: Van Loan discretization ----
function [Ad, Bd] = van_loan_discretization(A, B, h)
    % expm([A B; 0 0] h) = [Ad Bd; 0 I]
    n = size(A,1);
    m = size(B,2);
    M = zeros(n+m, n+m);
    M(1:n,1:n) = A;
    M(1:n,n+1:n+m) = B;
    E = expm(M*h);
    Ad = E(1:n,1:n);
    Bd = E(1:n,n+1:n+m);
end
      

6.5 Wolfram Mathematica (NDSolve, RK4 from scratch, MatrixExp)

File: Chapter7_Lesson5.nb

(* ::Package:: *)

(* Chapter7_Lesson5.nb
   Modern Control — Chapter 7, Lesson 5: Numerical Simulation of State Equations

   This notebook (in text form) demonstrates:
     - Continuous-time simulation with NDSolve
     - Fixed-step RK4 implemented from scratch
     - Exact ZOH discretization via MatrixExp and Van Loan augmentation
*)

Notebook[{
  Cell["Chapter 7 — Lesson 5: Numerical Simulation of State Equations", "Title"],

  Cell["Define an LTI state equation x'(t) = A x(t) + B u(t).", "Text"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"A", "=", RowBox[{"{", RowBox[{RowBox[{"{", RowBox[{"0", ",", "1"}], "}"}], ",",
                                        RowBox[{"{", RowBox[{"-2", ",", RowBox[{"-", "3"}]}], "}"}]}], "}"}]}], ";",
      RowBox[{"B", "=", RowBox[{"{", RowBox[{RowBox[{"{", "0", "}"}], ",", RowBox[{"{", "1", "}"}]}], "}"}]}], ";",
      RowBox[{"u", "[", "t_", "]"}], ":=", RowBox[{"Sin", "[", RowBox[{"2", " ", "t"}], "]"}], ";",
      RowBox[{"x0", "=", RowBox[{"{", RowBox[{"1", ",", "0"}], "}"}]}], ";",
      RowBox[{"t0", "=", "0"}], ";", RowBox[{"tf", "=", "10"}], ";", RowBox[{"h", "=", "0.01"}], ";"
    }]], "Input"],

  Cell["1) Reference continuous-time solution using NDSolve.", "Section"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"sol", "=", RowBox[{"NDSolveValue", "[",
        RowBox[{
          RowBox[{"{",
            RowBox[{
              RowBox[{
                RowBox[{"x", "'"}], "[", "t", "]"}], "\[Equal]",
                RowBox[{
                  RowBox[{"A", ".", RowBox[{"x", "[", "t", "]"}]}], "+",
                  RowBox[{"Flatten", "[", RowBox[{"B", " ", RowBox[{"u", "[", "t", "]"}]}], "]"}]
                }
              }],
              ",",
              RowBox[{
                RowBox[{"x", "[", "t0", "]"}], "\[Equal]", "x0"
              }]
            }], "}"}], ",",
          "x", ",", RowBox[{"{", RowBox[{"t", ",", "t0", ",", "tf"}], "}"}]
        }], "]"}]}], ";"
    }]], "Input"],

  Cell["2) Fixed-step RK4 implemented from scratch.", "Section"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"f", "[", RowBox[{"t_", ",", "x_"}], "]"}], ":=",
        RowBox[{
          RowBox[{"A", ".", "x"}], "+",
          RowBox[{"Flatten", "[", RowBox[{"B", " ", RowBox[{"u", "[", "t", "]"}]}], "]"}]
        }], ";",
      RowBox[{"rk4Step", "[", RowBox[{"t_", ",", "x_", ",", "h_"}], "]"}], ":=",
        RowBox[{"Module", "[",
          RowBox[{
            RowBox[{"{", RowBox[{"k1", ",", "k2", ",", "k3", ",", "k4"}], "}"}], ",",
            RowBox[{
              RowBox[{"k1", "=", RowBox[{"f", "[", RowBox[{"t", ",", "x"}], "]"}]}], ";",
              RowBox[{"k2", "=", RowBox[{"f", "[", RowBox[{"t", "+", RowBox[{"h", "/", "2"}], ",",
                                RowBox[{"x", "+", RowBox[{"h", " ", RowBox[{"k1", "/", "2"}]}]}]}], "]"}]}], ";",
              RowBox[{"k3", "=", RowBox[{"f", "[", RowBox[{"t", "+", RowBox[{"h", "/", "2"}], ",",
                                RowBox[{"x", "+", RowBox[{"h", " ", RowBox[{"k2", "/", "2"}]}]}]}], "]"}]}], ";",
              RowBox[{"k4", "=", RowBox[{"f", "[", RowBox[{"t", "+", "h", ",",
                                RowBox[{"x", "+", RowBox[{"h", " ", "k3"}]}]}], "]"}]}], ";",
              RowBox[{"x", "+", RowBox[{"h", " ", RowBox[{"(", RowBox[{"k1", "+", RowBox[{"2", " ", "k2"}], "+",
                               RowBox[{"2", " ", "k3"}], "+", "k4"}], ")"}], "/", "6"}]}]
            }
          }], "]"}], ";"
    }]], "Input"],

  Cell["3) Exact ZOH discretization (Van Loan).", "Section"],
  Cell[BoxData[
    RowBox[{
      RowBox[{"n", "=", RowBox[{"Length", "[", "A", "]"}]}], ";",
      RowBox[{"m", "=", RowBox[{"Dimensions", "[", "B", "]"}], "[", RowBox[{"[", "2", "]"}], "]"}], ";",
      RowBox[{"M", "=", RowBox[{"ArrayFlatten", "[",
        RowBox[{"{",
          RowBox[{
            RowBox[{"{", RowBox[{"A", ",", "B"}], "}"}], ",",
            RowBox[{"{", RowBox[{RowBox[{"ConstantArray", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"m", ",", "n"}], "}"}]}], "]"}], ",",
                             RowBox[{"ConstantArray", "[", RowBox[{"0", ",", RowBox[{"{", RowBox[{"m", ",", "m"}], "}"}]}], "]"}]}], "}"}]
          }], "}"}], "]"}]}], ";",
      RowBox[{"E", "=", RowBox[{"MatrixExp", "[", RowBox[{"M", " ", "h"}], "]"}]}], ";",
      RowBox[{"Ad", "=", RowBox[{"E", "[", RowBox[{"[", RowBox[{RowBox[{"1", ";;", "n"}], ",", RowBox[{"1", ";;", "n"}]}], "]"}], "]"}]}], ";",
      RowBox[{"Bd", "=", RowBox[{"E", "[", RowBox[{"[", RowBox[{RowBox[{"1", ";;", "n"}], ",", RowBox[{"n", "+", "1", ";;", "n", "+", "m"}]}], "]"}], "]"}]}], ";"
    }]], "Input"]
}]
      

Problem 1 (Euler local truncation error): Let \( \dot{\mathbf{x}}=\mathbf{f}(t,\mathbf{x}) \) with \( \mathbf{x}(t) \) twice continuously differentiable. Show that Euler’s local truncation error satisfies \( \boldsymbol{\tau}_{k+1}=\mathcal{O}(h^2) \).

Solution: By Taylor’s theorem,

\[ \mathbf{x}(t_{k+1}) = \mathbf{x}(t_k) + h\dot{\mathbf{x}}(t_k) + \frac{h^2}{2}\ddot{\mathbf{x}}(\xi_k) \]

for some \( \xi_k\in(t_k,t_{k+1}) \). Using \( \dot{\mathbf{x}}(t_k)=\mathbf{f}(t_k,\mathbf{x}(t_k)) \) and the definition of \( \boldsymbol{\tau}_{k+1} \) yields

\[ \boldsymbol{\tau}_{k+1} = \mathbf{x}(t_{k+1})-\Big(\mathbf{x}(t_k)+h\mathbf{f}(t_k,\mathbf{x}(t_k))\Big) = \frac{h^2}{2}\ddot{\mathbf{x}}(\xi_k) = \mathcal{O}(h^2). \]

Problem 2 (Exact ZOH discretization): Assume \( \mathbf{u}(t)=\mathbf{u}_k \) for \( t\in[t_k,t_{k+1}) \). Starting from the exact continuous-time solution, derive \( \mathbf{x}_{k+1}=\mathbf{A}_d\mathbf{x}_k+\mathbf{B}_d\mathbf{u}_k \) with \( \mathbf{A}_d=e^{\mathbf{A}h} \) and \( \mathbf{B}_d=\int_0^h e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \).

Solution: From \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_k)}\mathbf{x}(t_k)+\int_{t_k}^{t}e^{\mathbf{A}(t-\tau)}\mathbf{B}\mathbf{u}(\tau)\,d\tau \), set \( t=t_k+h \) and \( \mathbf{u}(\tau)=\mathbf{u}_k \). Change variables \( \sigma=t_k+h-\tau \), obtaining

\[ \mathbf{x}(t_{k+1}) = e^{\mathbf{A}h}\mathbf{x}(t_k) + \left(\int_0^h e^{\mathbf{A}\sigma}\mathbf{B}\,d\sigma\right)\mathbf{u}_k, \]

which is the claimed discrete-time recursion.

Problem 3 (Euler stability step bound for a scalar mode): For \( \dot{x}=\lambda x \) with \( \lambda < 0 \), find the range of \( h \) such that Euler does not diverge.

Solution: Euler gives \( x_{k+1}=(1+h\lambda)x_k \). Decay requires \( |1+h\lambda| < 1 \), i.e.

\[ -1 < 1+h\lambda < 1 \quad \Rightarrow \quad -2 < h\lambda < 0 \quad \Rightarrow \quad 0 < h < \frac{2}{|\lambda|}. \]

Problem 4 (Van Loan identity): Let \( \mathbf{M}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\\\mathbf{0}&\mathbf{0}\end{bmatrix} \). Show that the top-right block of \( e^{\mathbf{M}h} \) equals \( \int_0^h e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \).

Solution: Let \( \mathbf{Z}(t)=e^{\mathbf{M}t} \). Then \( \dot{\mathbf{Z}}=\mathbf{M}\mathbf{Z} \) with \( \mathbf{Z}(0)=\mathbf{I} \). Partition \( \mathbf{Z}(t)=\begin{bmatrix}\mathbf{Z}_{11}(t)&\mathbf{Z}_{12}(t)\\\mathbf{0}&\mathbf{I}\end{bmatrix} \). Block multiplication yields \( \dot{\mathbf{Z}}_{11}=\mathbf{A}\mathbf{Z}_{11} \) and \( \dot{\mathbf{Z}}_{12}=\mathbf{A}\mathbf{Z}_{12}+\mathbf{B} \) with \( \mathbf{Z}_{12}(0)=\mathbf{0} \). Variation of constants gives

\[ \mathbf{Z}_{12}(h)=\int_0^h e^{\mathbf{A}(h-\tau)}\mathbf{B}\,d\tau = \int_0^h e^{\mathbf{A}\sigma}\mathbf{B}\,d\sigma, \]

after the substitution \( \sigma=h-\tau \).

We formulated numerical simulation of LTI state equations as ODE integration, derived Euler and RK4 update laws, and proved Euler’s local error \( \mathcal{O}(h^2) \) and global error \( \mathcal{O}(h) \). We then exploited LTI structure to derive the exact ZOH discrete-time recursion with \( \mathbf{A}_d=e^{\mathbf{A}h} \) and \( \mathbf{B}_d=\int_0^h e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \), including the robust Van Loan augmented-exponential computation. The accompanying multi-language implementations illustrate both approaches in practice.