In previous lessons, you studied scalar ordinary differential equations (ODEs) and Laplace transforms. Realistic control systems (mechanical, electrical, mechatronic) are rarely scalar; they contain several coupled energy-storage elements. The natural mathematical representation is a vector of variables, coupled through a matrix of coefficients.

A typical linear, time-invariant (LTI) system of first-order ODEs can be written as

\[ \frac{d\mathbf{x}(t)}{dt} = \mathbf{A}\,\mathbf{x}(t) + \mathbf{B}\,\mathbf{u}(t), \quad \mathbf{x}(t) \in \mathbb{R}^n,\; \mathbf{u}(t) \in \mathbb{R}^m. \]

Here \( \mathbf{x}(t) \) is a vector collecting dynamic variables (positions, velocities, currents, temperatures, etc.), and \( \mathbf{A} \) is an \( n \times n \) matrix whose entries encode coupling between them. The structure of \( \mathbf{A} \) determines whether the system is stable, oscillatory, or unstable. Understanding this matrix through eigenvalues and eigenvectors is therefore central to system dynamics.

flowchart TD
  A["Physical laws (Newton, KVL/KCL, conservation)"] --> B["Scalar ODEs"]
  B --> C["Collect variables into vector x(t)"]
  C --> D["Write coupled system: dx/dt = A x + B u"]
  D --> E["Analyze A: eigenvalues lambda"]
  E --> F["Modes: exp(lambda * t)"]
  F --> G["Stability, oscillations, transients"]
        

A vector \( \mathbf{x} \in \mathbb{R}^n \) is an ordered collection of real numbers. In dynamic systems, it is typical to interpret components as physical variables:

\[ \mathbf{x}(t) = \begin{bmatrix} x_1(t) \\ x_2(t) \\ \vdots \\ x_n(t) \end{bmatrix}. \]

To quantify the “size” of a vector, we use norms. The most common is the Euclidean (2-)norm:

\[ \|\mathbf{x}\|_2 = \sqrt{x_1^2 + x_2^2 + \dots + x_n^2} = \sqrt{\mathbf{x}^\top \mathbf{x}}. \]

The associated inner product between two vectors \( \mathbf{x}, \mathbf{y} \in \mathbb{R}^n \) is

\[ \langle \mathbf{x}, \mathbf{y} \rangle = \mathbf{x}^\top \mathbf{y} = \sum_{i=1}^n x_i y_i. \]

In dynamic systems, norms of \( \mathbf{x}(t) \) are used to measure overall “energy” or “magnitude” of the state, and the inner product describes projections of trajectories onto modes.

A matrix \( \mathbf{A} \in \mathbb{R}^{n \times n} \) represents a linear map \( \mathbf{A} : \mathbb{R}^n \to \mathbb{R}^n \) via \( \mathbf{y} = \mathbf{A}\mathbf{x} \). Linearity means

\[ \mathbf{A}(\alpha \mathbf{x} + \beta \mathbf{y}) = \alpha \mathbf{A}\mathbf{x} + \beta \mathbf{A}\mathbf{y} \quad \forall\, \alpha,\beta \in \mathbb{R},\; \mathbf{x},\mathbf{y} \in \mathbb{R}^n. \]

Consider a second-order mechanical system (mass–spring–damper) already familiar from ODEs. By introducing velocity as an additional variable, we obtain a first-order vector ODE of the form

\[ \frac{d}{dt} \begin{bmatrix} x_1(t) \\ x_2(t) \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -\dfrac{k}{m} & -\dfrac{c}{m} \end{bmatrix} \begin{bmatrix} x_1(t) \\ x_2(t) \end{bmatrix} + \begin{bmatrix} 0 \\ \dfrac{1}{m} \end{bmatrix} u(t). \]

Here \( x_1 \) is displacement, \( x_2 \) is velocity, \( u(t) \) is the external force, and the system matrix \( \mathbf{A} \) captures the intrinsic dynamics. All future chapters will work heavily with such matrix representations.

Let \( \mathbf{A} \in \mathbb{R}^{n \times n} \). A nonzero vector \( \mathbf{v} \in \mathbb{C}^n \) is called an eigenvector of \( \mathbf{A} \) with eigenvalue \( \lambda \in \mathbb{C} \) if

\[ \mathbf{A}\mathbf{v} = \lambda \mathbf{v}, \quad \mathbf{v} \neq \mathbf{0}. \]

Rearranging,

\[ (\mathbf{A} - \lambda \mathbf{I})\mathbf{v} = \mathbf{0}. \]

A nontrivial solution \( \mathbf{v} \neq \mathbf{0} \) exists if and only if the matrix \( \mathbf{A} - \lambda \mathbf{I} \) is singular, i.e. its determinant is zero:

\[ \det(\mathbf{A} - \lambda \mathbf{I}) = 0. \]

The polynomial \( p_{\mathbf{A}}(\lambda) = \det(\mathbf{A} - \lambda \mathbf{I}) \) is the characteristic polynomial of \( \mathbf{A} \). Its roots are the eigenvalues of \( \mathbf{A} \), with multiplicities known as algebraic multiplicities. The dimension of the eigenspace \( \mathcal{E}_\lambda = \ker(\mathbf{A} - \lambda \mathbf{I}) \) is the geometric multiplicity.

A central case in dynamic systems is when \( \mathbf{A} \) is diagonalizable: there exists a basis of eigenvectors forming an invertible matrix \( \mathbf{S} = [\mathbf{v}_1 \ \dots \ \mathbf{v}_n] \) such that

\[ \mathbf{A} = \mathbf{S}\boldsymbol{\Lambda}\mathbf{S}^{-1}, \quad \boldsymbol{\Lambda} = \begin{bmatrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_n \end{bmatrix}. \]

When diagonalization is possible, analysis of the dynamic system becomes much simpler, as we will see next.

Consider the homogeneous system (no input)

\[ \frac{d\mathbf{x}(t)}{dt} = \mathbf{A}\mathbf{x}(t), \quad \mathbf{x}(0) = \mathbf{x}_0. \]

We seek solutions of the form \( \mathbf{x}(t) = e^{\lambda t}\mathbf{v} \). Substituting:

\[ \frac{d}{dt}(e^{\lambda t}\mathbf{v}) = \lambda e^{\lambda t}\mathbf{v} \quad\Rightarrow\quad \lambda e^{\lambda t}\mathbf{v} = \mathbf{A} e^{\lambda t}\mathbf{v} \quad\Rightarrow\quad \mathbf{A}\mathbf{v} = \lambda \mathbf{v}. \]

Thus exponential solutions correspond exactly to eigenpairs \( (\lambda,\mathbf{v}) \) of \( \mathbf{A} \). If \( \mathbf{A} \) is diagonalizable with eigenvalues \( \lambda_i \) and eigenvectors \( \mathbf{v}_i \), any initial condition can be written as

\[ \mathbf{x}_0 = \sum_{i=1}^n \alpha_i \mathbf{v}_i. \]

The corresponding solution is then

\[ \mathbf{x}(t) = \sum_{i=1}^n \alpha_i e^{\lambda_i t} \mathbf{v}_i. \]

Equivalently, introducing \( \mathbf{S} \) and \( \boldsymbol{\Lambda} \) as before:

\[ \mathbf{x}(t) = \mathbf{S} e^{\boldsymbol{\Lambda} t} \mathbf{S}^{-1} \mathbf{x}_0, \quad e^{\boldsymbol{\Lambda} t} = \begin{bmatrix} e^{\lambda_1 t} & 0 & \cdots & 0 \\ 0 & e^{\lambda_2 t} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & e^{\lambda_n t} \end{bmatrix}. \]

The matrix \( e^{\mathbf{A}t} = \mathbf{S} e^{\boldsymbol{\Lambda} t} \mathbf{S}^{-1} \) is called the matrix exponential and acts as the propagator of the dynamics.

More generally, for any square matrix:

\[ e^{\mathbf{A}t} = \sum_{k=0}^{\infty} \frac{(\mathbf{A}t)^k}{k!}, \quad t \in \mathbb{R}. \]

Diagonalization provides a practical way to evaluate this series in closed form when possible.

From Lesson 2, you know that for a scalar ODE with constant coefficients, the Laplace transform leads to a rational transfer function whose poles are roots of the characteristic polynomial. For the vector system

\[ \frac{d\mathbf{x}(t)}{dt} = \mathbf{A}\mathbf{x}(t), \quad \mathbf{x}(0) = \mathbf{x}_0, \]

taking the Laplace transform (with \( \mathcal{L}\{\mathbf{x}(t)\} = \mathbf{X}(s) \)) gives

\[ s\mathbf{X}(s) - \mathbf{x}_0 = \mathbf{A}\mathbf{X}(s) \quad\Rightarrow\quad (s\mathbf{I} - \mathbf{A})\mathbf{X}(s) = \mathbf{x}_0. \]

Thus,

\[ \mathbf{X}(s) = (s\mathbf{I} - \mathbf{A})^{-1} \mathbf{x}_0. \]

The poles of the entries of \( \mathbf{X}(s) \) occur at values of \( s \) where \( (s\mathbf{I} - \mathbf{A}) \) is singular, i.e. when

\[ \det(s\mathbf{I} - \mathbf{A}) = 0. \]

This is exactly the characteristic equation for eigenvalues of \( \mathbf{A} \), with \( s = \lambda \). Therefore, eigenvalues of the system matrix coincide with poles of the Laplace-domain representation of the homogeneous dynamics.

In later chapters, we will exploit this equivalence to move freely between time-domain matrix analysis and frequency-domain pole analysis.

For the homogeneous system \( \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} \), the equilibrium at \( \mathbf{x} = \mathbf{0} \) is asymptotically stable if \( \mathbf{x}(t) \to \mathbf{0} \) as \( t \to \infty \) for all initial conditions. Using the modal decomposition:

\[ \mathbf{x}(t) = \sum_{i=1}^n \alpha_i e^{\lambda_i t} \mathbf{v}_i. \]

For diagonalizable \( \mathbf{A} \), it is clear that \( \mathbf{x}(t) \to \mathbf{0} \) if and only if \( e^{\lambda_i t} \to 0 \) for all \( i \), which holds exactly when

\[ \operatorname{Re}(\lambda_i) < 0 \quad \text{for all } i. \]

If any eigenvalue has \( \operatorname{Re}(\lambda_i) > 0 \), there exists a component that grows exponentially, and the equilibrium is unstable. If \( \operatorname{Re}(\lambda_i) = 0 \) and the corresponding Jordan block is larger than \( 1 \), polynomial terms in \( t \) appear, leading to instability.

These eigenvalue-based stability conditions are the multi-variable analogue of the scalar condition that a pole must lie in the left half-plane for stability, already encountered in earlier lessons.

Not every matrix is diagonalizable. When geometric multiplicities are smaller than algebraic multiplicities, the best we can do is transform \( \mathbf{A} \) into its Jordan normal form:

\[ \mathbf{A} = \mathbf{S}\mathbf{J}\mathbf{S}^{-1}, \quad \mathbf{J} = \begin{bmatrix} J_1 & 0 & \cdots & 0 \\ 0 & J_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & J_r \end{bmatrix}, \]

where each \( J_k \) is a Jordan block of the form

\[ J_k = \begin{bmatrix} \lambda & 1 & 0 & \cdots & 0 \\ 0 & \lambda & 1 & \cdots & 0 \\ \vdots & & \ddots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda & 1 \\ 0 & 0 & \cdots & 0 & \lambda \end{bmatrix}. \]

For a Jordan block of size \( m \) with eigenvalue \( \lambda \), the contribution to the solution contains terms of the form

\[ t^k e^{\lambda t}, \quad k = 0,1,\dots,m-1. \]

Thus even when \( \operatorname{Re}(\lambda) = 0 \), higher-order Jordan blocks can lead to unbounded growth (through the powers of \( t \)). The diagonalizable case (all Jordan blocks of size 1) is therefore especially favorable in control engineering.

In practice, eigenvalues and eigenvectors are computed numerically. Modern libraries (for Python, C++, Java, MATLAB, Mathematica) implement variants of the QR algorithm and related methods. Two important numerical issues for dynamic systems are:

• Conditioning of eigenvalues: small perturbations of the matrix coefficients can lead to changes in eigenvalues; ill-conditioned matrices may exhibit large sensitivity.
• Scaling: choosing appropriate units of variables and time can significantly improve numerical behavior of eigenvalue routines and time integration algorithms.

When simulating or designing controllers, it is good practice to verify eigenvalues using more than one library or precision level if the system is close to marginal stability (eigenvalues near the imaginary axis).

We illustrate eigenvalue computation and modal solution for the 2D system

\[ \frac{d}{dt} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \mathbf{A} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad \mathbf{A} = \begin{bmatrix} 0 & 1 \\ -2 & -3 \end{bmatrix}, \quad \mathbf{x}(0) = \begin{bmatrix} 1 \\ 0 \end{bmatrix}. \]

This system can be interpreted as a mass–spring–damper model (with suitable units). We compute eigenvalues and the state trajectory using different languages.

10.1 Python (NumPy / SciPy)

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

# System matrix
A = np.array([[0.0, 1.0],
              [-2.0, -3.0]])

x0 = np.array([1.0, 0.0])  # initial state

# Eigenvalues and eigenvectors
lam, V = eig(A)
print("Eigenvalues:", lam)
print("Eigenvectors (columns of V):\n", V)

# Matrix exponential solution at time t
t = 2.0
Phi = expm(A * t)               # Phi = e^{A t}
x_t = Phi @ x0
print("State at t = 2:", x_t)

# Time-domain simulation using ODE solver
def f(t, x):
    return A @ x

t_span = (0.0, 5.0)
t_eval = np.linspace(0.0, 5.0, 200)
sol = solve_ivp(f, t_span, x0, t_eval=t_eval)

# sol.t contains time points, sol.y contains states
x1 = sol.y[0, :]
x2 = sol.y[1, :]
print("Final state (t = 5):", sol.y[:, -1])
      

The eigenvalues of \( \mathbf{A} \) are complex numbers with negative real part, indicating a stable oscillatory mode (underdamped second-order system).

10.2 C++ (Eigen Library)

#include <iostream>
#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>  // for matrix exponential

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

    Eigen::Vector2d x0(1.0, 0.0);

    // Eigenvalue decomposition
    Eigen::EigenSolver<Eigen::Matrix2d> es(A);
    std::cout << "Eigenvalues:\n" << es.eigenvalues() << std::endl;
    std::cout << "Eigenvectors (columns):\n" << es.eigenvectors() << std::endl;

    // Matrix exponential at t = 2
    double t = 2.0;
    Eigen::Matrix2d Phi = (A * t).exp();  // e^{A t}
    Eigen::Vector2d x_t = Phi * x0;
    std::cout << "State at t = 2:\n" << x_t << std::endl;

    // Simple explicit Euler discretization (for illustration)
    double h = 0.01;
    int N = static_cast<int>(5.0 / h);  // simulate to t = 5
    Eigen::Vector2d x = x0;
    for (int k = 0; k < N; ++k) {
        x = x + h * (A * x);
    }
    std::cout << "Approximate state at t = 5 (Euler):\n" << x << std::endl;

    return 0;
}
      

The Eigen library provides high-level linear algebra operations suitable for real-time control implementations where dynamic matrices and eigenvalue computations are needed.

10.3 Java (EJML)

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.dense.row.decomposition.eig.SwitchingEigenDecomposition_DDRM;

public class LinearSystem2D {
    public static void main(String[] args) {
        // System matrix A
        DMatrixRMaj A = new DMatrixRMaj(new double[][] {
            { 0.0,  1.0 },
            { -2.0, -3.0 }
        });

        // Eigenvalue decomposition
        SwitchingEigenDecomposition_DDRM eig =
            new SwitchingEigenDecomposition_DDRM(2);
        eig.decompose(A);

        for (int i = 0; i < 2; ++i) {
            System.out.println("Eigenvalue " + i + ": " + eig.getEigenvalue(i));
            DMatrixRMaj v = eig.getEigenVector(i);
            System.out.println("Eigenvector " + i + ": ");
            v.print();
        }

        // Explicit Euler integration
        double h = 0.01;
        int N = (int)(5.0 / h);
        DMatrixRMaj x = new DMatrixRMaj(2, 1);
        x.set(0, 0, 1.0);  // x1(0)
        x.set(1, 0, 0.0);  // x2(0)

        DMatrixRMaj Ax = new DMatrixRMaj(2, 1);
        for (int k = 0; k < N; ++k) {
            CommonOps_DDRM.mult(A, x, Ax);
            // x = x + h * A * x
            CommonOps_DDRM.addEquals(x, h, Ax);
        }

        System.out.println("Approximate state at t = 5 (Euler):");
        x.print();
    }
}
      

EJML (Efficient Java Matrix Library) is a convenient choice for implementing linear system simulations and eigenvalue computations in Java-based control software.

10.4 MATLAB / Simulink

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

x0 = [1; 0];

% Eigenvalues and eigenvectors
[V, D] = eig(A);      % columns of V are eigenvectors, D is diagonal of eigenvalues

disp('Eigenvalues:'), disp(diag(D))
disp('Eigenvectors (columns):'), disp(V)

% Matrix exponential solution at t = 2
t = 2;
Phi = expm(A*t);
x_t = Phi * x0;
disp('State at t = 2:'), disp(x_t)

% Time response using ODE solver
f = @(t, x) A*x;
[t_sim, x_sim] = ode45(f, [0 5], x0);

% Plot
figure; plot(t_sim, x_sim(:,1), 'LineWidth', 1.5); hold on;
plot(t_sim, x_sim(:,2), 'LineWidth', 1.5);
xlabel('t'); ylabel('states'); legend('x_1', 'x_2'); grid on;
title('Response of 2D linear system');

% Simulink note (conceptual):
% In Simulink, you can build an equivalent model by:
%  - Using an Integrator block for x_dot = A x
%  - Feeding the state vector x into a Matrix Gain block with gain A
%  - Closing the loop so that the Matrix Gain output drives the Integrator input
%  - Initial condition of the Integrator is x0
      

MATLAB has built-in high-level support for matrix exponentials and eigenvalues. Simulink provides block-diagram based simulation of the underlying vector ODE without manually coding the integrator.

10.5 Wolfram Mathematica

(* System matrix and initial state *)
A = { {0, 1},
     {-2, -3} };
x0 = {1, 0};

(* Eigenvalues and eigenvectors *)
evals = Eigenvalues[A];
evecs = Eigenvectors[A];
Print["Eigenvalues: ", evals];
Print["Eigenvectors: ", evecs];

(* Matrix exponential solution at t = 2 *)
t = 2;
Phi = MatrixExp[A t];
xT = Phi . x0;
Print["State at t = 2: ", xT];

(* Symbolic solution x(t) = exp(A t) x0 *)
Clear[tvar];
x[tvar_] := MatrixExp[A tvar].x0;
Simplify[x[tvar]]

(* Numerical time simulation *)
sol = NDSolve[{x'[t] == A.x[t], x[0] == x0}, x, {t, 0, 5}];
Plot[Evaluate[{x[t][[1]], x[t][[2]]} /. sol],
     {t, 0, 5}, PlotLegends -> {"x1", "x2"},
     AxesLabel -> {"t", "state"}]
      

Mathematica combines symbolic and numeric linear algebra, making it useful for deriving closed-form expressions for \( e^{\mathbf{A}t} \) and for exploring parameter dependencies analytically.

Problem 1 (Eigenvalues and Modal Solution):

Consider the system

\[ \frac{d}{dt} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ -4 & -4 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad \mathbf{x}(0) = \begin{bmatrix} 1 \\ 0 \end{bmatrix}. \]

(a) Compute the eigenvalues and eigenvectors of the system matrix. (b) Derive the explicit solution \( \mathbf{x}(t) \). (c) Determine whether the equilibrium is asymptotically stable.

Solution:

(a) The characteristic polynomial is

\[ p(\lambda) = \det \begin{bmatrix} -\lambda & 1 \\ -4 & -4-\lambda \end{bmatrix} = \lambda(4+\lambda)+4 = \lambda^2 + 4\lambda + 4 = (\lambda+2)^2. \]

Hence there is a repeated eigenvalue \( \lambda = -2 \) of algebraic multiplicity 2. The eigenvectors satisfy

\[ (\mathbf{A} + 2\mathbf{I})\mathbf{v} = \mathbf{0} \quad\Rightarrow\quad \begin{bmatrix} 2 & 1 \\ -4 & -2 \end{bmatrix} \begin{bmatrix} v_1 \\ v_2 \end{bmatrix} = \mathbf{0}. \]

The rows are multiples, so we have \( 2 v_1 + v_2 = 0 \), i.e., \( v_2 = -2 v_1 \). A convenient eigenvector is \( \mathbf{v}_1 = [1\ -2]^\top \). The eigenspace is one-dimensional, so \( \mathbf{A} \) is not diagonalizable.

(b) For a defective eigenvalue with a Jordan block of size 2, the general solution has the form

\[ \mathbf{x}(t) = e^{-2 t}\left( \alpha_1 \mathbf{v}_1 + \alpha_2 (\mathbf{v}_1 t + \mathbf{w}) \right), \]

where \( \mathbf{w} \) is a generalized eigenvector satisfying \( (\mathbf{A}+2\mathbf{I})\mathbf{w} = \mathbf{v}_1 \). One possible construction:

\[ (\mathbf{A}+2\mathbf{I})\mathbf{w} = \begin{bmatrix} 2 & 1 \\ -4 & -2 \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \end{bmatrix} = \begin{bmatrix} 1 \\ -2 \end{bmatrix}. \]

Solving the first row \( 2 w_1 + w_2 = 1 \), choose, for example, \( w_1 = 0 \) and \( w_2 = 1 \). Then \( \mathbf{w} = [0\ 1]^\top \) satisfies both rows.

The initial condition \( \mathbf{x}(0) = \mathbf{x}_0 = [1\ 0]^\top \) implies

\[ \mathbf{x}_0 = \alpha_1 \mathbf{v}_1 + \alpha_2 \mathbf{w}. \]

So

\[ \alpha_1 \begin{bmatrix} 1 \\ -2 \end{bmatrix} + \alpha_2 \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \end{bmatrix} \quad\Rightarrow\quad \alpha_1 = 1,\; -2\alpha_1 + \alpha_2 = 0 \Rightarrow \alpha_2 = 2. \]

Thus

\[ \mathbf{x}(t) = e^{-2 t}\left( \mathbf{v}_1 + 2(\mathbf{v}_1 t + \mathbf{w}) \right) = e^{-2 t} \left( (1+2t) \begin{bmatrix} 1 \\ -2 \end{bmatrix} + 2 \begin{bmatrix} 0 \\ 1 \end{bmatrix} \right). \]

(c) Despite the defective eigenvalue, \( \operatorname{Re}(\lambda) = -2 < 0 \), and the polynomial factor in \( t \) is multiplied by \( e^{-2 t} \), so \( \mathbf{x}(t) \to \mathbf{0} \) as \( t \to \infty \). The equilibrium is asymptotically stable.

Problem 2 (Diagonalization and Matrix Exponential):

Let \( \mathbf{A} \) be diagonalizable with \( \mathbf{A} = \mathbf{S}\boldsymbol{\Lambda}\mathbf{S}^{-1} \). Prove that \( e^{\mathbf{A}t} = \mathbf{S}e^{\boldsymbol{\Lambda}t}\mathbf{S}^{-1} \).

Solution:

By definition of the matrix exponential,

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

Since \( \mathbf{A} = \mathbf{S}\boldsymbol{\Lambda}\mathbf{S}^{-1} \), by induction \( \mathbf{A}^k = \mathbf{S}\boldsymbol{\Lambda}^k\mathbf{S}^{-1} \). Hence

\[ e^{\mathbf{A}t} = \sum_{k=0}^{\infty} \frac{t^k}{k!} \mathbf{S}\boldsymbol{\Lambda}^k\mathbf{S}^{-1} = \mathbf{S} \left( \sum_{k=0}^{\infty} \frac{t^k}{k!} \boldsymbol{\Lambda}^k \right) \mathbf{S}^{-1} = \mathbf{S} e^{\boldsymbol{\Lambda} t} \mathbf{S}^{-1}, \]

where the series for \( e^{\boldsymbol{\Lambda} t} \) reduces to exponentials of diagonal entries. This proves the identity.

Problem 3 (Stability Condition):

Let \( \mathbf{A} \in \mathbb{R}^{n \times n} \) be diagonalizable with eigenvalues \( \lambda_i \). Prove that \( \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} \) is asymptotically stable if and only if \( \operatorname{Re}(\lambda_i) < 0 \) for all \( i \).

Solution:

Using the modal decomposition,

\[ \mathbf{x}(t) = \mathbf{S} e^{\boldsymbol{\Lambda} t} \mathbf{S}^{-1} \mathbf{x}_0 = \sum_{i=1}^n \alpha_i e^{\lambda_i t} \mathbf{v}_i, \]

where \( \mathbf{v}_i \) are eigenvectors and \( \alpha_i \) depend on \( \mathbf{x}_0 \). As \( t \to \infty \), each term behaves like \( e^{\lambda_i t} \). If \( \operatorname{Re}(\lambda_i) < 0 \), then \( e^{\lambda_i t} \to 0 \), so all modes decay and \( \mathbf{x}(t) \to \mathbf{0} \). Conversely, if there exists \( j \) with \( \operatorname{Re}(\lambda_j) > 0 \), then for an initial condition with \( \alpha_j \neq 0 \), the corresponding mode grows unboundedly, so the equilibrium is unstable. If \( \operatorname{Re}(\lambda_j) = 0 \), the mode does not decay, violating asymptotic stability. Hence the condition is necessary and sufficient.

Problem 4 (Workflow for Analyzing a System Matrix):

Outline a systematic workflow to analyze a given system matrix \( \mathbf{A} \) from a dynamic system model, indicating steps to determine stability and modal structure.

Solution (conceptual flow):

flowchart TD
  S["Start with system matrix A"] --> C["Compute eigenvalues lambda_i"]
  C --> G["Check Re(lambda_i)"]
  G -->|all < 0| STABLE["Classify as \nasymptotically stable"]
  G -->|"some > 0"| UNSTABLE["Classify as unstable"]
  G -->|"some = 0"| BORDER["Borderline - inspect Jordan form"]
  BORDER --> J["Compute geometric multiplicities / \nJordan blocks"]
  J --> D["If any block size > 1, unstable; \nelse marginally stable"]
        

In practice, numerical routines approximate eigenvalues and optionally estimate condition numbers to assess sensitivity. For borderline cases, further analysis (e.g. Lyapunov methods) is recommended, which will appear in later chapters.

In this lesson, we connected linear algebra to dynamic systems modeling. Vectors and matrices provide a compact representation for coupled first-order ODEs describing multi-variable systems. Eigenvalues and eigenvectors characterize the natural modes of these systems, with diagonalization yielding a modal decomposition in which each mode evolves as \( e^{\lambda_i t} \). We showed how the matrix exponential solves homogeneous linear systems, how eigenvalues coincide with poles from the Laplace-domain viewpoint, and how stability is governed by the real parts of eigenvalues. Finally, we implemented these concepts in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica, preparing the computational foundation for more advanced modeling and control design in later chapters.