In the previous lesson, diagonalization gave a convenient expression for the state transition matrix: \( A=V\Lambda V^{-1} \) implies \( \Phi(t)=Ve^{\Lambda t}V^{-1} \). This requires a complete set of linearly independent eigenvectors. Some repeated-root systems do not have enough eigenvectors, so the diagonal modal expansion is not available.

Jordan form repairs this situation by using generalized eigenvectors. If \( A\in\mathbb{C}^{n\times n} \) has Jordan decomposition \( A=PJP^{-1} \), then

\[ \Phi(t)=e^{At}=P e^{Jt} P^{-1}. \]

The computational task is therefore reduced to computing exponentials of Jordan blocks rather than the full matrix at once.

flowchart TD
  A["Start with xdot = A x"] --> B["Check if A has enough eigenvectors"]
  B -->|yes| C["Use diagonal modal form"]
  B -->|no| D["Build Jordan chains"]
  D --> E["Form A = P J P^-1"]
  E --> F["Compute exp(J t) block by block"]
  F --> G["Phi(t) = P exp(J t) P^-1"]
  G --> H["State solution x(t) = Phi(t) x(0)"]
        

A Jordan block of size \( m \) associated with eigenvalue \( \lambda \) is

\[ J_m(\lambda)= \begin{bmatrix} \lambda & 1 & 0 & \cdots & 0\\ 0 & \lambda & 1 & \cdots & 0\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \cdots & \lambda & 1\\ 0 & 0 & \cdots & 0 & \lambda \end{bmatrix} =\lambda I_m+N_m. \]

The matrix \( N_m \) is nilpotent: after multiplying it by itself \( m \) times, the result is zero:

\[ N_m^m=0. \]

Since \( \lambda I_m \) commutes with \( N_m \), the exponential separates:

\[ e^{J_m(\lambda)t} =e^{(\lambda I_m+N_m)t} =e^{\lambda t}e^{N_m t}. \]

The nilpotent exponential is finite, not infinite:

\[ e^{N_m t}=I_m+tN_m+\frac{t^2}{2!}N_m^2+\cdots+ \frac{t^{m-1} }{(m-1)!}N_m^{m-1}. \]

flowchart TD
  JB["Jordan block J(lambda)"] --> S["Split as lambda I + N"]
  S --> C["lambda I commutes with N"]
  C --> E["exp(J t)=exp(lambda t) exp(N t)"]
  E --> T["N is nilpotent, so series terminates"]
  T --> M["Entries are exp(lambda t), t exp(lambda t), t^2/2 exp(lambda t), ..."]
        

The finite series above gives the explicit upper-triangular matrix

\[ e^{J_m(\lambda)t} = e^{\lambda t} \begin{bmatrix} 1 & t & \frac{t^2}{2!} & \cdots & \frac{t^{m-1} }{(m-1)!}\\ 0 & 1 & t & \cdots & \frac{t^{m-2} }{(m-2)!}\\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \cdots & 1 & t\\ 0 & 0 & \cdots & 0 & 1 \end{bmatrix}. \]

Equivalently, the entry in row \( i \) and column \( j \) is

\[ \left[e^{J_m(\lambda)t}\right]_{ij} = \begin{cases} e^{\lambda t}\dfrac{t^{j-i} }{(j-i)!}, & j-i\in\{0,1,\ldots,m-1\},\\ 0, & \text{otherwise}. \end{cases} \]

For example, a size-two defective block produces \( te^{\lambda t} \) terms:

\[ J_2(\lambda)= \begin{bmatrix}\lambda&1\\0&\lambda\end{bmatrix}, \qquad e^{J_2(\lambda)t}= e^{\lambda t} \begin{bmatrix}1&t\\0&1\end{bmatrix}. \]

Thus repeated defective eigenvalues create modal components that are not merely \( e^{\lambda t} \), but polynomial multiples of \( e^{\lambda t} \). In control interpretation, these terms affect transient growth even when the real part of \( \lambda \) is negative.

The matrix exponential is defined by the convergent power series

\[ e^{At}=\sum_{k=0}^{\infty}\frac{(At)^k}{k!}. \]

Suppose \( A=PJP^{-1} \). Then \( A^k=PJ^kP^{-1} \) for every nonnegative integer \( k \), by repeated cancellation of \( P^{-1}P \). Substituting into the series gives

\[ e^{At} = \sum_{k=0}^{\infty}\frac{P(Jt)^kP^{-1} }{k!} = P\left(\sum_{k=0}^{\infty}\frac{(Jt)^k}{k!}\right)P^{-1} = Pe^{Jt}P^{-1}. \]

If \( J \) is block diagonal, \( J=\operatorname{diag}(J_1,J_2,\ldots,J_q) \), then powers and exponentials preserve that block diagonal structure:

\[ e^{Jt}= \operatorname{diag}\left(e^{J_1t},e^{J_2t},\ldots,e^{J_qt}\right). \]

This proves that the full state transition matrix can be assembled from independent Jordan-block exponentials.

Consider the Jordan matrix

\[ J= \begin{bmatrix} 2&1&0\\ 0&2&0\\ 0&0&-1 \end{bmatrix} = \operatorname{diag}\left(J_2(2),J_1(-1)\right). \]

Its exponential is

\[ e^{Jt}= \begin{bmatrix} e^{2t}&te^{2t}&0\\ 0&e^{2t}&0\\ 0&0&e^{-t} \end{bmatrix}. \]

If the original state matrix is \( A=PJP^{-1} \), with

\[ P= \begin{bmatrix} 1&1&0\\ 0&1&1\\ 1&0&1 \end{bmatrix}, \]

then the transition matrix is simply

\[ \Phi(t)=P \begin{bmatrix} e^{2t}&te^{2t}&0\\ 0&e^{2t}&0\\ 0&0&e^{-t} \end{bmatrix} P^{-1}. \]

This expression is exact. The term \( te^{2t} \) comes from the generalized eigenvector associated with the repeated eigenvalue \( 2 \). If this eigenvalue were stable instead, for example \( -2 \), the term \( te^{-2t} \) would initially shape the transient but still decay asymptotically.

The symbolic Jordan-form computation of \( \Phi(t) \) proceeds as follows:

1. Find the eigenvalues of \( A \) and their algebraic multiplicities.
2. For each repeated eigenvalue, construct eigenvectors and generalized eigenvectors satisfying a chain relation such as \( (A-\lambda I)v_1=0 \) and \( (A-\lambda I)v_{k+1}=v_k \).
3. Assemble the Jordan basis \( P \) from all Jordan chains and form \( J=P^{-1}AP \).
4. Compute each block exponential using the finite nilpotent series.
5. Transform back: \( \Phi(t)=Pe^{Jt}P^{-1} \).

In numerical control software, one rarely computes Jordan form for large ill-conditioned systems. It is mainly a theoretical and symbolic tool. For numerical simulation, robust algorithms based on scaling, squaring, Padé approximation, and Schur decompositions are preferred. Nevertheless, Jordan form is essential for understanding the exact modal structure of defective systems.

The following file computes \( \Phi(t) \) from a supplied Jordan basis and validates it against scipy.linalg.expm.

Chapter8_Lesson4.py

# Chapter8_Lesson4.py
# Computing the state transition matrix Phi(t) using Jordan form.
# Requires: numpy, scipy (for validation only).

import numpy as np
from math import factorial
from scipy.linalg import expm

def jordan_block_exp(lam: complex, size: int, t: float) -> np.ndarray:
    """Return exp((lam I + N)t) for one Jordan block."""
    B = np.zeros((size, size), dtype=complex)
    for k in range(size):
        for j in range(k, size):
            power = j - k
            B[k, j] = np.exp(lam * t) * (t ** power) / factorial(power)
    return B

def block_diag(blocks):
    n = sum(B.shape[0] for B in blocks)
    M = np.zeros((n, n), dtype=complex)
    r = 0
    for B in blocks:
        m = B.shape[0]
        M[r:r+m, r:r+m] = B
        r += m
    return M

def phi_from_jordan(P: np.ndarray, blocks: list[tuple[complex, int]], t: float) -> np.ndarray:
    """Compute Phi(t)=P exp(Jt) P^{-1} from a supplied Jordan basis."""
    exp_blocks = [jordan_block_exp(lam, size, t) for lam, size in blocks]
    expJ = block_diag(exp_blocks)
    return P @ expJ @ np.linalg.inv(P)

if __name__ == "__main__":
    # Example: one size-2 Jordan block at lambda=2 and one size-1 block at lambda=-1.
    J = np.array([[2.0, 1.0, 0.0],
                  [0.0, 2.0, 0.0],
                  [0.0, 0.0, -1.0]], dtype=complex)

    # A nontrivial coordinate basis. The state matrix is A=P J P^{-1}.
    P = np.array([[1.0, 1.0, 0.0],
                  [0.0, 1.0, 1.0],
                  [1.0, 0.0, 1.0]], dtype=complex)
    A = P @ J @ np.linalg.inv(P)

    t = 0.40
    Phi_jordan = phi_from_jordan(P, [(2.0, 2), (-1.0, 1)], t)
    Phi_scipy = expm(A * t)

    np.set_printoptions(precision=6, suppress=True)
    print("A =")
    print(A)
    print("\nPhi(t) using Jordan form =")
    print(Phi_jordan)
    print("\nValidation with scipy.linalg.expm(A*t) =")
    print(Phi_scipy)
    print("\nFrobenius error =", np.linalg.norm(Phi_jordan - Phi_scipy, ord="fro"))

    x0 = np.array([1.0, -0.5, 0.25], dtype=complex)
    x_t = Phi_jordan @ x0
    print("\nx(t)=Phi(t)x(0) for x0=[1,-0.5,0.25]^T:")
    print(x_t)

The C++ version implements the necessary matrix operations from scratch: multiplication, Gauss-Jordan inverse, block construction, and Jordan block exponentiation.

Chapter8_Lesson4.cpp

// Chapter8_Lesson4.cpp
// Computing Phi(t) from a known Jordan form: Phi(t)=P exp(Jt) P^{-1}.
// Compile: g++ -std=c++17 Chapter8_Lesson4.cpp -o Chapter8_Lesson4

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

Matrix zeros(int r, int c) {
    return Matrix(r, std::vector<double>(c, 0.0));
}

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

Matrix multiply(const Matrix& A, const Matrix& B) {
    int r = static_cast<int>(A.size());
    int m = static_cast<int>(A[0].size());
    int c = static_cast<int>(B[0].size());
    Matrix C = zeros(r, c);
    for (int i = 0; i < r; ++i)
        for (int k = 0; k < m; ++k)
            for (int j = 0; j < c; ++j)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}

Matrix inverse(Matrix A) {
    int n = static_cast<int>(A.size());
    Matrix Aug = zeros(n, 2 * n);
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) Aug[i][j] = A[i][j];
        Aug[i][n + i] = 1.0;
    }

    for (int col = 0; col < n; ++col) {
        int pivot = col;
        for (int r = col + 1; r < n; ++r)
            if (std::fabs(Aug[r][col]) > std::fabs(Aug[pivot][col])) pivot = r;
        if (std::fabs(Aug[pivot][col]) < 1e-12) throw std::runtime_error("Singular matrix.");
        std::swap(Aug[col], Aug[pivot]);

        double p = Aug[col][col];
        for (int j = 0; j < 2 * n; ++j) Aug[col][j] /= p;

        for (int r = 0; r < n; ++r) {
            if (r == col) continue;
            double factor = Aug[r][col];
            for (int j = 0; j < 2 * n; ++j) Aug[r][j] -= factor * Aug[col][j];
        }
    }

    Matrix Inv = zeros(n, n);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            Inv[i][j] = Aug[i][n + j];
    return Inv;
}

double factorial(int n) {
    double f = 1.0;
    for (int i = 2; i <= n; ++i) f *= i;
    return f;
}

Matrix jordanBlockExp(double lambda, int size, double t) {
    Matrix B = zeros(size, size);
    double e = std::exp(lambda * t);
    for (int i = 0; i < size; ++i) {
        for (int j = i; j < size; ++j) {
            int p = j - i;
            B[i][j] = e * std::pow(t, p) / factorial(p);
        }
    }
    return B;
}

Matrix blockDiag(const std::vector<Matrix>& blocks) {
    int n = 0;
    for (const auto& B : blocks) n += static_cast<int>(B.size());
    Matrix M = zeros(n, n);
    int offset = 0;
    for (const auto& B : blocks) {
        int m = static_cast<int>(B.size());
        for (int i = 0; i < m; ++i)
            for (int j = 0; j < m; ++j)
                M[offset + i][offset + j] = B[i][j];
        offset += m;
    }
    return M;
}

void printMatrix(const Matrix& A, const std::string& name) {
    std::cout << name << " =\n";
    for (const auto& row : A) {
        for (double v : row) std::cout << std::setw(14) << std::setprecision(7) << v << " ";
        std::cout << "\n";
    }
}

int main() {
    Matrix P = {
        {1.0, 1.0, 0.0},
        {0.0, 1.0, 1.0},
        {1.0, 0.0, 1.0}
    };

    Matrix expJ = blockDiag({
        jordanBlockExp(2.0, 2, 0.40),
        jordanBlockExp(-1.0, 1, 0.40)
    });

    Matrix Phi = multiply(multiply(P, expJ), inverse(P));
    printMatrix(Phi, "Phi(0.40)");

    return 0;
}

The Java version mirrors the C++ implementation and is useful when students need a pure language-level implementation without depending on external numerical libraries.

Chapter8_Lesson4.java

// Chapter8_Lesson4.java
// Computing Phi(t) from a known Jordan form: Phi(t)=P exp(Jt) P^{-1}.
// Compile: javac Chapter8_Lesson4.java
// Run:     java Chapter8_Lesson4

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

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

    static double factorial(int n) {
        double f = 1.0;
        for (int i = 2; i <= n; i++) f *= i;
        return f;
    }

    static double[][] jordanBlockExp(double lambda, int size, double t) {
        double[][] B = zeros(size, size);
        double e = Math.exp(lambda * t);
        for (int i = 0; i < size; i++) {
            for (int j = i; j < size; j++) {
                int p = j - i;
                B[i][j] = e * Math.pow(t, p) / factorial(p);
            }
        }
        return B;
    }

    static double[][] blockDiag(double[][]... blocks) {
        int n = 0;
        for (double[][] B : blocks) n += B.length;
        double[][] M = zeros(n, n);
        int off = 0;
        for (double[][] B : blocks) {
            for (int i = 0; i < B.length; i++) {
                for (int j = 0; j < B.length; j++) {
                    M[off + i][off + j] = B[i][j];
                }
            }
            off += B.length;
        }
        return M;
    }

    static double[][] inverse(double[][] A) {
        int n = A.length;
        double[][] aug = zeros(n, 2 * n);

        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) aug[i][j] = A[i][j];
            aug[i][n + i] = 1.0;
        }

        for (int col = 0; col < n; col++) {
            int pivot = col;
            for (int r = col + 1; r < n; r++) {
                if (Math.abs(aug[r][col]) > Math.abs(aug[pivot][col])) pivot = r;
            }
            if (Math.abs(aug[pivot][col]) < 1e-12) {
                throw new RuntimeException("Singular matrix.");
            }

            double[] tmp = aug[col];
            aug[col] = aug[pivot];
            aug[pivot] = tmp;

            double p = aug[col][col];
            for (int j = 0; j < 2 * n; j++) aug[col][j] /= p;

            for (int r = 0; r < n; r++) {
                if (r == col) continue;
                double factor = aug[r][col];
                for (int j = 0; j < 2 * n; j++) aug[r][j] -= factor * aug[col][j];
            }
        }

        double[][] inv = zeros(n, n);
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) inv[i][j] = aug[i][n + j];
        }
        return inv;
    }

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

    public static void main(String[] args) {
        double[][] P = {
            {1.0, 1.0, 0.0},
            {0.0, 1.0, 1.0},
            {1.0, 0.0, 1.0}
        };

        double[][] expJ = blockDiag(
            jordanBlockExp(2.0, 2, 0.40),
            jordanBlockExp(-1.0, 1, 0.40)
        );

        double[][] Phi = multiply(multiply(P, expJ), inverse(P));
        printMatrix(Phi, "Phi(0.40)");
    }
}

MATLAB can compute the symbolic expression directly and validate it against expm. In Simulink, a zero-input continuous State-Space block with matrix \( A \) and initial condition \( x(0) \) has output trajectory \( x(t)=\Phi(t)x(0) \).

Chapter8_Lesson4.m

% Chapter8_Lesson4.m
% Computing Phi(t) using Jordan form in MATLAB.
% Requires Symbolic Math Toolbox for jordan(); otherwise use the supplied P,J.

clear; clc;

syms t real

J = [2 1 0;
     0 2 0;
     0 0 -1];

P = [1 1 0;
     0 1 1;
     1 0 1];

A = P * J / P;

% Exact symbolic transition matrix from the known Jordan form.
expJ = [exp(2*t), t*exp(2*t), 0;
        0,        exp(2*t),   0;
        0,        0,          exp(-t)];

Phi = simplify(P * expJ / P);

disp('A =');
disp(A);
disp('Phi(t) =');
pretty(Phi);

% Numerical validation against MATLAB expm.
tn = 0.40;
PhiJordan = double(subs(Phi, t, tn));
PhiExpm = expm(double(A) * tn);

disp('Phi(0.40) from Jordan form =');
disp(PhiJordan);
disp('Phi(0.40) from expm(A*t) =');
disp(PhiExpm);
disp('Frobenius error =');
disp(norm(PhiJordan - PhiExpm, 'fro'));

% Simulink connection:
% A Continuous State-Space block with matrix A and initial condition x0 evolves as
% x(t)=Phi(t)x0 when the input is zero. Use Phi above to verify the block output.

Mathematica is especially convenient for exact symbolic verification of the equality \( Pe^{Jt}P^{-1}=e^{At} \).

Chapter8_Lesson4.nb

(* Chapter8_Lesson4.nb *)
(* Computing Phi(t) using Jordan form in Wolfram Mathematica. *)

Clear["Global`*"];

J = { {2, 1, 0}, {0, 2, 0}, {0, 0, -1} };
P = { {1, 1, 0}, {0, 1, 1}, {1, 0, 1} };
A = Simplify[P . J . Inverse[P]];

expJ[t_] := { {Exp[2 t], t Exp[2 t], 0},
             {0, Exp[2 t], 0},
             {0, 0, Exp[-t]} };

PhiJordan[t_] := Simplify[P . expJ[t] . Inverse[P]];
PhiBuiltIn[t_] := MatrixExp[A t];

FullSimplify[PhiJordan[t] == PhiBuiltIn[t]]

N[PhiJordan[0.40] // MatrixForm]
N[Norm[PhiJordan[0.40] - PhiBuiltIn[0.40], "Frobenius"]]

x0 = {1, -1/2, 1/4};
xt = Simplify[PhiJordan[t] . x0]

Problem 1 (Exponential of a Size-Three Jordan Block): Let \( J_3(\lambda)=\lambda I_3+N_3 \). Derive \( e^{J_3(\lambda)t} \).

Solution: Since \( N_3^3=0 \),

\[ e^{J_3(\lambda)t} = e^{\lambda t} \left(I_3+tN_3+\frac{t^2}{2}N_3^2\right) = e^{\lambda t} \begin{bmatrix} 1&t&\frac{t^2}{2}\\ 0&1&t\\ 0&0&1 \end{bmatrix}. \]

Problem 2 (Transition Matrix from a Known Jordan Form): Suppose \( A=PJP^{-1} \), where \( J=\operatorname{diag}(J_2(-1),J_1(3)) \). Write \( \Phi(t) \).

Solution: The block exponential is

\[ e^{Jt}= \operatorname{diag} \left( e^{-t}\begin{bmatrix}1&t\\0&1\end{bmatrix}, e^{3t} \right). \]

Therefore, \( \Phi(t)=Pe^{Jt}P^{-1} \).

Problem 3 (Generalized Eigenvector Interpretation): For a Jordan chain \( v_1,v_2 \) satisfying \( (A-\lambda I)v_1=0 \) and \( (A-\lambda I)v_2=v_1 \), show that the state initialized at \( x(0)=v_2 \) contains a \( te^{\lambda t}v_1 \) component.

Solution: In Jordan coordinates, initializing the second vector of a size-two chain corresponds to \( z(0)=\begin{bmatrix}0&1\end{bmatrix}^T \). Since

\[ e^{J_2(\lambda)t} \begin{bmatrix}0\\1\end{bmatrix} = e^{\lambda t} \begin{bmatrix}t\\1\end{bmatrix}, \]

transforming back to physical coordinates gives \( x(t)=e^{\lambda t}(t v_1+v_2) \).

Problem 4 (Stability Insight): Does the polynomial factor \( t^k \) make a Jordan block unstable when \( \operatorname{Re}(\lambda) \) is negative?

Solution: No. Any finite polynomial grows slower than a decaying exponential. Thus \( t^k e^{\lambda t} \) tends to zero when \( \operatorname{Re}(\lambda) \) is negative. However, the polynomial factor can create large transient amplification before decay, especially when the Jordan block is large or the basis \( P \) is ill-conditioned.

Problem 5 (Why Jordan Form Is Not the Main Numerical Tool): Explain why Jordan form is useful theoretically but often avoided in numerical computation.

Solution: Jordan structure is extremely sensitive to small perturbations. A defective matrix can become diagonalizable after a tiny perturbation, and nearly defective matrices can have ill-conditioned eigenvector bases. Therefore, numerical software usually evaluates \( e^{At} \) using stable transformations such as the Schur decomposition combined with scaling-and-squaring and Padé approximation. Jordan form remains valuable for deriving exact modal terms and understanding generalized eigenvector dynamics.

Jordan form extends modal computation of the state transition matrix to matrices that lack a complete eigenvector basis. Each Jordan block \( J_m(\lambda)=\lambda I_m+N_m \) produces a finite polynomial-exponential expression because the nilpotent matrix \( N_m \) has a terminating exponential series. The complete formula \( \Phi(t)=Pe^{Jt}P^{-1} \) is exact and exposes the origin of terms such as \( te^{\lambda t} \), which are central to the dynamics of repeated defective modes.