In the previous lesson, stability of the homogeneous LTI system \( \dot{\mathbf{x} }(t)=A\mathbf{x}(t) \) was characterized by the eigenvalues of \( A \). This lesson goes one level deeper: it asks how each eigenvalue appears in the actual state trajectory and output response. A mode is a dynamical component of the free response associated with an eigenvalue and its eigenvector direction.

\[ \mathbf{x}(t)=e^{At}\mathbf{x}(0),\qquad \mathbf{y}(t)=C e^{At}\mathbf{x}(0) \]

If \( A \) is diagonalizable, the response is a sum of elementary terms of the form \( e^{\lambda_i t}\mathbf{v}_i \). Thus, eigenvalues do not only decide stability; they also decide the shape, speed, oscillation, and relative importance of the transient response.

flowchart TD
  A["State matrix A"] --> B["Compute eigenvalues and eigenvectors"]
  B --> C["Transform state: x = V z"]
  C --> D["Decoupled modal equations: z_i_dot = lambda_i z_i"]
  D --> E["Modal response: z_i(t) = exp(lambda_i t) z_i(0)"]
  E --> F["Reconstruct state and output"]
  F --> G["Identify slow, fast, oscillatory, and dominant modes"]
        

Assume \( A\in\mathbb{R}^{n\times n} \) has \( n \) linearly independent eigenvectors \( \mathbf{v}_1,\ldots,\mathbf{v}_n \). Define

\[ V=\begin{bmatrix}\mathbf{v}_1 & \mathbf{v}_2 & \cdots & \mathbf{v}_n\end{bmatrix},\qquad \Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n) \]

Since \( AV=V\Lambda \), we have the similarity relation \( A=V\Lambda V^{-1} \). Using the modal coordinate vector \( \mathbf{z}=V^{-1}\mathbf{x} \), the state equation becomes

\[ \dot{\mathbf{z} }=V^{-1}\dot{\mathbf{x} } =V^{-1}A\mathbf{x} =V^{-1}AV\mathbf{z} =\Lambda\mathbf{z}. \]

Therefore, each scalar modal coordinate satisfies an independent first-order differential equation:

\[ \dot{z}_i(t)=\lambda_i z_i(t),\qquad z_i(t)=e^{\lambda_i t}z_i(0). \]

Reconstructing the physical state gives the modal expansion

\[ \mathbf{x}(t)=V e^{\Lambda t}V^{-1}\mathbf{x}(0) =\sum_{i=1}^n \mathbf{v}_i e^{\lambda_i t} \left(\mathbf{w}_i^T\mathbf{x}(0)\right), \]

where \( \mathbf{w}_i^T \) is the \( i \)-th row of \( V^{-1} \). The scalar \( \mathbf{w}_i^T\mathbf{x}(0) \) measures how strongly the initial condition excites the \( i \)-th mode.

Proof: matrix exponential under a similarity transformation

The series definition of the matrix exponential gives

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

Since \( A=V\Lambda V^{-1} \), repeated multiplication yields \( A^k=V\Lambda^kV^{-1} \). Thus

\[ e^{At} =\sum_{k=0}^\infty \frac{V\Lambda^kV^{-1}t^k}{k!} =V\left(\sum_{k=0}^\infty \frac{(\Lambda t)^k}{k!}\right)V^{-1} =Ve^{\Lambda t}V^{-1}. \]

State modes are not always equally visible in a measured output. For \( \mathbf{y}=C\mathbf{x} \), substituting the modal expansion gives

\[ \mathbf{y}(t) =\sum_{i=1}^n \left(C\mathbf{v}_i\right) e^{\lambda_i t} \left(\mathbf{w}_i^T\mathbf{x}(0)\right). \]

The factor \( C\mathbf{v}_i \) is the output direction of the mode. The factor \( \mathbf{w}_i^T\mathbf{x}(0) \) is its excitation by the initial condition. Their product is the modal contribution seen in the output.

\[ \text{output modal amplitude}_i = \left(C\mathbf{v}_i\right) \left(\mathbf{w}_i^T\mathbf{x}(0)\right). \]

A mode can be internally present but absent from a particular output if \( C\mathbf{v}_i=\mathbf{0} \). A mode can also be present in the system but not excited by a specific initial condition if \( \mathbf{w}_i^T\mathbf{x}(0)=0 \). This is a key reason why dominant modes are not identified by eigenvalues alone; they also depend on modal amplitudes and on the output being examined.

For a real state matrix, complex eigenvalues occur in conjugate pairs. If \( \lambda=\sigma+j\omega \) and \( \bar{\lambda}=\sigma-j\omega \), their combined real response has an exponentially scaled sinusoidal form:

\[ e^{\sigma t} \left[ \mathbf{a}\cos(\omega t)+ \mathbf{b}\sin(\omega t) \right]. \]

The real part \( \sigma \) controls the envelope: if \( \sigma < 0 \), the oscillation decays; if \( \sigma=0 \), the oscillation has constant envelope; if \( \sigma > 0 \), the oscillation grows. The imaginary part \( \omega \) controls the damped angular frequency.

\[ T_d=\frac{2\pi}{|\omega|},\qquad \text{oscillation frequency}=\frac{|\omega|}{2\pi}. \]

For a real eigenvalue \( \lambda \), the associated elementary term is non-oscillatory: \( e^{\lambda t}\mathbf{v} \). A negative real eigenvalue gives monotone decay along the corresponding eigenvector direction; a positive real eigenvalue gives monotone growth.

The dominant modes of a stable system are usually the modes whose eigenvalues have real parts closest to the imaginary axis. Define the spectral abscissa:

\[ \alpha(A)=\max_i \operatorname{Re}(\lambda_i). \]

If the system is asymptotically stable, then \( \alpha(A)<0 \). The modes with \( \operatorname{Re}(\lambda_i)=\alpha(A) \) decay slowest. Fast stable modes have more negative real parts and disappear earlier in the response.

Suppose \( \lambda_d \) is a dominant eigenvalue and \( \operatorname{Re}(\lambda_j)<\operatorname{Re}(\lambda_d) \). For nonzero coefficients, the magnitude ratio of the two exponential envelopes satisfies

\[ \frac{|e^{\lambda_j t}|}{|e^{\lambda_d t}|} = e^{\left(\operatorname{Re}(\lambda_j) -\operatorname{Re}(\lambda_d)\right)t}. \]

Since \( \operatorname{Re}(\lambda_j)-\operatorname{Re}(\lambda_d)<0 \), this ratio approaches zero for large time. Thus, after the fast transients have decayed, the output is often well approximated by the slowest visible and excited modes.

flowchart TD
  S["Start with eigenvalues and modal amplitudes"] --> A["Discard modes with zero excitation"]
  A --> B["Discard modes invisible in selected output"]
  B --> C["Sort remaining modes by real part"]
  C --> D["Largest real part gives slowest decay or fastest growth"]
  D --> E["Check complex pairs for oscillatory dominance"]
  E --> F["Use dominant modes for approximate transient model"]
        

For a stable complex dominant pair \( \lambda_{1,2}=\sigma\pm j\omega \), a common engineering approximation is

\[ T_s\approx \frac{4}{|\sigma|},\qquad T_d=\frac{2\pi}{|\omega|}, \]

where \( T_s \) is the approximate two-percent settling time based on the exponential envelope and \( T_d \) is the damped oscillation period.

If \( A \) is not diagonalizable, modal decomposition uses Jordan blocks. Detailed Jordan-chain construction is developed later in the course, but the key dynamical fact is already important for stability. For a Jordan block

\[ J=\lambda I+N,\qquad N^m=0, \]

where \( N \) is nilpotent. Since \( \lambda I \) and \( N \) commute,

\[ e^{Jt}=e^{(\lambda I+N)t} =e^{\lambda t}e^{Nt} = e^{\lambda t} \left(I+Nt+\frac{N^2t^2}{2!}+\cdots+ \frac{N^{m-1}t^{m-1} }{(m-1)!}\right). \]

Therefore, defective repeated eigenvalues produce terms such as \( t e^{\lambda t} \), \( t^2 e^{\lambda t} \), and higher polynomially weighted exponentials. For continuous-time stability, these polynomial factors are harmless if \( \operatorname{Re}(\lambda)<0 \), but they destroy Lyapunov stability if \( \operatorname{Re}(\lambda)=0 \) and the corresponding Jordan block has size larger than one.

If \( A=V\Lambda V^{-1} \), then

\[ \|e^{At}\| = \|Ve^{\Lambda t}V^{-1}\| \leq \|V\|\,\|V^{-1}\|\,e^{\alpha(A)t}. \]

The factor \( \kappa(V)=\|V\|\|V^{-1}\| \) is the condition number of the eigenvector matrix. If \( V \) is ill-conditioned, stable eigenvalues may coexist with large transient amplification before eventual decay. This is common in nonnormal systems, where eigenvectors are far from orthogonal. Thus, eigenvalues determine asymptotic stability, but modal geometry influences finite-time behavior.

In practical control analysis, this warning matters because a system may satisfy the eigenvalue stability criterion yet still exhibit large temporary excursions, actuator saturation, or large internal state magnitudes before the asymptotic decay becomes visible.

Modern-control software typically computes modes using dense or sparse eigenvalue routines. For small and medium systems, direct eigenvalue decompositions are sufficient. For large systems, Krylov and Arnoldi methods are often used to compute only the dominant eigenvalues.

• Python: numpy.linalg.eig, scipy.linalg.eig, scipy.linalg.expm, scipy.sparse.linalg.eigs, and the optional python-control package.
• C++: Eigen, Armadillo, Blaze, and LAPACK bindings for eigenvalue and matrix-exponential computations.
• Java: EJML, ojAlgo, and Apache Commons Math for numerical linear algebra.
• MATLAB/Simulink: eig, expm, ss, initial, and the State-Space block in Simulink.
• Wolfram Mathematica: Eigensystem, MatrixExp, JordanDecomposition, and symbolic simplification for exact modal expressions.

The Python program computes eigenvalues, left modal rows, modal amplitudes, full matrix-exponential response, and a dominant-mode approximation.

Chapter9_Lesson3.py

# Chapter9_Lesson3.py
# Modes, modal decomposition, and dominant modes for continuous-time LTI systems.
#
# Required libraries:
#   pip install numpy scipy matplotlib
#
# Optional control-engineering library:
#   pip install control
#
# This file studies x_dot = A x, y = C x by comparing:
#   1) full matrix-exponential response,
#   2) exact modal reconstruction using eigenvectors,
#   3) dominant-mode approximation.

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


def modal_decomposition(A: np.ndarray):
    """Return eigenvalues, right eigenvectors, and left modal rows.

    If A is diagonalizable, A = V Lambda V^{-1}.
    The rows of W = V^{-1} are left modal coordinate functionals.
    Modal coordinate z_i(0) = w_i x(0).
    """
    lambdas, V = eig(A)
    W = np.linalg.inv(V)
    return lambdas, V, W


def full_response(A: np.ndarray, x0: np.ndarray, t_grid: np.ndarray) -> np.ndarray:
    """Compute x(t)=exp(A t)x0 for each t."""
    return np.array([expm(A * t) @ x0 for t in t_grid])


def modal_response(lambdas: np.ndarray, V: np.ndarray, W: np.ndarray,
                   x0: np.ndarray, t_grid: np.ndarray,
                   selected_modes=None) -> np.ndarray:
    """Reconstruct the state from selected modes.

    selected_modes=None uses all modes. For real systems with complex-conjugate
    pairs, the final real part is the physical state.
    """
    n = len(lambdas)
    if selected_modes is None:
        selected_modes = list(range(n))

    z0 = W @ x0
    xs = []
    for t in t_grid:
        x = np.zeros(V.shape[0], dtype=complex)
        for i in selected_modes:
            x += V[:, i] * np.exp(lambdas[i] * t) * z0[i]
        xs.append(np.real_if_close(x).real)
    return np.array(xs)


def dominant_modes_by_real_part(lambdas: np.ndarray, keep: int):
    """Return indices of modes with largest real parts."""
    return [int(i) for i in np.argsort(np.real(lambdas))[::-1][:keep]]


def modal_table(lambdas: np.ndarray, V: np.ndarray, W: np.ndarray,
                C: np.ndarray, x0: np.ndarray):
    """Build a readable modal table containing excitation and output residue."""
    z0 = W @ x0
    rows = []
    for i, lam in enumerate(lambdas):
        output_residue = C @ V[:, i] * z0[i]
        rows.append({
            "mode": i,
            "lambda": lam,
            "real_part": np.real(lam),
            "imag_part": np.imag(lam),
            "initial_modal_amplitude": z0[i],
            "output_residue": output_residue[0],
        })
    return rows


def main():
    # A stable third-order system:
    #   lambda_1,2 = -0.25 +/- 1.50 j  (slow oscillatory dominant pair)
    #   lambda_3   = -3.00             (fast real mode)
    A = np.array([
        [-0.25,  1.50, 0.0],
        [-1.50, -0.25, 0.0],
        [ 0.00,  0.00, -3.0],
    ], dtype=float)

    C = np.array([[1.0, 0.0, 0.4]])
    x0 = np.array([1.0, -0.2, 2.0])
    t_grid = np.linspace(0.0, 20.0, 800)

    lambdas, V, W = modal_decomposition(A)

    print("Eigenvalues:")
    for lam in lambdas:
        print(f"  {lam.real:+.4f} {lam.imag:+.4f}j")

    print("\nModal participation/output residue table:")
    for row in modal_table(lambdas, V, W, C, x0):
        print(
            f"mode={row['mode']}, lambda={row['lambda']:.4g}, "
            f"z0={row['initial_modal_amplitude']:.4g}, "
            f"output residue={row['output_residue']:.4g}"
        )

    x_full = full_response(A, x0, t_grid)
    x_modal = modal_response(lambdas, V, W, x0, t_grid)

    dominant = dominant_modes_by_real_part(lambdas, keep=2)
    x_dom = modal_response(lambdas, V, W, x0, t_grid, selected_modes=dominant)

    y_full = (C @ x_full.T).ravel()
    y_modal = (C @ x_modal.T).ravel()
    y_dom = (C @ x_dom.T).ravel()

    print("\nMax reconstruction error between expm and all-mode reconstruction:")
    print(np.max(np.abs(y_full - y_modal)))

    print("\nDominant mode indices:", dominant)
    print("Dominant eigenvalues:", lambdas[dominant])

    plt.figure(figsize=(9, 5))
    plt.plot(t_grid, y_full, label="full output y(t)")
    plt.plot(t_grid, y_dom, "--", label="dominant-mode approximation")
    plt.xlabel("time [s]")
    plt.ylabel("output")
    plt.title("Modal decomposition and dominant-mode approximation")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.show()

    plt.figure(figsize=(9, 5))
    plt.semilogy(t_grid, np.abs(y_full - y_dom) + 1e-14)
    plt.xlabel("time [s]")
    plt.ylabel("|full - dominant|")
    plt.title("Error caused by neglecting fast modes")
    plt.grid(True)
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

The C++ implementation uses the Eigen library. It computes the complex eigen-decomposition and reconstructs the output from all modes and from only the dominant modes.

Chapter9_Lesson3.cpp

// Chapter9_Lesson3.cpp
// Modes, modal decomposition, and dominant modes for x_dot = A x.
//
// Dependency: Eigen C++ template library
// Ubuntu/Debian: sudo apt install libeigen3-dev
// Compile:
//   g++ -std=c++17 Chapter9_Lesson3.cpp -I /usr/include/eigen3 -O2 -o Chapter9_Lesson3
//
// The program compares full modal reconstruction with a dominant-mode
// approximation for a third-order continuous-time LTI system.

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <complex>
#include <iomanip>
#include <iostream>
#include <vector>
#include <algorithm>

using Complex = std::complex<double>;
using CMatrix = Eigen::MatrixXcd;
using CVector = Eigen::VectorXcd;

CVector modalState(
    const Eigen::VectorXcd& lambda,
    const CMatrix& V,
    const CMatrix& W,
    const CVector& x0,
    double t,
    const std::vector<int>& selectedModes
) {
    CVector z0 = W * x0;
    CVector x = CVector::Zero(x0.size());

    for (int idx : selectedModes) {
        x += V.col(idx) * std::exp(lambda(idx) * t) * z0(idx);
    }
    return x;
}

int main() {
    Eigen::Matrix3d A;
    A << -0.25,  1.50, 0.0,
         -1.50, -0.25, 0.0,
          0.00,  0.00, -3.0;

    Eigen::RowVector3d C;
    C << 1.0, 0.0, 0.4;

    Eigen::Vector3d x0real;
    x0real << 1.0, -0.2, 2.0;

    Eigen::ComplexEigenSolver<Eigen::Matrix3d> solver(A);
    Eigen::VectorXcd lambda = solver.eigenvalues();
    CMatrix V = solver.eigenvectors();
    CMatrix W = V.inverse();

    CVector x0 = x0real.cast<Complex>();

    std::cout << "Eigenvalues:\n";
    for (int i = 0; i < lambda.size(); ++i) {
        std::cout << "  lambda_" << i << " = "
                  << lambda(i).real() << " + "
                  << lambda(i).imag() << "j\n";
    }

    std::vector<int> allModes = {0, 1, 2};

    std::vector<int> dominant = allModes;
    std::sort(dominant.begin(), dominant.end(), [&](int a, int b) {
        return lambda(a).real() > lambda(b).real();
    });
    dominant.resize(2);

    std::cout << "\nDominant mode indices:";
    for (int idx : dominant) {
        std::cout << " " << idx;
    }
    std::cout << "\n\n";

    Eigen::RowVectorXcd Cc = C.cast<Complex>();

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "t, y_full, y_dominant, absolute_error\n";

    for (int k = 0; k <= 40; ++k) {
        double t = 0.5 * k;

        CVector xFull = modalState(lambda, V, W, x0, t, allModes);
        CVector xDom = modalState(lambda, V, W, x0, t, dominant);

        Complex yFull = (Cc * xFull)(0);
        Complex yDom = (Cc * xDom)(0);
        double error = std::abs(yFull.real() - yDom.real());

        std::cout << t << ", "
                  << yFull.real() << ", "
                  << yDom.real() << ", "
                  << error << "\n";
    }

    return 0;
}

This Java implementation is dependency-free and uses an analytically known modal structure: one stable complex-conjugate pair and one fast real mode. For general matrices, use EJML, ojAlgo, or Apache Commons Math.

Chapter9_Lesson3.java

// Chapter9_Lesson3.java
// Analytic modal decomposition for a real second-order oscillatory mode
// plus one fast real mode.
//
// Compile:
//   javac Chapter9_Lesson3.java
// Run:
//   java Chapter9_Lesson3
//
// For general eigenvalue computations in Java, libraries such as EJML,
// ojAlgo, or Apache Commons Math can be used. This file stays dependency-free
// by using a system whose modes are known analytically.

public class Chapter9_Lesson3 {
    static class State {
        double x1;
        double x2;
        double x3;

        State(double x1, double x2, double x3) {
            this.x1 = x1;
            this.x2 = x2;
            this.x3 = x3;
        }
    }

    static State fullAnalyticResponse(double sigma, double omega, double fastPole,
                                      State x0, double t) {
        double eSlow = Math.exp(sigma * t);
        double c = Math.cos(omega * t);
        double s = Math.sin(omega * t);

        // Dynamics:
        // x1_dot = sigma*x1 + omega*x2
        // x2_dot = -omega*x1 + sigma*x2
        double x1 = eSlow * (x0.x1 * c + x0.x2 * s);
        double x2 = eSlow * (-x0.x1 * s + x0.x2 * c);
        double x3 = Math.exp(fastPole * t) * x0.x3;

        return new State(x1, x2, x3);
    }

    static State dominantApproximation(double sigma, double omega, State x0, double t) {
        double eSlow = Math.exp(sigma * t);
        double c = Math.cos(omega * t);
        double s = Math.sin(omega * t);

        double x1 = eSlow * (x0.x1 * c + x0.x2 * s);
        double x2 = eSlow * (-x0.x1 * s + x0.x2 * c);

        // The fast state is neglected in the dominant-mode model.
        return new State(x1, x2, 0.0);
    }

    static double output(State x) {
        return x.x1 + 0.4 * x.x3;
    }

    public static void main(String[] args) {
        double sigma = -0.25;
        double omega = 1.50;
        double fastPole = -3.00;

        State x0 = new State(1.0, -0.2, 2.0);

        System.out.println("Eigenvalues:");
        System.out.println("  lambda_1,2 = " + sigma + " +/- " + omega + "j");
        System.out.println("  lambda_3   = " + fastPole);
        System.out.println();
        System.out.println("t, y_full, y_dominant, absolute_error");

        for (int k = 0; k <= 40; k++) {
            double t = 0.5 * k;

            State full = fullAnalyticResponse(sigma, omega, fastPole, x0, t);
            State dom = dominantApproximation(sigma, omega, x0, t);

            double yFull = output(full);
            double yDom = output(dom);
            double error = Math.abs(yFull - yDom);

            System.out.printf("%.2f, %.8f, %.8f, %.8e%n", t, yFull, yDom, error);
        }

        double settlingTimeApprox = 4.0 / Math.abs(sigma);
        double dampedPeriod = 2.0 * Math.PI / omega;

        System.out.println();
        System.out.printf("Dominant-mode settling-time approximation: %.4f s%n",
                          settlingTimeApprox);
        System.out.printf("Dominant oscillation period: %.4f s%n", dampedPeriod);
    }
}

The MATLAB script computes modal coordinates and uses Control System Toolbox functions. It also includes an optional Simulink model creation block using a State-Space block.

Chapter9_Lesson3.m

% Chapter9_Lesson3.m
% Modes, modal decomposition, and dominant modes for x_dot = A x, y = C x.
%
% MATLAB functions used:
%   eig, expm, ss, initial
%
% Optional Simulink section:
%   creates a minimal State-Space block model if Simulink is available.

clear; clc; close all;

A = [-0.25  1.50  0.0;
     -1.50 -0.25  0.0;
      0.00  0.00 -3.0];

B = zeros(3,1);
C = [1.0 0.0 0.4];
D = 0;

x0 = [1.0; -0.2; 2.0];
t = linspace(0,20,800);

[V,Lambda] = eig(A);
lambda = diag(Lambda);
W = inv(V);
z0 = W*x0;

disp('Eigenvalues:');
disp(lambda);

% Sort modes by largest real part.
[~,idx] = sort(real(lambda),'descend');
dominant = idx(1:2);

Xfull = zeros(length(t),3);
Xmodal = zeros(length(t),3);
Xdom = zeros(length(t),3);

for k = 1:length(t)
    tk = t(k);

    Xfull(k,:) = (expm(A*tk)*x0).';

    xModal = V*diag(exp(lambda*tk))*z0;
    Xmodal(k,:) = real(xModal).';

    xDom = zeros(3,1);
    for m = 1:length(dominant)
        i = dominant(m);
        xDom = xDom + V(:,i)*exp(lambda(i)*tk)*z0(i);
    end
    Xdom(k,:) = real(xDom).';
end

yFull = (C*Xfull.').';
yModal = (C*Xmodal.').';
yDom = (C*Xdom.').';

fprintf('Max all-mode reconstruction error: %.3e\n', max(abs(yFull-yModal)));
fprintf('Dominant modes: ');
fprintf('%d ', dominant);
fprintf('\n');

figure;
plot(t,yFull,'LineWidth',1.5); hold on;
plot(t,yDom,'--','LineWidth',1.5);
grid on;
xlabel('time [s]');
ylabel('output');
title('Full response and dominant-mode approximation');
legend('full y(t)','dominant-mode approximation');

figure;
semilogy(t,abs(yFull-yDom)+1e-14,'LineWidth',1.5);
grid on;
xlabel('time [s]');
ylabel('|full - dominant|');
title('Error caused by neglected fast mode');

% Control System Toolbox representation.
sys = ss(A,B,C,D);
[yInitial,tInitial,xInitial] = initial(sys,x0,t);

figure;
plot(tInitial,yInitial,'LineWidth',1.5);
grid on;
xlabel('time [s]');
ylabel('output');
title('Initial-condition response using ss/initial');

% Optional: create a simple Simulink State-Space model.
% This section is safe to skip on systems without Simulink.
if exist('new_system','file') == 2
    model = 'Chapter9_Lesson3_Simulink';
    if bdIsLoaded(model)
        close_system(model,0);
    end

    new_system(model);
    open_system(model);

    add_block('simulink/Sources/Constant',[model '/zero input'],...
        'Value','0','Position',[60 80 110 110]);

    add_block('simulink/Continuous/State-Space',[model '/State-Space'],...
        'A','[-0.25 1.5 0; -1.5 -0.25 0; 0 0 -3]',...
        'B','[0;0;0]',...
        'C','[1 0 0.4]',...
        'D','0',...
        'X0','[1;-0.2;2]',...
        'Position',[180 65 300 125]);

    add_block('simulink/Sinks/Scope',[model '/Scope'],...
        'Position',[390 75 430 115]);

    add_line(model,'zero input/1','State-Space/1');
    add_line(model,'State-Space/1','Scope/1');

    set_param(model,'StopTime','20');
    save_system(model);
    fprintf('Simulink model created: %s.slx\n',model);
end

Mathematica is useful for exact modal formulas, symbolic eigen-analysis, and closed-form matrix exponentials.

Chapter9_Lesson3.nb

(* Chapter9_Lesson3.nb *)
(* Wolfram Mathematica / Wolfram Language code for modal decomposition
   and dominant-mode approximation. Paste into a Mathematica notebook or
   open this file as plain text in Mathematica. *)

ClearAll["Global`*"];

A = { {-0.25, 1.50, 0.0},
     {-1.50, -0.25, 0.0},
     {0.0, 0.0, -3.0} };

Cmat = { {1.0, 0.0, 0.4} };
x0 = {1.0, -0.2, 2.0};

{vals, vecsRaw} = Eigensystem[A];
V = Transpose[vecsRaw];
W = Inverse[V];
z0 = W.x0;

Print["Eigenvalues:"];
Print[N[vals]];

dominant = Take[Ordering[Re[vals], All, Greater], 2];
Print["Dominant mode indices: ", dominant];

xFull[t_] := MatrixExp[A t].x0;

xModal[t_] := Re[V.DiagonalMatrix[Exp[vals t]].z0];

xDominant[t_] := Module[{x = ConstantArray[0, 3]},
   Do[
    x = x + V[[All, i]] Exp[vals[[i]] t] z0[[i]],
    {i, dominant}
   ];
   Re[x]
];

yFull[t_] := First[Cmat.xFull[t]];
yModal[t_] := First[Cmat.xModal[t]];
yDominant[t_] := First[Cmat.xDominant[t]];

Print["All-mode reconstruction error at t=5:"];
Print[N[Abs[yFull[5] - yModal[5]]]];

Plot[
 {yFull[t], yDominant[t]},
 {t, 0, 20},
 PlotLegends -> {"full output", "dominant-mode approximation"},
 AxesLabel -> {"time [s]", "output"},
 PlotLabel -> "Modal decomposition and dominant modes",
 PlotRange -> All
]

Plot[
 Abs[yFull[t] - yDominant[t]] + 10^-14,
 {t, 0, 20},
 ScalingFunctions -> "Log",
 AxesLabel -> {"time [s]", "|full - dominant|"},
 PlotLabel -> "Error from neglected fast mode",
 PlotRange -> All
]

Problem 1 (Diagonalizable Modal Expansion): Let \( A\in\mathbb{R}^{n\times n} \) have \( n \) linearly independent eigenvectors. Prove that the homogeneous solution can be written as \( \mathbf{x}(t)=\sum_{i=1}^n c_i e^{\lambda_i t}\mathbf{v}_i \).

Solution: Since the eigenvectors form a basis, write \( \mathbf{x}(0)=\sum_{i=1}^n c_i\mathbf{v}_i \). The matrix exponential is linear, so

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

Since \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \), the exponential maps each eigenvector direction to itself:

\[ e^{At}\mathbf{v}_i = \sum_{k=0}^\infty \frac{A^k t^k}{k!}\mathbf{v}_i = \sum_{k=0}^\infty \frac{\lambda_i^k t^k}{k!}\mathbf{v}_i = e^{\lambda_i t}\mathbf{v}_i. \]

Substitution gives the desired modal expansion.

Problem 2 (Dominant Stable Mode): A stable system has eigenvalues \( -5 \), \( -2 \), and \( -0.4 \). Which mode is dominant and why?

Solution: The dominant eigenvalue is \( -0.4 \) because it has the largest real part. Its exponential envelope \( e^{-0.4t} \) decays more slowly than \( e^{-2t} \) and \( e^{-5t} \). For example,

\[ \frac{e^{-2t} }{e^{-0.4t} }=e^{-1.6t},\qquad \frac{e^{-5t} }{e^{-0.4t} }=e^{-4.6t}. \]

Both ratios decay to zero, so the \( -0.4 \) mode remains visible longest if it is excited and observable.

Problem 3 (Complex Dominant Pair): A system has dominant eigenvalues \( \lambda_{1,2}=-0.2\pm 3j \). Estimate the approximate settling time and oscillation period.

Solution: The exponential envelope is \( e^{-0.2t} \), so the two-percent settling-time approximation gives

\[ T_s\approx \frac{4}{0.2}=20. \]

The damped oscillation period is

\[ T_d=\frac{2\pi}{3}\approx 2.094. \]

Problem 4 (Invisible Mode in the Output): Suppose \( A\mathbf{v}_i=\lambda_i\mathbf{v}_i \) and \( C\mathbf{v}_i=\mathbf{0} \). What is the contribution of this mode to \( \mathbf{y}(t)=C\mathbf{x}(t) \)?

Solution: The output contribution of mode \( i \) is

\[ \left(C\mathbf{v}_i\right)e^{\lambda_i t} \left(\mathbf{w}_i^T\mathbf{x}(0)\right). \]

If \( C\mathbf{v}_i=\mathbf{0} \), this term is identically zero for every initial condition. The mode may still exist internally, but it is not visible in the selected output channel.

Problem 5 (Jordan Block and Stability): Consider \( J=\begin{bmatrix}\lambda & 1\\0 & \lambda\end{bmatrix} \). Compute \( e^{Jt} \) and explain what happens when \( \lambda=0 \).

Solution: Write \( J=\lambda I+N \), where \( N=\begin{bmatrix}0 & 1\\0 & 0\end{bmatrix} \) and \( N^2=0 \). Therefore,

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

If \( \lambda=0 \), then \( e^{Jt}=\begin{bmatrix}1 & t\\0 & 1\end{bmatrix} \), which is unbounded as time increases. Thus, a repeated eigenvalue on the imaginary axis with a nontrivial Jordan block is not Lyapunov stable.

Modal decomposition expresses the free response of an LTI system as a sum of elementary dynamical components. For diagonalizable systems, modal coordinates decouple the dynamics into scalar equations \( \dot{z}_i=\lambda_i z_i \). Dominant modes are the modes that remain most important in the response, usually because their eigenvalues have real parts closest to the imaginary axis and because their modal amplitudes are excited and visible in the output. Repeated defective eigenvalues introduce polynomially weighted exponential terms, while ill-conditioned eigenvectors can cause transient amplification even when all eigenvalues are stable.