In real digital control and simulation pipelines, the computer updates the control input only at sampling instants \( t_k = kT_s \), where \( T_s \) is the sampling period. Between two updates, the actuator typically holds the last command constant. This is the zero-order hold assumption.

If the continuous-time plant is modeled by \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \) and \( \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t) \), then under ZOH the sampled system is not an approximation if we discretize correctly: it has an exact discrete-time representation at the sample instants.

This is crucial because the discrete model used for analysis/design should reflect the true sampled behavior of the continuous system, not merely a numerical integration approximation.

flowchart TD
  A["Continuous plant: xdot = A x + B u"] --> B["Sampler with period Ts"]
  B --> C["Controller computes u[k]"]
  C --> D["Hold block keeps u(t) constant on each interval"]
  D --> E["State evolves continuously during kTs to (k+1)Ts"]
  E --> F["Exact sampled recursion: x[k+1] = Ad x[k] + Bd u[k]"]
        

Let \( \mathbf{x}_k \equiv \mathbf{x}(kT_s) \) and \( \mathbf{u}_k \equiv \mathbf{u}(kT_s) \). Under a zero-order hold,

\[ \mathbf{u}(t)=\mathbf{u}_k,\qquad kT_s \le t < (k+1)T_s. \]

Over one sampling interval, the continuous system becomes a linear ODE with constant input:

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}_k,\qquad t\in[kT_s,(k+1)T_s). \]

Using the variation-of-constants formula and shifting the local interval variable \( \sigma = t-kT_s \), we obtain:

\[ \mathbf{x}((k+1)T_s)=e^{\mathbf{A}T_s}\mathbf{x}(kT_s) +\int_{0}^{T_s} e^{\mathbf{A}(T_s-\sigma)}\mathbf{B}\,d\sigma\;\mathbf{u}_k. \]

By the change of variable \( \tau=T_s-\sigma \), the integral can be written in the standard form:

\[ \mathbf{x}_{k+1}=\mathbf{A}_d\mathbf{x}_k+\mathbf{B}_d\mathbf{u}_k, \qquad \mathbf{A}_d=e^{\mathbf{A}T_s}, \qquad \mathbf{B}_d=\int_{0}^{T_s} e^{\mathbf{A}\tau}\mathbf{B}\,d\tau. \]

The output equation at the sample instants is \( \mathbf{y}_k=\mathbf{C}\mathbf{x}_k+\mathbf{D}\mathbf{u}_k \).

For completeness, we prove the exact result. Consider the local-time system on one sampling interval:

\[ \frac{d}{d\tau}\mathbf{x}(\tau)=\mathbf{A}\mathbf{x}(\tau)+\mathbf{B}\mathbf{u}_k, \qquad \mathbf{x}(0)=\mathbf{x}_k, \qquad 0\le \tau < T_s. \]

Multiply by \( e^{-\mathbf{A}\tau} \) from the left and use the product rule:

\[ \frac{d}{d\tau}\left(e^{-\mathbf{A}\tau}\mathbf{x}(\tau)\right)=e^{-\mathbf{A}\tau}\mathbf{B}\mathbf{u}_k. \]

Integrating from \( 0 \) to \( T_s \) gives

\[ e^{-\mathbf{A}T_s}\mathbf{x}(T_s)-\mathbf{x}(0)=\int_0^{T_s} e^{-\mathbf{A}\tau}\mathbf{B}\,d\tau\;\mathbf{u}_k. \]

Multiplying both sides by \( e^{\mathbf{A}T_s} \) yields

\[ \mathbf{x}(T_s)=e^{\mathbf{A}T_s}\mathbf{x}(0)+\int_0^{T_s} e^{\mathbf{A}(T_s-\tau)}\mathbf{B}\,d\tau\;\mathbf{u}_k, \]

which is exactly the one-step sampled recursion after reindexing by \( k \). Therefore, the ZOH discretization is mathematically exact at sample times.

The matrix \( \mathbf{A}_d=e^{\mathbf{A}T_s} \) is direct. The challenging term is usually \( \mathbf{B}_d \).

If \( \mathbf{A} \) is nonsingular, then

\[ \mathbf{B}_d =\left(\int_0^{T_s} e^{\mathbf{A}\tau}\,d\tau\right)\mathbf{B} =\mathbf{A}^{-1}\left(e^{\mathbf{A}T_s}-\mathbf{I}\right)\mathbf{B} =\mathbf{A}^{-1}\left(\mathbf{A}_d-\mathbf{I}\right)\mathbf{B}. \]

However, this formula fails numerically or formally when \( \mathbf{A} \) is singular. A robust method uses the augmented matrix

\[ \mathbf{M}= \begin{bmatrix} \mathbf{A} & \mathbf{B}\\ \mathbf{0} & \mathbf{0} \end{bmatrix}. \]

Then

\[ e^{\mathbf{M}T_s} = \begin{bmatrix} \mathbf{A}_d & \mathbf{B}_d\\ \mathbf{0} & \mathbf{I} \end{bmatrix}. \]

This identity follows from the block-triangular series expansion and is the preferred computational approach in scientific software.

flowchart TD
  A["Given A, B, Ts"] --> B["Build augmented matrix M = [[A,B],[0,0]]"]
  B --> C["Compute exp(M*Ts) with stable matrix-exponential routine"]
  C --> D["Extract top-left block as Ad"]
  C --> E["Extract top-right block as Bd"]
  D --> F["Use x[k+1] = Ad x[k] + Bd u[k]"]
  E --> F
        

Scalar first-order system. Let \( \dot{x}(t)=-ax(t)+bu(t) \) with \( a > 0 \). Then \( A=-a \) and \( B=b \), so

\[ A_d=e^{-aT_s}, \qquad B_d=\int_0^{T_s} e^{-a\tau} b\,d\tau =\frac{b}{a}\left(1-e^{-aT_s}\right). \]

Hence the exact discrete model is

\[ x_{k+1}=e^{-aT_s}x_k+\frac{b}{a}\left(1-e^{-aT_s}\right)u_k. \]

Second-order mass–spring–damper model. For \( m\ddot{q}+c\dot{q}+kq=bu \), choose state \( \mathbf{x}=[q\;\; \dot{q}]^T \), then

\[ \dot{\mathbf{x} }= \begin{bmatrix} 0 & 1\\ -k/m & -c/m \end{bmatrix}\mathbf{x} + \begin{bmatrix} 0\\ b/m \end{bmatrix}u. \]

The exact sampled model is obtained directly from the formulas above. In this lesson's code, we compare that exact model against forward Euler to highlight modeling error when \( T_s \) is not very small.

Since Lesson 3 introduced the z-transform, we now connect the exact state discretization to the discrete transfer function. For the exact sampled state model,

\[ \mathbf{x}_{k+1}=\mathbf{A}_d\mathbf{x}_k+\mathbf{B}_d\mathbf{u}_k,\qquad \mathbf{y}_k=\mathbf{C}\mathbf{x}_k+\mathbf{D}\mathbf{u}_k, \]

the discrete-time transfer function is

\[ \mathbf{G}_d(z)=\mathbf{C}\left(z\mathbf{I}-\mathbf{A}_d\right)^{-1}\mathbf{B}_d+\mathbf{D}. \]

For a strictly proper SISO continuous plant \( G(s) \) with ZOH at the input, the pulse-transfer relationship can also be written as a sampled convolution. A common expression is

\[ G_d(z)=\mathcal{Z}\!\left\{ g_h(kT_s) \right\},\qquad g_h(t)=\mathcal{L}^{-1}\!\left\{ G(s)\frac{1-e^{-sT_s} }{s} \right\}, \]

where \( g_h(t) \) is the hold-equivalent impulse response sequence generator. In practice, the state-space exact discretization is often numerically cleaner and generalizes immediately to MIMO systems.

A major theoretical advantage of exact discretization is that modal information is preserved at sampling instants. If \( \lambda_i(\mathbf{A}) \) is a continuous-time eigenvalue, then the corresponding discrete-time eigenvalue is

\[ \mu_i=e^{\lambda_i(\mathbf{A})T_s}. \]

Therefore,

\[ \mathrm{Re}\!\left(\lambda_i(\mathbf{A})\right) < 0 \;\Rightarrow\; |\mu_i|=|e^{\lambda_i(\mathbf{A})T_s}| < 1. \]

Thus asymptotic stability of the continuous LTI system implies Schur stability of the exact ZOH discrete model. By contrast, approximate discretization schemes (Euler, backward Euler, trapezoidal, etc.) impose their own numerical stability regions and may distort the sampled dynamics.

This distinction is important: \( \mathbf{A}_d=e^{\mathbf{A}T_s} \) is a modeling conversion, while Euler/RK are integration methods. They serve different purposes.

For forward Euler, one writes \( \mathbf{x}_{k+1}\approx (\mathbf{I}+T_s\mathbf{A})\mathbf{x}_k + T_s\mathbf{B}\mathbf{u}_k \). This is a first-order approximation to the exact ZOH model:

\[ e^{\mathbf{A}T_s} =\mathbf{I}+T_s\mathbf{A}+\frac{T_s^2}{2!}\mathbf{A}^2+\frac{T_s^3}{3!}\mathbf{A}^3+\cdots \]

\[ \mathbf{B}_d =\left(T_s\mathbf{I}+\frac{T_s^2}{2!}\mathbf{A}+\frac{T_s^3}{3!}\mathbf{A}^2+\cdots\right)\mathbf{B}. \]

Hence Euler drops all higher-order terms. For small \( T_s \), this may be acceptable. For larger sampling periods or faster dynamics, it can significantly shift poles and time response.

In control engineering workflows, exact ZOH conversion is typically used for plant discretization, while numerical integrators are used for simulation of more general systems (especially nonlinear systems from Chapter 15).

The following implementations all use the same lesson example (mass–spring–damper) and compute exact ZOH discretization. Where possible, we also compare against Euler discretization.

Python libraries: numpy, scipy.linalg.expm, matplotlib (and optionally python-control for c2d comparisons).

Code: Chapter16_Lesson4.py

# Chapter16_Lesson4.py
# Continuous–Discrete Conversions: Zero-Order Hold, Exact Discretization
# Python implementation (NumPy/SciPy/Matplotlib)

import numpy as np
from scipy.linalg import expm
import matplotlib.pyplot as plt


def exact_discretize(A: np.ndarray, B: np.ndarray, Ts: float):
    """
    Exact ZOH discretization via augmented matrix exponential.
    Works for singular or nonsingular A.
    x[k+1] = Ad x[k] + Bd u[k]
    """
    n, m = B.shape
    M = np.zeros((n + m, n + m))
    M[:n, :n] = A
    M[:n, n:] = B
    Md = expm(M * Ts)
    Ad = Md[:n, :n]
    Bd = Md[:n, n:]
    return Ad, Bd


def euler_discretize(A: np.ndarray, B: np.ndarray, Ts: float):
    """Forward-Euler (approximate) discretization for comparison."""
    n = A.shape[0]
    Ad = np.eye(n) + Ts * A
    Bd = Ts * B
    return Ad, Bd


def simulate_discrete(Ad, Bd, C, D, u_seq, x0):
    n_steps = len(u_seq)
    x = np.array(x0, dtype=float).reshape(-1, 1)
    ys = []
    xs = [x.flatten()]
    for k in range(n_steps):
        u = np.array([[u_seq[k]]], dtype=float)
        y = C @ x + D @ u
        ys.append(float(y))
        x = Ad @ x + Bd @ u
        xs.append(x.flatten())
    return np.array(xs), np.array(ys)


def main():
    # Example: mass-spring-damper (position/velocity states)
    # m xdd + c xd + k x = b u
    m, c, k, b = 1.0, 0.6, 4.0, 1.0
    A = np.array([[0.0, 1.0],
                  [-k / m, -c / m]])
    B = np.array([[0.0],
                  [b / m]])
    C = np.array([[1.0, 0.0]])
    D = np.array([[0.0]])

    Ts = 0.1  # sample period [s]

    Ad_exact, Bd_exact = exact_discretize(A, B, Ts)
    Ad_euler, Bd_euler = euler_discretize(A, B, Ts)

    print("A =\n", A)
    print("B =\n", B)
    print("\nExact ZOH discretization:")
    print("Ad_exact =\n", Ad_exact)
    print("Bd_exact =\n", Bd_exact)

    print("\nForward Euler discretization (approximate):")
    print("Ad_euler =\n", Ad_euler)
    print("Bd_euler =\n", Bd_euler)

    # Input sequence: unit step after k=5
    N = 120
    u_seq = np.zeros(N)
    u_seq[5:] = 1.0

    x0 = np.array([0.0, 0.0])

    xs_exact, y_exact = simulate_discrete(Ad_exact, Bd_exact, C, D, u_seq, x0)
    xs_euler, y_euler = simulate_discrete(Ad_euler, Bd_euler, C, D, u_seq, x0)

    t = np.arange(N) * Ts

    plt.figure(figsize=(9, 4.8))
    plt.plot(t, y_exact, label="Exact ZOH")
    plt.plot(t, y_euler, "--", label="Forward Euler")
    plt.step(t, u_seq, where="post", label="Input u[k]", alpha=0.6)
    plt.xlabel("Time [s]")
    plt.ylabel("Output / Input")
    plt.title("Continuous-to-Discrete Conversion: Exact ZOH vs Euler")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

    # Stability check: map continuous eigenvalues to discrete eigenvalues
    lam_c = np.linalg.eigvals(A)
    lam_d = np.linalg.eigvals(Ad_exact)
    print("\nContinuous eigenvalues:", lam_c)
    print("Discrete eigenvalues (exact):", lam_d)
    print("Magnitudes |lambda_d|:", np.abs(lam_d))


if __name__ == "__main__":
    main()

C++ libraries: Eigen and unsupported/Eigen/MatrixFunctions for the matrix exponential. This is the standard scientific-computing path in modern C++.

Code: Chapter16_Lesson4.cpp

// Chapter16_Lesson4.cpp
// Continuous–Discrete Conversions: Zero-Order Hold, Exact Discretization
// C++ implementation using Eigen (and unsupported MatrixFunctions for exp())

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

struct DiscreteModel {
    Eigen::MatrixXd Ad;
    Eigen::MatrixXd Bd;
};

DiscreteModel exactDiscretizeZOH(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, double Ts) {
    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 Md = (M * Ts).exp();  // matrix exponential
    DiscreteModel dm;
    dm.Ad = Md.block(0, 0, n, n);
    dm.Bd = Md.block(0, n, n, m);
    return dm;
}

DiscreteModel eulerDiscretize(const Eigen::MatrixXd& A, const Eigen::MatrixXd& B, double Ts) {
    DiscreteModel dm;
    dm.Ad = Eigen::MatrixXd::Identity(A.rows(), A.cols()) + Ts * A;
    dm.Bd = Ts * B;
    return dm;
}

int main() {
    // Mass-spring-damper example
    const double m = 1.0, c = 0.6, k = 4.0, b = 1.0;
    Eigen::MatrixXd A(2, 2);
    Eigen::MatrixXd B(2, 1);
    Eigen::MatrixXd C(1, 2);
    Eigen::MatrixXd D(1, 1);

    A << 0.0, 1.0,
        -k / m, -c / m;
    B << 0.0,
         b / m;
    C << 1.0, 0.0;
    D << 0.0;

    const double Ts = 0.1;

    DiscreteModel exact = exactDiscretizeZOH(A, B, Ts);
    DiscreteModel euler = eulerDiscretize(A, B, Ts);

    std::cout << "Ad (exact ZOH):\n" << exact.Ad << "\n\n";
    std::cout << "Bd (exact ZOH):\n" << exact.Bd << "\n\n";
    std::cout << "Ad (Euler):\n" << euler.Ad << "\n\n";
    std::cout << "Bd (Euler):\n" << euler.Bd << "\n\n";

    // Simulate with a step input (u[k] = 0 for k<5, then 1)
    const int N = 120;
    Eigen::VectorXd xExact = Eigen::VectorXd::Zero(2);
    Eigen::VectorXd xEuler = Eigen::VectorXd::Zero(2);

    std::ofstream csv("Chapter16_Lesson4_output.csv");
    csv << "k,t,u,y_exact,y_euler,x1_exact,x2_exact\n";

    for (int kstep = 0; kstep < N; ++kstep) {
        double u = (kstep >= 5) ? 1.0 : 0.0;

        double yExact = (C * xExact)(0, 0) + D(0, 0) * u;
        double yEuler = (C * xEuler)(0, 0) + D(0, 0) * u;

        csv << kstep << "," << (kstep * Ts) << "," << u << ","
            << yExact << "," << yEuler << ","
            << xExact(0) << "," << xExact(1) << "\n";

        xExact = exact.Ad * xExact + exact.Bd * u;
        xEuler = euler.Ad * xEuler + euler.Bd * u;
    }

    csv.close();
    std::cout << "Simulation data written to Chapter16_Lesson4_output.csv\n";
    return 0;
}

Java libraries: for production use, prefer EJML or Apache Commons Math. Below we provide a from-scratch educational implementation (dense matrices + truncated series for \( e^{\mathbf{M}T_s} \)) to make the exact-discretization logic transparent.

Code: Chapter16_Lesson4.java

// Chapter16_Lesson4.java
// Continuous–Discrete Conversions: Zero-Order Hold, Exact Discretization
// Pure Java implementation (from-scratch matrix utilities and expm series for small systems)
// For production, you may use EJML or Apache Commons Math.

import java.util.Arrays;

public class Chapter16_Lesson4 {

    static class DiscreteModel {
        double[][] Ad;
        double[][] Bd;
        DiscreteModel(double[][] Ad, double[][] Bd) { this.Ad = Ad; this.Bd = Bd; }
    }

    static double[][] zeros(int r, int c) {
        return new double[r][c];
    }

    static double[][] eye(int n) {
        double[][] I = zeros(n, n);
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] copy(double[][] A) {
        double[][] B = zeros(A.length, A[0].length);
        for (int i = 0; i < A.length; i++) System.arraycopy(A[i], 0, B[i], 0, A[0].length);
        return B;
    }

    static double[][] add(double[][] A, double[][] B) {
        int r = A.length, c = A[0].length;
        double[][] R = zeros(r, c);
        for (int i = 0; i < r; i++)
            for (int j = 0; j < c; j++)
                R[i][j] = A[i][j] + B[i][j];
        return R;
    }

    static double[][] scale(double[][] A, double s) {
        int r = A.length, c = A[0].length;
        double[][] R = zeros(r, c);
        for (int i = 0; i < r; i++)
            for (int j = 0; j < c; j++)
                R[i][j] = s * A[i][j];
        return R;
    }

    static double[][] mul(double[][] A, double[][] B) {
        int r = A.length, n = A[0].length, c = B[0].length;
        double[][] R = zeros(r, c);
        for (int i = 0; i < r; i++) {
            for (int k = 0; k < n; k++) {
                double aik = A[i][k];
                for (int j = 0; j < c; j++) {
                    R[i][j] += aik * B[k][j];
                }
            }
        }
        return R;
    }

    static double[][] blockAugment(double[][] A, double[][] B) {
        int n = A.length, m = B[0].length;
        double[][] M = zeros(n + m, n + m);
        for (int i = 0; i < n; i++)
            for (int j = 0; j < n; j++)
                M[i][j] = A[i][j];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                M[i][n + j] = B[i][j];
        return M;
    }

    // Exponential by truncated Taylor series; good for small Ts and small matrices (teaching/demo)
    static double[][] expmSeries(double[][] A, int terms) {
        int n = A.length;
        double[][] result = eye(n);
        double[][] power = eye(n);
        double factorial = 1.0;

        for (int k = 1; k <= terms; k++) {
            power = mul(power, A);
            factorial *= k;
            result = add(result, scale(power, 1.0 / factorial));
        }
        return result;
    }

    static DiscreteModel exactDiscretizeZOH(double[][] A, double[][] B, double Ts) {
        int n = A.length;
        int m = B[0].length;

        double[][] M = blockAugment(A, B);
        double[][] Md = expmSeries(scale(M, Ts), 40);

        double[][] Ad = zeros(n, n);
        double[][] Bd = zeros(n, m);

        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) Ad[i][j] = Md[i][j];
            for (int j = 0; j < m; j++) Bd[i][j] = Md[i][n + j];
        }
        return new DiscreteModel(Ad, Bd);
    }

    static DiscreteModel eulerDiscretize(double[][] A, double[][] B, double Ts) {
        int n = A.length;
        double[][] Ad = add(eye(n), scale(A, Ts));
        double[][] Bd = scale(B, Ts);
        return new DiscreteModel(Ad, Bd);
    }

    static double[] matVec(double[][] A, double[] x) {
        int r = A.length, c = A[0].length;
        double[] y = new double[r];
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) y[i] += A[i][j] * x[j];
        }
        return y;
    }

    static double[] vecAdd(double[] a, double[] b) {
        double[] r = new double[a.length];
        for (int i = 0; i < a.length; i++) r[i] = a[i] + b[i];
        return r;
    }

    static double[] scaleVec(double[] a, double s) {
        double[] r = new double[a.length];
        for (int i = 0; i < a.length; i++) r[i] = s * a[i];
        return r;
    }

    static void printMatrix(String name, double[][] A) {
        System.out.println(name + " =");
        for (double[] row : A) System.out.println(Arrays.toString(row));
        System.out.println();
    }

    public static void main(String[] args) {
        // Mass-spring-damper example
        double m = 1.0, c = 0.6, k = 4.0, b = 1.0;
        double[][] A = {
            {0.0, 1.0},
            {-k / m, -c / m}
        };
        double[][] B = {
            {0.0},
            {b / m}
        };
        double[][] C = {
            {1.0, 0.0}
        };
        double[][] D = {
            {0.0}
        };
        double Ts = 0.1;

        DiscreteModel exact = exactDiscretizeZOH(A, B, Ts);
        DiscreteModel euler = eulerDiscretize(A, B, Ts);

        printMatrix("Ad_exact", exact.Ad);
        printMatrix("Bd_exact", exact.Bd);
        printMatrix("Ad_euler", euler.Ad);
        printMatrix("Bd_euler", euler.Bd);

        int N = 60;
        double[] xExact = {0.0, 0.0};
        double[] xEuler = {0.0, 0.0};

        System.out.println("k,t,u,yExact,yEuler");
        for (int kStep = 0; kStep < N; kStep++) {
            double u = (kStep >= 5) ? 1.0 : 0.0;

            double yExact = matVec(C, xExact)[0] + D[0][0] * u;
            double yEuler = matVec(C, xEuler)[0] + D[0][0] * u;
            System.out.printf("%d,%.3f,%.1f,%.6f,%.6f%n", kStep, kStep * Ts, u, yExact, yEuler);

            xExact = vecAdd(matVec(exact.Ad, xExact), scaleVec(new double[]{exact.Bd[0][0], exact.Bd[1][0]}, u));
            xEuler = vecAdd(matVec(euler.Ad, xEuler), scaleVec(new double[]{euler.Bd[0][0], euler.Bd[1][0]}, u));
        }
    }
}

MATLAB/Simulink tools: expm, ss, c2d (Control System Toolbox), and standard plotting. The script also states a Simulink workflow for verifying sampled outputs.

Code: Chapter16_Lesson4.m

% Chapter16_Lesson4.m
% Continuous–Discrete Conversions: Zero-Order Hold, Exact Discretization
% MATLAB / Simulink-oriented implementation

clear; clc; close all;

% Mass-spring-damper example
m = 1.0; c = 0.6; k = 4.0; b = 1.0;
A = [0 1; -k/m -c/m];
B = [0; b/m];
C = [1 0];
D = 0;

Ts = 0.1;

% Exact ZOH discretization via augmented matrix exponential
n = size(A,1);
m_in = size(B,2);
M = [A B; zeros(m_in,n) zeros(m_in,m_in)];
Md = expm(M*Ts);
Ad_exact = Md(1:n,1:n);
Bd_exact = Md(1:n,n+1:n+m_in);

% Closed-form when A is nonsingular: Bd = A^{-1}(Ad - I)B
if rank(A) == n
    Bd_closed = A \ ((Ad_exact - eye(n))*B);
else
    Bd_closed = NaN(size(B));
end

% Toolbox method (Control System Toolbox): c2d with ZOH
sysc = ss(A,B,C,D);
sysd = c2d(sysc, Ts, 'zoh');

Ad_c2d = sysd.A;
Bd_c2d = sysd.B;

disp('Ad_exact ='); disp(Ad_exact);
disp('Bd_exact ='); disp(Bd_exact);
disp('Ad_c2d =');  disp(Ad_c2d);
disp('Bd_c2d =');  disp(Bd_c2d);

fprintf('||Ad_exact - Ad_c2d||_F = %.3e\n', norm(Ad_exact - Ad_c2d, 'fro'));
fprintf('||Bd_exact - Bd_c2d||_F = %.3e\n', norm(Bd_exact - Bd_c2d, 'fro'));

% Compare exact ZOH and forward Euler discretization
Ad_euler = eye(n) + Ts*A;
Bd_euler = Ts*B;

N = 120;
u = zeros(N,1);
u(6:end) = 1;  % unit step after k=5
t = (0:N-1)'*Ts;

x_exact = zeros(n,1);
x_euler = zeros(n,1);
y_exact = zeros(N,1);
y_euler = zeros(N,1);

for kstep = 1:N
    y_exact(kstep) = C*x_exact + D*u(kstep);
    y_euler(kstep) = C*x_euler + D*u(kstep);

    x_exact = Ad_exact*x_exact + Bd_exact*u(kstep);
    x_euler = Ad_euler*x_euler + Bd_euler*u(kstep);
end

figure;
plot(t, y_exact, 'LineWidth', 1.5); hold on;
plot(t, y_euler, '--', 'LineWidth', 1.3);
stairs(t, u, ':', 'LineWidth', 1.0);
grid on;
xlabel('Time [s]');
ylabel('Output / Input');
title('Exact ZOH vs Forward Euler Discretization');
legend('Exact ZOH','Forward Euler','u[k]', 'Location','best');

% Simulink note (manual workflow):
% 1) Build a continuous State-Space block with (A,B,C,D).
% 2) Add a Zero-Order Hold block with sample time Ts at the input.
% 3) Compare against a Discrete State-Space block using (Ad_exact,Bd_exact,C,D).
% 4) The sampled outputs should match at t = k*Ts (up to numerical tolerances).

Wolfram Mathematica tools: built-in MatrixExp, symbolic/numeric matrix operations, and plotting. The code uses the same augmented-matrix method and discrete recursion.

Code: Chapter16_Lesson4.nb

(* Chapter16_Lesson4.nb *)
(* Continuous–Discrete Conversions: Zero-Order Hold, Exact Discretization *)
(* Wolfram Mathematica / Wolfram Language implementation *)

ClearAll["Global`*"];

(* Mass-spring-damper example *)
m = 1.0; c = 0.6; k = 4.0; b = 1.0;
A = { {0., 1.}, {-k/m, -c/m} };
B = { {0.}, {b/m} };
C = { {1., 0.} };
D = { {0.} };
Ts = 0.1;

n = Length[A];
mIn = Dimensions[B][[2]];

(* Exact ZOH discretization via augmented matrix exponential *)
M = ArrayFlatten[{ {A, B}, {ConstantArray[0., {mIn, n}], ConstantArray[0., {mIn, mIn}]} }];
Md = MatrixExp[M Ts];

AdExact = Md[[1 ;; n, 1 ;; n]];
BdExact = Md[[1 ;; n, n + 1 ;; n + mIn]];

Print["AdExact = ", MatrixForm[AdExact]];
Print["BdExact = ", MatrixForm[BdExact]];

(* If A is nonsingular, closed-form Bd *)
If[Det[A] != 0,
  BdClosed = LinearSolve[A, (AdExact - IdentityMatrix[n]).B];
  Print["BdClosed = ", MatrixForm[BdClosed]];
];

(* Forward Euler (approximate) *)
AdEuler = IdentityMatrix[n] + Ts A;
BdEuler = Ts B;

(* Discrete simulation with a step input *)
Nsteps = 80;
uSeq = Join[ConstantArray[0., 5], ConstantArray[1., Nsteps - 5]];

xExact = {0., 0.};
xEuler = {0., 0.};

yExact = Table[0., {Nsteps}];
yEuler = Table[0., {Nsteps}];

Do[
  yExact[[k]] = (C . xExact)[[1]];
  yEuler[[k]] = (C . xEuler)[[1]];

  xExact = Flatten[AdExact . xExact + Flatten[BdExact] uSeq[[k]]];
  xEuler = Flatten[AdEuler . xEuler + Flatten[BdEuler] uSeq[[k]]];
, {k, 1, Nsteps}];

t = Range[0, Nsteps - 1] Ts;

ListLinePlot[
  {
    Transpose[{t, yExact}],
    Transpose[{t, yEuler}],
    Transpose[{t, uSeq}]
  },
  PlotLegends -> {"Exact ZOH", "Forward Euler", "u[k]"},
  AxesLabel -> {"t [s]", "value"},
  PlotRange -> All,
  GridLines -> Automatic,
  ImageSize -> Large
]

(* Optional (newer versions):
   sysc = StateSpaceModel[{A, B, C, D}];
   sysd = ToDiscreteTimeModel[sysc, Ts, Method -> "ZeroOrderHold"];
*)

Problem 1 (Exact ZOH for a First-Order Plant): Consider \( \dot{x}(t)=-3x(t)+2u(t) \) under ZOH with sample period \( T_s=0.2 \) s. Derive the exact discrete-time model \( x_{k+1}=A_d x_k + B_d u_k \).

Solution: Here \( A=-3 \) and \( B=2 \). Therefore

\[ A_d=e^{-3T_s}=e^{-0.6}, \qquad B_d=\frac{2}{3}\left(1-e^{-3T_s}\right)=\frac{2}{3}\left(1-e^{-0.6}\right). \]

Numerically, \( A_d\approx 0.5488 \) and \( B_d\approx 0.3008 \). Hence

\[ x_{k+1}\approx 0.5488\,x_k + 0.3008\,u_k. \]

Problem 2 (Augmented-Matrix Identity): Show that if \( \mathbf{M}=\begin{bmatrix}\mathbf{A}&\mathbf{B}\\ \mathbf{0}&\mathbf{0}\end{bmatrix} \), then the top-right block of \( e^{\mathbf{M}T_s} \) equals \( \int_0^{T_s} e^{\mathbf{A}\tau}\mathbf{B}\,d\tau \).

Solution: Expand the exponential as a power series:

\[ e^{\mathbf{M}T_s}=\mathbf{I}+\mathbf{M}T_s+\frac{(\mathbf{M}T_s)^2}{2!}+\cdots \]

Because \( \mathbf{M} \) is block upper triangular, every power \( \mathbf{M}^r \) has top-right block \( \mathbf{A}^{r-1}\mathbf{B} \) for \( r\ge 1 \). Thus the top-right block of the series is

\[ T_s\mathbf{B}+\frac{T_s^2}{2!}\mathbf{A}\mathbf{B}+\frac{T_s^3}{3!}\mathbf{A}^2\mathbf{B}+\cdots = \left(\int_0^{T_s} e^{\mathbf{A}\tau}\,d\tau\right)\mathbf{B} = \int_0^{T_s} e^{\mathbf{A}\tau}\mathbf{B}\,d\tau. \]

This proves the extraction rule for \( \mathbf{B}_d \).

Problem 3 (Continuous vs Discrete Stability Mapping): Let a continuous system have eigenvalues \( -1 \) and \( -4 \). With \( T_s=0.1 \), find the exact discrete eigenvalues and verify Schur stability.

Solution: Exact mapping gives

\[ \mu_1=e^{-1\cdot 0.1}=e^{-0.1}\approx 0.9048,\qquad \mu_2=e^{-4\cdot 0.1}=e^{-0.4}\approx 0.6703. \]

Since \( |\mu_1|<1 \) and \( |\mu_2|<1 \), the exact discrete model is Schur stable.

Problem 4 (Euler Approximation Error Term): Starting from the exact formula for \( \mathbf{A}_d \), show that Euler's \( \mathbf{A}_{d,\mathrm{Euler} }=\mathbf{I}+T_s\mathbf{A} \) has local matrix truncation error of order \( O(T_s^2) \).

Solution: By the matrix exponential series,

\[ \mathbf{A}_d=e^{\mathbf{A}T_s}=\mathbf{I}+T_s\mathbf{A}+\frac{T_s^2}{2}\mathbf{A}^2+O(T_s^3). \]

Therefore

\[ \mathbf{A}_d-\mathbf{A}_{d,\mathrm{Euler} } = \frac{T_s^2}{2}\mathbf{A}^2+O(T_s^3), \]

so the leading neglected term is proportional to \( T_s^2 \).

Problem 5 (Singular \( \mathbf{A} \) Case): Consider the double integrator \( \dot{\mathbf{x} }=\begin{bmatrix}0&1\\0&0\end{bmatrix}\mathbf{x}+\begin{bmatrix}0\\1\end{bmatrix}u \). Compute the exact ZOH matrices for sample period \( T_s \).

Solution: Since \( \mathbf{A}^2=\mathbf{0} \), the exponential truncates:

\[ \mathbf{A}_d=e^{\mathbf{A}T_s}=\mathbf{I}+T_s\mathbf{A} = \begin{bmatrix} 1 & T_s\\ 0 & 1 \end{bmatrix}. \]

Also,

\[ e^{\mathbf{A}\tau}=\mathbf{I}+\tau\mathbf{A} = \begin{bmatrix} 1 & \tau\\ 0 & 1 \end{bmatrix}, \quad \mathbf{B}= \begin{bmatrix} 0\\ 1 \end{bmatrix}. \]

Hence

\[ \mathbf{B}_d=\int_0^{T_s} \begin{bmatrix} \tau\\ 1 \end{bmatrix}d\tau = \begin{bmatrix} T_s^2/2\\ T_s \end{bmatrix}. \]

This example shows why the augmented-matrix method is convenient: it works even when \( \mathbf{A}^{-1} \) does not exist.

We derived the exact zero-order-hold discretization of continuous-time LTI systems, proved the formulas for \( \mathbf{A}_d \) and \( \mathbf{B}_d \), and showed a robust block-matrix exponential method that works for singular and nonsingular systems alike. We also connected exact discretization to z-domain transfer functions and to stability mapping via \( \mu=e^{\lambda T_s} \). Finally, we implemented the conversion pipeline in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.