In previous lessons, the state transition matrix was defined by \( \Phi(t)=e^{At} \), and several symbolic routes were developed: eigen-decomposition for diagonalizable matrices and Jordan form for repeated or defective eigenvalues. For a continuous-time LTI state equation

\[ \dot{\mathbf{x} }(t)=A\mathbf{x}(t),\qquad \mathbf{x}(0)=\mathbf{x}_0 , \]

the exact homogeneous solution is \( \mathbf{x}(t)=\Phi(t)\mathbf{x}_0 \). Numerically, however, the phrase “compute \( e^{At} \)” can mean several different tasks: computing the whole dense matrix, computing only \( e^{At}\mathbf{v} \), discretizing a continuous model, checking stability-sensitive simulations, or verifying a controller design.

\[ e^{At} = I + At + \frac{(At)^2}{2!}+ \frac{(At)^3}{3!}+\cdots \]

The infinite series is mathematically fundamental, but a naive truncation can be inefficient or inaccurate when \( \|At\| \) is large, when \( A \) is badly scaled, or when cancellation occurs. Practical algorithms therefore combine approximation theory, matrix factorizations, scaling, and error checks.

flowchart TD
  A["Problem: need Phi(t) or Phi(t) times vector"] --> B["Small dense A?"]
  B -->|yes| C["Use scaling and squaring \nwith Pade or Schur method"]
  B -->|no| D["Large sparse A?"]
  D -->|yes| E["Compute action exp(A t) v \nusing Krylov method"]
  D -->|no| F["Check structure: diagonal, \ntriangular, block, nilpotent"]
  C --> G["Verify residual, semigroup, inverse, conditioning"]
  E --> G
  F --> G
  G --> H["Use Phi(t) in simulation or exact discretization"]
        

The state transition matrix for an LTI system satisfies four identities that are useful as numerical diagnostics:

\[ \Phi(0)=I,\qquad \frac{d}{dt}\Phi(t)=A\Phi(t),\qquad \Phi(t+s)=\Phi(t)\Phi(s),\qquad \Phi(t)^{-1}=\Phi(-t). \]

Since \( A \) commutes with every polynomial in \( A \), it also commutes with \( e^{At} \). Hence \( A\Phi(t)=\Phi(t)A \). This implies that \( \Phi(t) \) solves both left and right matrix differential equations:

\[ \dot{\Phi}(t)=A\Phi(t)=\Phi(t)A. \]

A computed approximation \( \widehat{\Phi}(t) \) should be checked through residuals such as

\[ r_{\mathrm{diff} }(t)= \left\|\frac{\widehat{\Phi}(t+\varepsilon)- \widehat{\Phi}(t-\varepsilon)}{2\varepsilon} -A\widehat{\Phi}(t)\right\|, \]

\[ r_{\mathrm{semi} }(t,s)= \left\|\widehat{\Phi}(t)\widehat{\Phi}(s)- \widehat{\Phi}(t+s)\right\|, \qquad r_{\mathrm{inv} }(t)= \left\|\widehat{\Phi}(t)\widehat{\Phi}(-t)-I\right\|. \]

These tests do not prove correctness, but they detect many failures caused by overflow, insufficient approximation order, or inconsistent time scaling.

If \( A=V\Lambda V^{-1} \), then \( e^{At}=Ve^{\Lambda t}V^{-1} \). This is elegant, but numerically risky when \( V \) is ill conditioned. A relative perturbation in \( A \) can be magnified by approximately \( \kappa(V)=\|V\|\|V^{-1}\| \). Thus a matrix can have stable eigenvalues while its modal representation is extremely sensitive.

\[ \frac{\|\delta(e^{At})\|}{\|e^{At}\|} \;\hbox{can be large when}\; \kappa(V) \gg 1. \]

Jordan form is even more fragile. Exact Jordan block structure changes discontinuously under small perturbations: a repeated eigenvalue may split into nearby distinct eigenvalues, and a defective matrix may become diagonalizable. Hence Jordan form is valuable for theoretical interpretation, but not usually recommended as a numerical algorithm.

A more stable dense approach uses the real or complex Schur decomposition \( A=QTQ^* \), where \( Q \) is unitary/orthogonal and \( T \) is triangular or quasi-triangular. Then

\[ e^{At} = Q e^{Tt} Q^*. \]

Since \( \kappa(Q)=1 \) in the 2-norm, Schur-based algorithms avoid the instability of arbitrary eigenvector matrices.

The standard dense matrix exponential algorithm in many numerical libraries is based on scaling and squaring with a rational Padé approximation. The central identity is

\[ e^M = \left(e^{M/2^s}\right)^{2^s}, \]

where \( s \) is chosen so that \( \|M/2^s\| \) is small enough for an accurate rational approximation. For \( M=At \), choose \( s \) such that

\[ \left\|\frac{At}{2^s}\right\| \le \theta_m, \]

where \( \theta_m \) is a method-dependent threshold for a Padé approximant of order \( m \). The diagonal Padé approximant has the form

\[ r_m(Z)=p_m(Z)q_m(Z)^{-1}, \]

\[ p_m(Z)=\sum_{k=0}^m \frac{(2m-k)!m!}{(2m)!k!(m-k)!}Z^k, \qquad q_m(Z)=p_m(-Z). \]

The approximation is then squared \( s \) times:

\[ R_0=r_m\!\left(\frac{At}{2^s}\right),\qquad R_{j+1}=R_j^2,\quad j=0,1,\dots,s-1,\qquad \widehat{\Phi}(t)=R_s. \]

For teaching, truncated Taylor series after scaling is easier to code. For production, Padé is preferred because it often reaches higher accuracy with fewer matrix multiplications.

flowchart TD
  A["Input: A and time t"] --> B["Form M = A t"]
  B --> C["Estimate matrix norm"]
  C --> D["Choose scale s so norm(M / 2^s) is small"]
  D --> E["Approximate exp(M / 2^s) by Pade or Taylor"]
  E --> F["Square result s times"]
  F --> G["Return Phi(t)"]
  G --> H["Run checks: Phi(0), semigroup, inverse, residual"]
        

In digital simulation and computer control, a continuous-time system

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

is often converted to a discrete-time model under zero-order hold: \( \mathbf{u}(t)=\mathbf{u}_k \) for \( kh \le t < (k+1)h \). Then

\[ \mathbf{x}_{k+1}=A_d\mathbf{x}_k+B_d\mathbf{u}_k, \qquad A_d=e^{Ah}, \]

\[ B_d=\int_0^h e^{A\tau}B\,d\tau. \]

If \( A \) is nonsingular, one may write

\[ B_d=A^{-1}(A_d-I)B, \]

but this formula is poor when \( A \) is singular or ill conditioned. A safer method computes both \( A_d \) and \( B_d \) from one block exponential:

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

This identity follows by expanding the block exponential in a power series and observing that the upper-right block accumulates precisely the convolution integral.

For large systems, forming the full matrix \( \Phi(t) \in \mathbb{R}^{n\times n} \) may be unnecessary and expensive. If only \( \mathbf{x}(t)=e^{At}\mathbf{x}_0 \) is required, one computes the action of the exponential on a vector. Krylov methods approximate \( e^{At}\mathbf{v} \) in a subspace

\[ \mathcal{K}_m(A,\mathbf{v})= \operatorname{span}\left\{ \mathbf{v},A\mathbf{v},A^2\mathbf{v},\dots,A^{m-1}\mathbf{v} \right\}. \]

With an orthonormal basis \( V_m \) generated by the Arnoldi process,

\[ A V_m \approx V_m H_m, \]

where \( H_m \) is small and upper Hessenberg. Then

\[ e^{At}\mathbf{v}\approx \|\mathbf{v}\| V_m e^{H_m t}\mathbf{e}_1. \]

This reduces the hard computation from a large \( n\times n \) exponential to a small \( m\times m \) exponential, while using only matrix-vector products with \( A \).

Let \( u \) denote machine roundoff. A stable algorithm should produce the exact exponential of a nearby matrix:

\[ \widehat{e^A}=e^{A+\Delta A}, \qquad \frac{\|\Delta A\|}{\|A\|}=O(u). \]

This is a backward-error statement. The forward error can still be large if the exponential is ill conditioned. The Fréchet derivative of the matrix exponential at \( A \) in direction \( E \) is

\[ L_{\exp}(A,E)=\int_0^1 e^{(1-\theta)A}E e^{\theta A}\,d\theta. \]

It measures local sensitivity:

\[ e^{A+E}-e^A=L_{\exp}(A,E)+O(\|E\|^2). \]

Practical rules are:

• Avoid computing \( \Phi(t) \) by explicit Jordan form in floating point arithmetic.
• Use library functions such as scipy.linalg.expm, MATLAB expm, Wolfram MatrixExp, or Eigen's matrix-function module for dense matrices.
• For large sparse systems, prefer expm_multiply-style algorithms over forming \( e^{At} \).
• Scale state variables when entries of \( A \) have very different units or magnitudes.
• Check \( \Phi(t+s)=\Phi(t)\Phi(s) \), \( \Phi(t)^{-1}=\Phi(-t) \), and residual consistency when results look suspicious.

The Python file uses scipy.linalg.expm as the reliable dense routine, implements an educational scaling-and-squaring Taylor method, and computes exact zero-order-hold discretization by a Van Loan block exponential.

Chapter8_Lesson5.py

"""
Chapter8_Lesson5.py

Numerical computation of the state transition matrix Phi(t)=exp(A t).

This file demonstrates:
1. SciPy's robust dense matrix exponential.
2. A compact from-scratch scaling-and-squaring Taylor implementation.
3. Exact continuous-to-discrete conversion using the Van Loan block matrix.
4. Practical verification tests: semigroup, inverse, residual, and conditioning.
"""

import numpy as np
from numpy.linalg import norm, inv, cond
from scipy.linalg import expm


def expm_taylor_scaling_squaring(M, order=40):
    """
    Compute exp(M) using scaling and squaring with a truncated Taylor series.

    This is educational. For production dense problems, prefer scipy.linalg.expm,
    which uses a high-quality scaling-and-squaring Pade algorithm.
    """
    n = M.shape[0]
    m_norm = norm(M, ord=np.inf)
    s = max(0, int(np.ceil(np.log2(m_norm)))) if m_norm > 0 else 0
    A_scaled = M / (2 ** s)

    E = np.eye(n)
    term = np.eye(n)
    for k in range(1, order + 1):
        term = term @ A_scaled / k
        E = E + term

    for _ in range(s):
        E = E @ E
    return E


def exact_c2d(A, B, h):
    """
    Exact zero-order-hold discretization:
        x_{k+1} = Ad x_k + Bd u_k
    where
        Ad = exp(A h),
        Bd = integral_0^h exp(A tau) B d tau.

    The Van Loan block exponential computes both without explicitly inverting A.
    """
    n, m = B.shape
    Z = np.zeros((m, m))
    M = np.block([[A, B],
                  [np.zeros((m, n)), Z]])
    EM = expm(M * h)
    Ad = EM[:n, :n]
    Bd = EM[:n, n:n + m]
    return Ad, Bd


def residual_check(A, t, Phi, eps=1e-6):
    """
    Approximate differential residual:
        d/dt exp(A t) - A exp(A t) = 0.
    """
    Phi_plus = expm(A * (t + eps))
    Phi_minus = expm(A * (t - eps))
    dPhi_dt = (Phi_plus - Phi_minus) / (2 * eps)
    return norm(dPhi_dt - A @ Phi, ord=np.inf)


def main():
    A = np.array([[0.0, 1.0, 0.0],
                  [0.0, 0.0, 1.0],
                  [-2.0, -3.0, -4.0]])

    B = np.array([[0.0],
                  [0.0],
                  [1.0]])

    t = 0.5
    h = 0.1

    Phi_scipy = expm(A * t)
    Phi_scratch = expm_taylor_scaling_squaring(A * t, order=45)

    print("A =")
    print(A)

    print("\nPhi(t) from scipy.linalg.expm:")
    print(Phi_scipy)

    print("\nPhi(t) from scratch scaling-and-squaring Taylor:")
    print(Phi_scratch)

    print("\nInfinity-norm difference:")
    print(norm(Phi_scipy - Phi_scratch, ord=np.inf))

    print("\nSemigroup test ||Phi(t)Phi(s)-Phi(t+s)||_inf:")
    s = 0.25
    semigroup_error = norm(expm(A * t) @ expm(A * s) - expm(A * (t + s)), ord=np.inf)
    print(semigroup_error)

    print("\nInverse test ||Phi(t)Phi(-t)-I||_inf:")
    inverse_error = norm(expm(A * t) @ expm(-A * t) - np.eye(A.shape[0]), ord=np.inf)
    print(inverse_error)

    print("\nDifferential residual ||dPhi/dt - A Phi||_inf:")
    print(residual_check(A, t, Phi_scipy))

    print("\nCondition number cond(Phi(t)):")
    print(cond(Phi_scipy))

    Ad, Bd = exact_c2d(A, B, h)
    print("\nExact zero-order-hold discretization with h=0.1:")
    print("Ad =")
    print(Ad)
    print("Bd =")
    print(Bd)

    # Optional comparison for nonsingular A:
    Bd_formula = inv(A) @ (Ad - np.eye(A.shape[0])) @ B
    print("\nBd comparison with A^{-1}(Ad-I)B:")
    print(norm(Bd - Bd_formula, ord=np.inf))


if __name__ == "__main__":
    main()

The C++ implementation uses Eigen. The dense matrix exponential is available through Eigen's unsupported MatrixFunctions module. A compact Taylor scaling-and-squaring implementation is also included.

Chapter8_Lesson5.cpp

/*
Chapter8_Lesson5.cpp

Dense computation of Phi(t)=exp(A t) in C++ using Eigen.

Requires Eigen with the unsupported MatrixFunctions module:
  g++ -std=c++17 Chapter8_Lesson5.cpp -I /path/to/eigen -O2 -o Chapter8_Lesson5

The unsupported module provides matrix exponential through:
  #include <unsupported/Eigen/MatrixFunctions>
*/

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

using Eigen::MatrixXd;

MatrixXd expm_taylor_scaling_squaring(const MatrixXd& M, int order = 40) {
    const int n = static_cast<int>(M.rows());
    double inf_norm = M.cwiseAbs().rowwise().sum().maxCoeff();
    int s = 0;
    if (inf_norm > 0.0) {
        s = std::max(0, static_cast<int>(std::ceil(std::log2(inf_norm))));
    }

    MatrixXd A_scaled = M / std::pow(2.0, s);
    MatrixXd E = MatrixXd::Identity(n, n);
    MatrixXd term = MatrixXd::Identity(n, n);

    for (int k = 1; k <= order; ++k) {
        term = (term * A_scaled) / static_cast<double>(k);
        E += term;
    }

    for (int i = 0; i < s; ++i) {
        E = E * E;
    }

    return E;
}

int main() {
    MatrixXd A(3, 3);
    A << 0.0, 1.0, 0.0,
         0.0, 0.0, 1.0,
        -2.0,-3.0,-4.0;

    MatrixXd B(3, 1);
    B << 0.0, 0.0, 1.0;

    double t = 0.5;
    double h = 0.1;

    MatrixXd Phi_eigen = (A * t).exp();
    MatrixXd Phi_scratch = expm_taylor_scaling_squaring(A * t, 45);

    std::cout << "A =\n" << A << "\n\n";

    std::cout << "Phi(t) from Eigen MatrixFunctions:\n" << Phi_eigen << "\n\n";
    std::cout << "Phi(t) from scratch scaling-and-squaring Taylor:\n" << Phi_scratch << "\n\n";
    std::cout << "Infinity-norm difference:\n"
              << (Phi_eigen - Phi_scratch).cwiseAbs().rowwise().sum().maxCoeff()
              << "\n\n";

    double s = 0.25;
    MatrixXd semigroup =
        (A * t).exp() * (A * s).exp() - (A * (t + s)).exp();

    std::cout << "Semigroup error ||Phi(t)Phi(s)-Phi(t+s)||_inf:\n"
              << semigroup.cwiseAbs().rowwise().sum().maxCoeff() << "\n\n";

    MatrixXd inverse_error =
        (A * t).exp() * (-A * t).exp() - MatrixXd::Identity(3, 3);

    std::cout << "Inverse error ||Phi(t)Phi(-t)-I||_inf:\n"
              << inverse_error.cwiseAbs().rowwise().sum().maxCoeff() << "\n\n";

    // Van Loan block exponential for exact zero-order-hold discretization.
    MatrixXd M = MatrixXd::Zero(4, 4);
    M.block(0, 0, 3, 3) = A;
    M.block(0, 3, 3, 1) = B;

    MatrixXd EM = (M * h).exp();
    MatrixXd Ad = EM.block(0, 0, 3, 3);
    MatrixXd Bd = EM.block(0, 3, 3, 1);

    std::cout << "Exact ZOH discretization with h=0.1\n";
    std::cout << "Ad =\n" << Ad << "\n\n";
    std::cout << "Bd =\n" << Bd << "\n";

    return 0;
}

Java numerical libraries differ in how much direct support they provide for matrix functions. To keep the lesson self-contained, the Java file implements the exponential from scratch using arrays and scaling-and-squaring Taylor approximation.

Chapter8_Lesson5.java

/*
Chapter8_Lesson5.java

Self-contained Java implementation of Phi(t)=exp(A t) using
scaling-and-squaring with a truncated Taylor series.

For production-grade scientific computing, consider libraries such as EJML,
Apache Commons Math, or ojAlgo for matrix storage and factorizations; however,
many Java libraries do not expose a direct matrix exponential, so this file
keeps the exponential implementation explicit for teaching purposes.

Compile:
  javac Chapter8_Lesson5.java
Run:
  java Chapter8_Lesson5
*/

public class Chapter8_Lesson5 {
    static double[][] identity(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

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

    static double[][] subtract(double[][] A, double[][] B) {
        int n = A.length, m = A[0].length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                C[i][j] = A[i][j] - B[i][j];
        return C;
    }

    static double[][] scalarMultiply(double[][] A, double alpha) {
        int n = A.length, m = A[0].length;
        double[][] B = new double[n][m];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                B[i][j] = alpha * A[i][j];
        return B;
    }

    static double[][] multiply(double[][] A, double[][] B) {
        int n = A.length, p = A[0].length, m = B[0].length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++) {
            for (int k = 0; k < p; k++) {
                for (int j = 0; j < m; j++) {
                    C[i][j] += A[i][k] * B[k][j];
                }
            }
        }
        return C;
    }

    static double infNorm(double[][] A) {
        double max = 0.0;
        for (int i = 0; i < A.length; i++) {
            double rowSum = 0.0;
            for (int j = 0; j < A[0].length; j++) {
                rowSum += Math.abs(A[i][j]);
            }
            max = Math.max(max, rowSum);
        }
        return max;
    }

    static double[][] expmTaylorScalingSquaring(double[][] M, int order) {
        int n = M.length;
        double mNorm = infNorm(M);
        int scale = 0;
        if (mNorm > 0.0) {
            scale = Math.max(0, (int)Math.ceil(Math.log(mNorm) / Math.log(2.0)));
        }

        double[][] AScaled = scalarMultiply(M, 1.0 / Math.pow(2.0, scale));
        double[][] E = identity(n);
        double[][] term = identity(n);

        for (int k = 1; k <= order; k++) {
            term = scalarMultiply(multiply(term, AScaled), 1.0 / k);
            E = add(E, term);
        }

        for (int i = 0; i < scale; i++) {
            E = multiply(E, E);
        }

        return E;
    }

    static void printMatrix(String name, double[][] A) {
        System.out.println(name + " =");
        for (double[] row : A) {
            for (double x : row) {
                System.out.printf("%14.8f ", x);
            }
            System.out.println();
        }
        System.out.println();
    }

    public static void main(String[] args) {
        double[][] A = {
            { 0.0, 1.0, 0.0 },
            { 0.0, 0.0, 1.0 },
            {-2.0,-3.0,-4.0 }
        };

        double t = 0.5;
        double s = 0.25;

        double[][] At = scalarMultiply(A, t);
        double[][] As = scalarMultiply(A, s);
        double[][] A_t_plus_s = scalarMultiply(A, t + s);

        double[][] PhiT = expmTaylorScalingSquaring(At, 45);
        double[][] PhiS = expmTaylorScalingSquaring(As, 45);
        double[][] PhiTplusS = expmTaylorScalingSquaring(A_t_plus_s, 45);

        printMatrix("A", A);
        printMatrix("Phi(t)", PhiT);

        double[][] semigroupError = subtract(multiply(PhiT, PhiS), PhiTplusS);
        System.out.println("Semigroup error ||Phi(t)Phi(s)-Phi(t+s)||_inf = "
                           + infNorm(semigroupError));

        double[][] PhiMinusT = expmTaylorScalingSquaring(scalarMultiply(A, -t), 45);
        double[][] inverseError = subtract(multiply(PhiT, PhiMinusT), identity(3));
        System.out.println("Inverse error ||Phi(t)Phi(-t)-I||_inf = "
                           + infNorm(inverseError));
    }
}

MATLAB provides expm for matrix exponentials and, with the Control System Toolbox, c2d for continuous-to-discrete conversion. The file below also performs the Van Loan block exponential explicitly, which is the mathematical object behind exact ZOH discretization.

Chapter8_Lesson5.m

% Chapter8_Lesson5.m
%
% Numerical computation of the state transition matrix Phi(t)=expm(A*t).
%
% Demonstrates:
% 1. MATLAB expm for dense matrix exponentials.
% 2. Educational scaling-and-squaring Taylor implementation.
% 3. Exact continuous-to-discrete conversion using a Van Loan block matrix.
% 4. Semigroup and inverse verification tests.

clear; clc;

A = [ 0  1  0;
      0  0  1;
     -2 -3 -4 ];

B = [0; 0; 1];

t = 0.5;
s = 0.25;
h = 0.1;

Phi_builtin = expm(A*t);
Phi_scratch = local_expm_taylor_scaling_squaring(A*t, 45);

disp('A =');
disp(A);

disp('Phi(t) from MATLAB expm =');
disp(Phi_builtin);

disp('Phi(t) from scratch scaling-and-squaring Taylor =');
disp(Phi_scratch);

disp('Infinity-norm difference =');
disp(norm(Phi_builtin - Phi_scratch, inf));

semigroup_error = norm(expm(A*t)*expm(A*s) - expm(A*(t+s)), inf);
disp('Semigroup error ||Phi(t)Phi(s)-Phi(t+s)||_inf =');
disp(semigroup_error);

inverse_error = norm(expm(A*t)*expm(-A*t) - eye(size(A)), inf);
disp('Inverse error ||Phi(t)Phi(-t)-I||_inf =');
disp(inverse_error);

% Exact zero-order-hold discretization:
% x_{k+1} = Ad x_k + Bd u_k
n = size(A,1);
m = size(B,2);
M = [A B; zeros(m,n) zeros(m,m)];
EM = expm(M*h);
Ad = EM(1:n, 1:n);
Bd = EM(1:n, n+1:n+m);

disp('Exact ZOH discretization with h=0.1:');
disp('Ad =');
disp(Ad);
disp('Bd =');
disp(Bd);

% If the Control System Toolbox is available, c2d gives the same result.
if exist('ss', 'file') == 2 && exist('c2d', 'file') == 2
    C = eye(n);
    D = zeros(n,m);
    sysc = ss(A,B,C,D);
    sysd = c2d(sysc, h, 'zoh');
    disp('Control System Toolbox comparison: norm(sysd.A - Ad, inf) =');
    disp(norm(sysd.A - Ad, inf));
end

function E = local_expm_taylor_scaling_squaring(M, order)
    n = size(M,1);
    mnorm = norm(M, inf);
    if mnorm > 0
        scale = max(0, ceil(log2(mnorm)));
    else
        scale = 0;
    end

    A_scaled = M / 2^scale;
    E = eye(n);
    term = eye(n);

    for k = 1:order
        term = term * A_scaled / k;
        E = E + term;
    end

    for i = 1:scale
        E = E * E;
    end
end

In Simulink, the same continuous-time model can be represented with a State-Space block using matrices \( A,B,C,D \). For exact sampled-data simulation, use MATLAB's c2d(sys,h,'zoh') to generate the discrete matrices and place them in a Discrete State-Space block.

Mathematica exposes symbolic and numerical matrix exponentials through MatrixExp. The notebook file demonstrates \( \Phi(t) \), the semigroup and inverse checks, Van Loan discretization, and a compact educational scaling-and-squaring Taylor routine.

Chapter8_Lesson5.nb

Notebook[{
Cell["Chapter8_Lesson5.nb", "Title"],
Cell["Numerical computation of the state transition matrix Phi(t)=MatrixExp[A t].", "Text"],
Cell[BoxData[
 RowBox[{
  RowBox[{"A", "=", "{ {0,1,0},{0,0,1},{-2,-3,-4} }"}], ";",
  RowBox[{"B", "=", "{ {0},{0},{1} }"}], ";",
  RowBox[{"t", "=", "0.5"}], ";",
  RowBox[{"h", "=", "0.1"}], ";"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"Phi", "=", "N[MatrixExp[A t]]"}], ";",
  RowBox[{"MatrixForm[Phi]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"SemigroupError", "=", 
   "Norm[N[MatrixExp[A t].MatrixExp[A 0.25]-MatrixExp[A (t+0.25)]], Infinity]"}],
  ";", "SemigroupError"}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"InverseError", "=", 
   "Norm[N[MatrixExp[A t].MatrixExp[-A t]-IdentityMatrix[3]], Infinity]"}],
  ";", "InverseError"}]], "Input"],
Cell["Van Loan block exponential for exact zero-order-hold discretization.", "Text"],
Cell[BoxData[
 RowBox[{
  RowBox[{"M", "=", "ArrayFlatten[{ {A,B},{ConstantArray[0,{1,3}],ConstantArray[0,{1,1}]} }]"}], ";",
  RowBox[{"EM", "=", "N[MatrixExp[M h]]"}], ";",
  RowBox[{"Ad", "=", "EM[[1;;3,1;;3]]"}], ";",
  RowBox[{"Bd", "=", "EM[[1;;3,4;;4]]"}], ";",
  RowBox[{"{MatrixForm[Ad],MatrixForm[Bd]}"}]}]], "Input"],
Cell["Educational scaling-and-squaring Taylor implementation.", "Text"],
Cell[BoxData[
 RowBox[{
  RowBox[{"ExpmTaylorScalingSquaring[M_, order_:40] := Module[{n, mnorm, scale, As, E, term, k},",
  "n=Length[M];",
  "mnorm=Norm[M, Infinity];",
  "scale=If[mnorm>0, Max[0, Ceiling[Log[2,mnorm]]], 0];",
  "As=N[M/2^scale];",
  "E=IdentityMatrix[n];",
  "term=IdentityMatrix[n];",
  "For[k=1,k<=order,k++, term=term.As/k; E=E+term];",
  "Do[E=E.E,{scale}]; E]"}]}]], "Input"],
Cell[BoxData[
 RowBox[{
  RowBox[{"PhiScratch", "=", "ExpmTaylorScalingSquaring[A t,45]"}], ";",
  RowBox[{"Norm[PhiScratch-Phi,Infinity]"}]}]], "Input"]
}]

Consider

\[ A= \begin{bmatrix} 0&1&0\\ 0&0&1\\ -2&-3&-4 \end{bmatrix}, \qquad B= \begin{bmatrix} 0\\0\\1 \end{bmatrix}. \]

This companion-type matrix represents a third-order homogeneous equation. For \( t=0.5 \), a reliable library routine computes \( \Phi(0.5)=e^{0.5A} \). To validate the result, one should compute

\[ \left\|e^{0.5A}e^{0.25A}-e^{0.75A}\right\|, \qquad \left\|e^{0.5A}e^{-0.5A}-I\right\|. \]

If both are close to machine precision in a compatible norm, the result is consistent with the semigroup and inverse properties. Large errors indicate a scaling problem, insufficient approximation order, overflow, or an implementation bug.

Problem 1 (Series and Scaling): Let \( M=At \). Explain why computing \( e^M \) by a direct Taylor series can be inaccurate when \( \|M\| \) is large, and show how scaling and squaring reduces the difficulty.

Solution: A direct Taylor approximation \( \sum_{k=0}^N M^k/k! \) may require large \( N \) if \( \|M\| \) is large. Intermediate powers can become large before cancellation produces the final result, which increases roundoff error. Scaling uses

\[ e^M=\left(e^{M/2^s}\right)^{2^s}, \]

and chooses \( s \) so that \( \|M/2^s\| \) is moderate. The exponential of the scaled matrix is easier to approximate, and repeated squaring restores the original time scale.

Problem 2 (Semigroup Verification): Suppose a program returns approximations \( \widehat{\Phi}(t) \). Define a numerical test based on the semigroup property and explain its interpretation.

Solution: For selected values \( t,s \), compute

\[ r_{\mathrm{semi} }(t,s)= \left\|\widehat{\Phi}(t)\widehat{\Phi}(s)- \widehat{\Phi}(t+s)\right\|. \]

A small value indicates consistency with the mathematical identity \( \Phi(t+s)=\Phi(t)\Phi(s) \). A large value does not identify the exact source of error, but it signals that at least one computed transition matrix is unreliable or inconsistently scaled.

Problem 3 (Exact Discretization): Derive the exact zero-order-hold matrices \( A_d \) and \( B_d \) for \( \dot{\mathbf{x} }=A\mathbf{x}+B\mathbf{u} \) when \( \mathbf{u}(t)=\mathbf{u}_k \) on \( kh\le t < (k+1)h \).

Solution: Integrating over one sample interval gives

\[ \mathbf{x}((k+1)h) =e^{Ah}\mathbf{x}(kh) +\int_0^h e^{A\tau}B\,d\tau\;\mathbf{u}_k. \]

Therefore

\[ A_d=e^{Ah},\qquad B_d=\int_0^h e^{A\tau}B\,d\tau. \]

The Van Loan block exponential computes both:

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

Problem 4 (Ill-Conditioned Eigenvectors): Let \( A=V\Lambda V^{-1} \). Why can \( Ve^{\Lambda t}V^{-1} \) be a poor numerical way to compute \( e^{At} \)?

Solution: If \( V \) is ill conditioned, small floating-point errors in the eigenvectors, eigenvalues, or matrix multiplications are amplified by \( \kappa(V)=\|V\|\|V^{-1}\| \). The formula is mathematically correct, but it can have large forward error. Schur-based methods are preferred because the Schur vectors are orthogonal/unitary, so the similarity transformation does not amplify errors in the same way.

Problem 5 (Krylov Action Approximation): Explain why a Krylov method is useful when the system dimension is large and only \( e^{At}\mathbf{x}_0 \) is needed.

Solution: Forming the full \( n\times n \) matrix exponential costs memory and computation that may be unnecessary. A Krylov method builds a low dimensional subspace \( \mathcal{K}_m(A,\mathbf{x}_0) \) and approximates

\[ e^{At}\mathbf{x}_0 \approx \|\mathbf{x}_0\|V_m e^{H_m t}\mathbf{e}_1. \]

where \( H_m \) is small. The expensive exponential is computed only for \( H_m \), while the large system is accessed through matrix-vector products.

Numerical computation of the state transition matrix requires more than applying a formal formula. Direct Taylor series and Jordan form are important theoretically, but robust software typically uses scaling-and-squaring, rational approximation, Schur decompositions, or Krylov action methods. For sampled-data models, the Van Loan block exponential computes both \( A_d \) and \( B_d \) reliably. Practical work should always include residual, semigroup, inverse, and conditioning checks.