In classical (transfer-function) analysis, the homogeneous response is characterized by poles and modes. In state space, the analogous object is the \( n\times n \) operator that maps an initial state \( \mathbf{x}(t_0) \) to \( \mathbf{x}(t) \) for the dynamics \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \). That operator is the matrix exponential \( e^{\mathbf{A}(t-t_0)} \).

Our target statement (proved rigorously below) is: for constant \( \mathbf{A}\in\mathbb{R}^{n\times n} \) and any initial condition \( \mathbf{x}(t_0)=\mathbf{x}_0 \), the unique solution is

\[ \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0. \]

flowchart TD
  A["Given: xdot = A x, initial state x(t0)=x0"] --> B["Define exp(At) by power series"]
  B --> C["Prove: d/dt exp(A(t-t0)) = A exp(A(t-t0))"]
  C --> D["Verify: x(t)=exp(A(t-t0)) x0 satisfies ODE + initial condition"]
  D --> E["Uniqueness: any other solution must coincide"]
  E --> F["Compute exp(At): eigen/Jordan/Cayley-Hamilton or numerical expm"]
        

The remainder of this lesson focuses exclusively on the homogeneous equation. Inputs and forced responses are deferred to Chapter 7, Lesson 2, where convolution-integral structure is introduced.

Consider the homogeneous LTI state equation: \( \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t) \), where \( \mathbf{x}(t)\in\mathbb{R}^n \) and \( \mathbf{A}\in\mathbb{R}^{n\times n} \) is constant. For a fixed \( t_0 \), define the map \( \mathcal{T}(t,t_0):\mathbb{R}^n \rarr \mathbb{R}^n \) by \( \mathbf{x}(t)=\mathcal{T}(t,t_0)\mathbf{x}(t_0) \).

Linearity of the differential equation implies linearity of the evolution map: for scalars \( \alpha,\beta\in\mathbb{R} \) and vectors \( \mathbf{x}_{0,1},\mathbf{x}_{0,2} \),

\[ \mathcal{T}(t,t_0)(\alpha \mathbf{x}_{0,1}+\beta \mathbf{x}_{0,2}) =\alpha \mathcal{T}(t,t_0)\mathbf{x}_{0,1}+\beta \mathcal{T}(t,t_0)\mathbf{x}_{0,2}. \]

Therefore, \( \mathcal{T}(t,t_0) \) must be representable as a matrix: there exists \( \mathbf{M}(t,t_0)\in\mathbb{R}^{n\times n} \) such that \( \mathbf{x}(t)=\mathbf{M}(t,t_0)\mathbf{x}_0 \). The matrix exponential will be shown to be exactly this matrix: \( \mathbf{M}(t,t_0)=e^{\mathbf{A}(t-t_0)} \).

From Chapter 3 (matrix exponential basics), we recall the power-series definition. For any square matrix \( \mathbf{A} \), define \( e^{\mathbf{A}t} \) by the absolutely convergent series

\[ e^{\mathbf{A}t} \triangleq \sum_{k=0}^{\infty}\frac{(\mathbf{A}t)^k}{k!} = \mathbf{I} + \mathbf{A}t + \frac{\mathbf{A}^2 t^2}{2!} + \cdots . \]

Convergence holds for all real \( t \) and all finite matrices \( \mathbf{A} \), because the scalar exponential series is absolutely convergent and matrix norms bound each term. For any submultiplicative norm \( \|\cdot\| \),

\[ \left\|\sum_{k=0}^{\infty}\frac{(\mathbf{A}t)^k}{k!}\right\| \le \sum_{k=0}^{\infty}\frac{\|\mathbf{A}\|^k |t|^k}{k!} = e^{\|\mathbf{A}\||t|} < \infty. \]

Two immediate properties (proved below using only the series definition) are:

• Identity at zero: \( e^{\mathbf{A}0}=\mathbf{I} \).
• Similarity invariance: if \( \mathbf{A}=\mathbf{T}\mathbf{B}\mathbf{T}^{-1} \), then \( e^{\mathbf{A}t}=\mathbf{T}e^{\mathbf{B}t}\mathbf{T}^{-1} \).

Theorem 1 (Homogeneous LTI solution). Let \( \mathbf{A}\in\mathbb{R}^{n\times n} \) be constant and \( \mathbf{x}(t) \) satisfy \( \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t) \) with \( \mathbf{x}(t_0)=\mathbf{x}_0 \). Then the unique solution is

\[ \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0. \]

Proof (existence by verification + uniqueness by a zero-solution argument).

Step 1: Differentiate \( e^{\mathbf{A}(t-t_0)} \). Using the series definition and term-by-term differentiation (justified by uniform convergence on bounded \( t \)-intervals), we have:

\[ \frac{d}{dt}e^{\mathbf{A}(t-t_0)} = \frac{d}{dt}\sum_{k=0}^{\infty}\frac{\mathbf{A}^k (t-t_0)^k}{k!} = \sum_{k=1}^{\infty}\frac{\mathbf{A}^k k (t-t_0)^{k-1}}{k!} = \sum_{k=1}^{\infty}\frac{\mathbf{A}^k (t-t_0)^{k-1}}{(k-1)!}. \]

Re-index with \( j=k-1 \), giving

\[ \frac{d}{dt}e^{\mathbf{A}(t-t_0)} = \sum_{j=0}^{\infty}\frac{\mathbf{A}^{j+1}(t-t_0)^j}{j!} = \mathbf{A}\sum_{j=0}^{\infty}\frac{\mathbf{A}^{j}(t-t_0)^j}{j!} = \mathbf{A}e^{\mathbf{A}(t-t_0)}. \]

Thus, \( \frac{d}{dt}e^{\mathbf{A}(t-t_0)}=\mathbf{A}e^{\mathbf{A}(t-t_0)} \).

Step 2: Verify the proposed state \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \). Differentiate:

\[ \dot{\mathbf{x}}(t)=\frac{d}{dt}\left(e^{\mathbf{A}(t-t_0)}\mathbf{x}_0\right) = \left(\frac{d}{dt}e^{\mathbf{A}(t-t_0)}\right)\mathbf{x}_0 = \mathbf{A}e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 = \mathbf{A}\mathbf{x}(t). \]

Also, at \( t=t_0 \):

\[ \mathbf{x}(t_0)=e^{\mathbf{A}(t_0-t_0)}\mathbf{x}_0=e^{\mathbf{A}0}\mathbf{x}_0=\mathbf{I}\mathbf{x}_0=\mathbf{x}_0. \]

Therefore the proposed expression satisfies the ODE and initial condition.

Step 3: Uniqueness. Suppose \( \mathbf{x}_1(t) \) and \( \mathbf{x}_2(t) \) both satisfy the same IVP. Let \( \mathbf{z}(t)=\mathbf{x}_1(t)-\mathbf{x}_2(t) \). Then \( \dot{\mathbf{z}}(t)=\mathbf{A}\mathbf{z}(t) \) and \( \mathbf{z}(t_0)=\mathbf{0} \). Consider the scalar function \( v(t)=\|\mathbf{z}(t)\| \) for any norm compatible with the induced matrix norm. Using standard estimates for linear systems (from Chapter 3 existence/uniqueness foundations),

\[ \frac{d}{dt}\|\mathbf{z}(t)\| \le \|\dot{\mathbf{z}}(t)\| = \|\mathbf{A}\mathbf{z}(t)\| \le \|\mathbf{A}\|\;\|\mathbf{z}(t)\|. \]

Grönwall’s inequality yields \( \|\mathbf{z}(t)\| \le \|\mathbf{z}(t_0)\|e^{\|\mathbf{A}\||t-t_0|}=0 \), hence \( \mathbf{z}(t)\equiv \mathbf{0} \) and \( \mathbf{x}_1(t)\equiv \mathbf{x}_2(t) \). This proves uniqueness. □

The following properties are repeatedly used throughout modern control (stability, controllability/observability, canonical forms). We prove them now using only the series definition, without introducing future chapter machinery.

Proposition 1 (Semigroup / composition property for fixed \( \mathbf{A} \)).

For all \( t,s\in\mathbb{R} \),

\[ e^{\mathbf{A}t}e^{\mathbf{A}s} = e^{\mathbf{A}(t+s)}. \]

Proof. Expand both exponentials:

\[ e^{\mathbf{A}t}e^{\mathbf{A}s} = \left(\sum_{k=0}^{\infty}\frac{\mathbf{A}^k t^k}{k!}\right)\left(\sum_{m=0}^{\infty}\frac{\mathbf{A}^m s^m}{m!}\right) = \sum_{k=0}^{\infty}\sum_{m=0}^{\infty}\frac{\mathbf{A}^{k+m} t^k s^m}{k!m!}. \]

Group terms by \( n=k+m \):

\[ \sum_{n=0}^{\infty}\mathbf{A}^{n}\left(\sum_{k=0}^{n}\frac{t^k s^{n-k}}{k!(n-k)!}\right) = \sum_{n=0}^{\infty}\mathbf{A}^n \frac{(t+s)^n}{n!} = e^{\mathbf{A}(t+s)}, \]

using the binomial identity \( \sum_{k=0}^{n}\frac{t^k s^{n-k}}{k!(n-k)!}=\frac{(t+s)^n}{n!} \). □

Corollary 1 (Inverse and nonsingularity).

\( e^{\mathbf{A}t} \) is nonsingular for all \( t \), and \( (e^{\mathbf{A}t})^{-1}=e^{-\mathbf{A}t} \).

Proof. By Proposition 1, \( e^{\mathbf{A}t}e^{-\mathbf{A}t}=e^{\mathbf{A}0}=\mathbf{I} \). □

Proposition 2 (Similarity invariance).

If \( \mathbf{A}=\mathbf{T}\mathbf{B}\mathbf{T}^{-1} \) with invertible \( \mathbf{T} \), then

\[ e^{\mathbf{A}t}=\mathbf{T}e^{\mathbf{B}t}\mathbf{T}^{-1}. \]

Proof. Use the series and the identity \( \mathbf{A}^k=(\mathbf{T}\mathbf{B}\mathbf{T}^{-1})^k=\mathbf{T}\mathbf{B}^k\mathbf{T}^{-1} \):

\[ e^{\mathbf{A}t}=\sum_{k=0}^{\infty}\frac{\mathbf{A}^k t^k}{k!} =\sum_{k=0}^{\infty}\frac{\mathbf{T}\mathbf{B}^k\mathbf{T}^{-1}t^k}{k!} =\mathbf{T}\left(\sum_{k=0}^{\infty}\frac{\mathbf{B}^k t^k}{k!}\right)\mathbf{T}^{-1} =\mathbf{T}e^{\mathbf{B}t}\mathbf{T}^{-1}. \]

□

Similarity invariance is the key reason eigenstructure and canonical forms are effective for solution analysis: compute in a convenient coordinate basis, then map back.

This section provides analytic formulas that are essential for hand analysis and proofs. Numerical algorithms are addressed afterward (still within this lesson, but treated separately).

6.1 Diagonalizable case

If \( \mathbf{A} \) is diagonalizable, then there exists an invertible \( \mathbf{V} \) and diagonal \( \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1,\dots,\lambda_n) \) such that \( \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \). By similarity invariance,

\[ e^{\mathbf{A}t}=\mathbf{V}e^{\boldsymbol{\Lambda}t}\mathbf{V}^{-1}, \quad e^{\boldsymbol{\Lambda}t}=\operatorname{diag}(e^{\lambda_1 t},\dots,e^{\lambda_n t}). \]

Consequently, the state decomposes into modal components (in coordinates \( \mathbf{z}=\mathbf{V}^{-1}\mathbf{x} \)) with independent scalar dynamics \( \dot{z}_i=\lambda_i z_i \), yielding \( z_i(t)=e^{\lambda_i(t-t_0)}z_i(t_0) \).

6.2 Jordan form case (non-diagonalizable)

If \( \mathbf{A} \) is not diagonalizable, then there exists an invertible \( \mathbf{T} \) such that \( \mathbf{A}=\mathbf{T}\mathbf{J}\mathbf{T}^{-1} \), where \( \mathbf{J} \) is Jordan. Again, similarity invariance gives \( e^{\mathbf{A}t}=\mathbf{T}e^{\mathbf{J}t}\mathbf{T}^{-1} \).

For a single Jordan block \( \mathbf{J}=\lambda \mathbf{I}+\mathbf{N} \), where \( \mathbf{N} \) is nilpotent (strictly upper triangular) and \( \mathbf{N}^p=\mathbf{0} \) for some \( p \), observe that \( \lambda \mathbf{I} \) commutes with \( \mathbf{N} \), so

\[ e^{(\lambda \mathbf{I}+\mathbf{N})t} = e^{\lambda t} e^{\mathbf{N}t} = e^{\lambda t}\sum_{k=0}^{p-1}\frac{\mathbf{N}^k t^k}{k!}. \]

This explicitly shows why repeated eigenvalues can generate polynomial-in-time factors multiplying exponentials (terms like \( t e^{\lambda t} \), \( t^2 e^{\lambda t} \), etc.).

6.3 Cayley–Hamilton reduction: \( e^{\mathbf{A}t} \) is a polynomial in \( \mathbf{A} \)

Let the characteristic polynomial of \( \mathbf{A} \) be \( p(\lambda)=\lambda^n + a_{n-1}\lambda^{n-1}+\cdots + a_0 \). By Cayley–Hamilton, \( p(\mathbf{A})=\mathbf{0} \), i.e.

\[ \mathbf{A}^n = -a_{n-1}\mathbf{A}^{n-1}-\cdots-a_1\mathbf{A}-a_0\mathbf{I}. \]

Hence any power \( \mathbf{A}^k \) for \( k\ge n \) can be reduced to a linear combination of \( \mathbf{I},\mathbf{A},\dots,\mathbf{A}^{n-1} \). Therefore, for each fixed \( t \) there exist scalar functions \( \alpha_0(t),\dots,\alpha_{n-1}(t) \) such that

\[ e^{\mathbf{A}t}=\sum_{k=0}^{n-1}\alpha_k(t)\mathbf{A}^k. \]

This statement is not merely theoretical; it underpins explicit algorithms (e.g., Putzer-type constructions) for hand computation in small dimensions without full eigenvector machinery.

In practical control computations, \( e^{\mathbf{A}t} \) is evaluated numerically. A critical caution: directly truncating the series \( \sum_{k=0}^{K}\frac{(\mathbf{A}t)^k}{k!} \) can be inaccurate for large \( \|\mathbf{A}t\| \). High-quality libraries typically implement scaling-and-squaring combined with rational (Padé-type) approximations.

A stable conceptual workflow for computation is:

flowchart TD
  S["Inputs: A, t, accuracy target"] --> N1["Compute norm estimate of A*t"]
  N1 -->|"small"| M1["Use truncated series or low-order rational approx"]
  N1 -->|"moderate/large"| M2["Scale: A*t = (A*t)/2^s"]
  M2 --> M3["Approx exp((A*t)/2^s) with Pade rational"]
  M3 --> M4["Square result s times"]
  M1 --> OUT["Return exp(A*t)"]
  M4 --> OUT
  OUT --> X["State: x(t)=exp(A*(t-t0))*x0"]

For this course, you should treat library functions (MATLAB \( \texttt{expm} \), SciPy \( \texttt{expm} \), etc.) as the reference implementation unless you are explicitly tasked with implementing the method from scratch. Nevertheless, implementing a basic method is pedagogically useful (and included in the coding sections).

Recommended libraries for modern control workflows in Python:

• NumPy for array/matrix operations.
• SciPy for numerically stable \( e^{\mathbf{A}t} \) via scipy.linalg.expm.
• SymPy for symbolic verification in small dimensions.
• python-control (the control package) for state-space objects (used more heavily in later chapters).

8.1 Compute \( x(t) \) using SciPy expm

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

A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])
t0 = 0.0
x0 = np.array([1.0, 0.0])

def x_of_t(t: float) -> np.ndarray:
    return expm(A * (t - t0)) @ x0

ts = np.linspace(0.0, 5.0, 6)
X = np.vstack([x_of_t(t) for t in ts])
print("Times:", ts)
print("States:\n", X)
      

8.2 From-scratch truncated series (educational)

This implementation is mathematically transparent but can be numerically fragile for large \( \| \mathbf{A}t \| \).

import numpy as np

def expm_series(A: np.ndarray, t: float, K: int = 30) -> np.ndarray:
    n = A.shape[0]
    I = np.eye(n)
    At = A * t
    term = I.copy()
    S = I.copy()
    for k in range(1, K + 1):
        term = term @ (At / k)   # term = (At^k)/k!
        S = S + term
    return S

A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])
t = 1.0
E_lib = expm(A * t)
E_ser = expm_series(A, t, K=40)
print("Relative error:", np.linalg.norm(E_lib - E_ser) / np.linalg.norm(E_lib))
      

8.3 Eigen-decomposition method (works when diagonalizable and well-conditioned)

import numpy as np

def expm_via_eig(A: np.ndarray, t: float) -> np.ndarray:
    # Works best if A is diagonalizable and eigenvector matrix is well-conditioned.
    w, V = np.linalg.eig(A)
    Vinv = np.linalg.inv(V)
    return V @ np.diag(np.exp(w * t)) @ Vinv

A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])
t = 1.0
E_eig = expm_via_eig(A, t)
E_lib = expm(A * t)
print("Relative error:", np.linalg.norm(E_lib - E_eig) / np.linalg.norm(E_lib))
      

Recommended C++ libraries:

• Eigen for linear algebra (widely used in control/robotics). Eigen also supports matrix functions in its unsupported modules.
• Armadillo is another option, but Eigen is more common in embedded/control toolchains.

9.1 Eigen-based computation (Matrix exponential)

The following uses Eigen’s matrix exponential implementation (typically via unsupported/Eigen/MatrixFunctions). Ensure your build includes Eigen’s unsupported module.

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

int main() {
    using Eigen::Matrix2d;
    using Eigen::Vector2d;

    Matrix2d A;
    A << 0.0, 1.0,
        -2.0, -3.0;

    double t0 = 0.0;
    double t  = 1.0;

    Vector2d x0;
    x0 << 1.0, 0.0;

    Matrix2d E = (A * (t - t0)).exp(); // matrix exponential
    Vector2d x = E * x0;

    std::cout << "exp(A*(t-t0)):\n" << E << "\n";
    std::cout << "x(t):\n" << x << "\n";
    return 0;
}
      

9.2 From-scratch truncated series (educational)

#include <iostream>
#include <Eigen/Dense>

Eigen::MatrixXd expm_series(const Eigen::MatrixXd& A, double t, int K=30) {
    const int n = A.rows();
    Eigen::MatrixXd I = Eigen::MatrixXd::Identity(n, n);
    Eigen::MatrixXd At = A * t;

    Eigen::MatrixXd term = I;
    Eigen::MatrixXd S = I;

    for (int k = 1; k <= K; ++k) {
        term = term * (At / double(k)); // term = At^k / k!
        S += term;
    }
    return S;
}

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

    double t = 1.0;
    Eigen::Matrix2d E = expm_series(A, t, 40);
    std::cout << "Series exp(A t):\n" << E << "\n";
    return 0;
}
      

Recommended Java libraries:

• EJML (Efficient Java Matrix Library) for fast dense linear algebra.
• Apache Commons Math for decompositions (eigen, etc.).

Not all Java math libraries ship a production-grade matrix exponential. For a course implementation, a carefully written truncated series (optionally with scaling) is acceptable for moderate norms. Below is a simple, readable implementation using EJML.

10.1 EJML truncated-series expm

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.dense.row.NormOps_DDRM;

public class MatrixExponentialLesson {

    // Educational series: exp(A t) ~= sum_{k=0}^K (A t)^k / k!
    public static DMatrixRMaj expmSeries(DMatrixRMaj A, double t, int K) {
        int n = A.numRows;
        DMatrixRMaj I = CommonOps_DDRM.identity(n);

        DMatrixRMaj At = new DMatrixRMaj(n, n);
        CommonOps_DDRM.scale(t, A, At);

        DMatrixRMaj term = I.copy();   // (At^0)/0! = I
        DMatrixRMaj S = I.copy();

        DMatrixRMaj tmp = new DMatrixRMaj(n, n);

        for (int k = 1; k <= K; k++) {
            // term = term * (At / k)
            CommonOps_DDRM.mult(term, At, tmp);
            CommonOps_DDRM.scale(1.0 / k, tmp, term);
            CommonOps_DDRM.addEquals(S, term);
        }
        return S;
    }

    public static void main(String[] args) {
        // A = [[0, 1], [-2, -3]]
        DMatrixRMaj A = new DMatrixRMaj(new double[][]{
                {0.0, 1.0},
                {-2.0, -3.0}
        });

        double t = 1.0;
        DMatrixRMaj E = expmSeries(A, t, 50);

        System.out.println("||A t||_F = " + NormOps_DDRM.normF(A) * t);
        System.out.println("exp(A t) approx:");
        E.print();
    }
}
      

Practical note: for larger norms, incorporate scaling: compute \( s \) such that \( \| \mathbf{A}t/2^s \| \) is moderate, approximate \( e^{\mathbf{A}t/2^s} \), then square \( s \) times. This exactly mirrors the conceptual workflow in Section 7.

MATLAB provides the canonical high-quality implementation expm. For the homogeneous response, compute \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \).

11.1 MATLAB: expm-based state solution

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

ts = linspace(0,5,6);
X = zeros(2, numel(ts));

for i = 1:numel(ts)
    t = ts(i);
    X(:,i) = expm(A*(t-t0)) * x0;
end

disp('Times:'); disp(ts);
disp('States:'); disp(X);

% Cross-check by ODE solver (same equation xdot = A x)
f = @(t,x) A*x;
[t_ode, x_ode] = ode45(f, [0 5], x0);

% Optional: compare final state
x_expm_final = expm(A*(5-t0)) * x0;
x_ode_final = x_ode(end,:)';
disp('Final state (expm):'); disp(x_expm_final);
disp('Final state (ode45):'); disp(x_ode_final);
      

11.2 Simulink: homogeneous response using a State-Space block

In Simulink, you can simulate \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) by setting the input to zero. Use a State-Space block with matrices: \( \mathbf{A} \) as given, \( \mathbf{B}=\mathbf{0} \), choose \( \mathbf{C}=\mathbf{I} \), \( \mathbf{D}=\mathbf{0} \), and set the initial condition to \( \mathbf{x}_0 \).

The script below programmatically builds a minimal model (useful for reproducible labs).

% Programmatic Simulink build for xdot = A x (u(t)=0)
A = [0 1; -2 -3];
B = [0; 0];
C = eye(2);
D = [0; 0];
x0 = [1; 0];

mdl = 'homogeneous_lti_expm_demo';
new_system(mdl); open_system(mdl);

add_block('simulink/Sources/Constant', [mdl '/u0'], 'Value', '0');
add_block('simulink/Continuous/State-Space', [mdl '/SS'], ...
    'A', 'A', 'B', 'B', 'C', 'C', 'D', 'D', 'X0', 'x0');
add_block('simulink/Sinks/To Workspace', [mdl '/x_out'], 'VariableName', 'x_sim');

set_param([mdl '/u0'], 'Position', [30 40 60 70]);
set_param([mdl '/SS'], 'Position', [120 30 260 90]);
set_param([mdl '/x_out'], 'Position', [320 40 390 70]);

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

set_param(mdl, 'StopTime', '5');
save_system(mdl);

sim(mdl);

% x_sim will contain the state trajectory (as configured by To Workspace)
      

Mathematica supports both symbolic and numeric matrix exponentials via MatrixExp. This is particularly effective for verifying analytic expressions for small matrices.

A = { {0, 1}, {-2, -3} };
t0 = 0;
x0 = {1, 0};

x[t_] := MatrixExp[A (t - t0)].x0

(* Evaluate at a few times *)
Table[{tt, x[tt]}, {tt, 0, 5, 1}]

(* Verify ODE: x'(t) == A x(t) *)
Simplify[D[x[t], t] == A.x[t]]

(* Symbolic exp(A t) *)
Simplify[MatrixExp[A t]]
      

The following problems mirror the style and difficulty commonly found in standard state-space/control texts. Solutions are provided in full detail to reinforce theorem-level understanding.

Problem 1 (Verification of the solution operator): Let \( \mathbf{A}\in\mathbb{R}^{n\times n} \) be constant and define \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \). Prove directly (using the power-series definition) that \( \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t) \) and \( \mathbf{x}(t_0)=\mathbf{x}_0 \).

Solution:

This is exactly Theorem 1, Steps 1–2. From the series definition, \( \frac{d}{dt}e^{\mathbf{A}(t-t_0)}=\mathbf{A}e^{\mathbf{A}(t-t_0)} \), hence \( \dot{\mathbf{x}}(t)=\mathbf{A}e^{\mathbf{A}(t-t_0)}\mathbf{x}_0=\mathbf{A}\mathbf{x}(t) \). Also \( \mathbf{x}(t_0)=e^{\mathbf{A}0}\mathbf{x}_0=\mathbf{x}_0 \).

Problem 2 (Semigroup property): Prove that for fixed \( \mathbf{A} \) and all \( s,t\in\mathbb{R} \), \( e^{\mathbf{A}t}e^{\mathbf{A}s}=e^{\mathbf{A}(t+s)} \), and deduce \( (e^{\mathbf{A}t})^{-1}=e^{-\mathbf{A}t} \).

Solution:

Proposition 1 proves the identity by multiplying the series and applying the binomial identity termwise. Setting \( s=-t \) yields \( e^{\mathbf{A}t}e^{-\mathbf{A}t}=e^{\mathbf{A}0}=\mathbf{I} \), so the inverse exists and equals \( e^{-\mathbf{A}t} \).

Problem 3 (Diagonal matrix case): Let \( \mathbf{A}=\operatorname{diag}(a_1,\dots,a_n) \). Compute \( e^{\mathbf{A}t} \) explicitly and show that each state component evolves independently.

Solution:

Since \( \mathbf{A}^k=\operatorname{diag}(a_1^k,\dots,a_n^k) \), the series gives

\[ e^{\mathbf{A}t}=\sum_{k=0}^{\infty}\frac{\mathbf{A}^k t^k}{k!} =\operatorname{diag}\left(\sum_{k=0}^{\infty}\frac{(a_1 t)^k}{k!},\dots,\sum_{k=0}^{\infty}\frac{(a_n t)^k}{k!}\right) =\operatorname{diag}(e^{a_1 t},\dots,e^{a_n t}). \]

Therefore \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \) implies \( x_i(t)=e^{a_i(t-t_0)}x_i(t_0) \) independently for each \( i \).

Problem 4 (Jordan block computation): Let \( \mathbf{J}=\lambda \mathbf{I}+\mathbf{N} \) where \( \mathbf{N}=\begin{bmatrix}0 & 1 \\ 0 & 0\end{bmatrix} \). Compute \( e^{\mathbf{J}t} \) explicitly.

Solution:

Here \( \mathbf{N}^2=\mathbf{0} \), so \( e^{\mathbf{N}t}=\mathbf{I}+\mathbf{N}t \). Also \( \lambda \mathbf{I} \) commutes with \( \mathbf{N} \), hence

\[ e^{\mathbf{J}t}=e^{(\lambda \mathbf{I}+\mathbf{N})t}=e^{\lambda t}e^{\mathbf{N}t} =e^{\lambda t}(\mathbf{I}+\mathbf{N}t) = e^{\lambda t}\begin{bmatrix}1 & t \\ 0 & 1\end{bmatrix}. \]

Problem 5 (Oscillatory system): Let \( \mathbf{A}=\begin{bmatrix}0 & 1 \\ -1 & 0\end{bmatrix} \). Compute \( e^{\mathbf{A}t} \) and show that it is a rotation matrix.

Solution:

Note that \( \mathbf{A}^2=-\mathbf{I} \), \( \mathbf{A}^3=-\mathbf{A} \), etc. Separate the exponential series into even and odd terms:

\[ e^{\mathbf{A}t}=\sum_{k=0}^{\infty}\frac{\mathbf{A}^k t^k}{k!} = \sum_{m=0}^{\infty}\frac{\mathbf{A}^{2m} t^{2m}}{(2m)!}+\sum_{m=0}^{\infty}\frac{\mathbf{A}^{2m+1} t^{2m+1}}{(2m+1)!}. \]

Since \( \mathbf{A}^{2m}=(-1)^m\mathbf{I} \) and \( \mathbf{A}^{2m+1}=(-1)^m\mathbf{A} \),

\[ e^{\mathbf{A}t}=\left(\sum_{m=0}^{\infty}\frac{(-1)^m t^{2m}}{(2m)!}\right)\mathbf{I} +\left(\sum_{m=0}^{\infty}\frac{(-1)^m t^{2m+1}}{(2m+1)!}\right)\mathbf{A} = (\cos t)\mathbf{I}+(\sin t)\mathbf{A}. \]

Substitute \( \mathbf{A} \) to obtain

\[ e^{\mathbf{A}t}= \begin{bmatrix} \cos t & \sin t \\ -\sin t & \cos t \end{bmatrix}, \]

which is exactly a planar rotation matrix (orthogonal with determinant 1).

We proved that the unique solution of the homogeneous LTI state equation \( \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} \) with \( \mathbf{x}(t_0)=\mathbf{x}_0 \) is \( \mathbf{x}(t)=e^{\mathbf{A}(t-t_0)}\mathbf{x}_0 \), derived core properties of the matrix exponential directly from its power series, and developed analytic computation routes via diagonalization, Jordan structure, and Cayley–Hamilton reduction. Multi-language implementations demonstrated both production-grade library usage and from-scratch educational approximations.

In the next lesson (Chapter 7, Lesson 2), we extend these ideas to the forced response with inputs, where the homogeneous operator becomes the kernel of the convolution integral.