Consider the continuous-time LTI system \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \), \( \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \), where \( \mathbf{x}\in\mathbb{C}^n \), \( \mathbf{u}\in\mathbb{C}^m \), and \( \mathbf{y}\in\mathbb{C}^p \). A modal form is a coordinate representation in which the homogeneous state dynamics are separated into independent scalar modes. The central assumption in this lesson is that \( \mathbf{A} \) has \( n \) distinct eigenvalues \( \lambda_1,\dots,\lambda_n \).

\[ \mathbf{A}\mathbf{v}_i=\lambda_i\mathbf{v}_i,\quad \lambda_i\neq\lambda_j\;\text{for}\;i\neq j,\quad \mathbf{V}=\begin{bmatrix}\mathbf{v}_1&\mathbf{v}_2&\cdots&\mathbf{v}_n\end{bmatrix}. \]

If \( \mathbf{V} \) is nonsingular, define modal coordinates by \( \mathbf{x}=\mathbf{V}\mathbf{z} \). Then each component \( z_i \) is the amplitude of the \( i \)-th eigenvector direction in the state.

flowchart TD
  A["Start with physical model: xdot = A x + B u, y = C x + D u"] --> B["Compute distinct eigenvalues lambda_i"]
  B --> C["Compute right eigenvectors v_i"]
  C --> D["Build V = [v1 v2 ... vn]"]
  D --> E["Change coordinates: x = V z"]
  E --> F["Modal model: zdot = Lambda z + Bm u"]
  F --> G["Output model: y = Cm z + D u"]
  G --> H["Interpret each scalar mode separately"]
        

The diagonal modal form rests on a fundamental linear algebra theorem: eigenvectors associated with distinct eigenvalues are linearly independent. Therefore, a matrix with \( n \) distinct eigenvalues has an eigenvector basis.

\[ \mathbf{A}\mathbf{V} =\begin{bmatrix}\mathbf{A}\mathbf{v}_1&\cdots&\mathbf{A}\mathbf{v}_n\end{bmatrix} =\begin{bmatrix}\lambda_1\mathbf{v}_1&\cdots&\lambda_n\mathbf{v}_n\end{bmatrix} =\mathbf{V}\boldsymbol{\Lambda}, \]

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

Proof of independence. Suppose \( \sum_{i=1}^r \alpha_i\mathbf{v}_i=\mathbf{0} \) is a shortest nontrivial dependence among eigenvectors associated with distinct eigenvalues. Multiplying by \( \mathbf{A} \) and also multiplying the original relation by \( \lambda_r \) gives

\[ \sum_{i=1}^r \alpha_i\lambda_i\mathbf{v}_i=\mathbf{0},\qquad \sum_{i=1}^r \alpha_i\lambda_r\mathbf{v}_i=\mathbf{0}. \]

Subtracting the two equations eliminates the last vector:

\[ \sum_{i=1}^{r-1} \alpha_i(\lambda_i-\lambda_r)\mathbf{v}_i= \mathbf{0}. \]

Since \( \lambda_i-\lambda_r\neq 0 \) for \( i=1,\dots,r-1 \), this is a shorter nontrivial dependence, contradicting the minimality of the original dependence. Hence the eigenvectors are linearly independent and \( \mathbf{V} \) is invertible.

With \( \mathbf{x}=\mathbf{V}\mathbf{z} \) and constant \( \mathbf{V} \), differentiation gives \( \dot{\mathbf{x} }=\mathbf{V}\dot{\mathbf{z} } \). Substitute this into the physical-coordinate state equation:

\[ \mathbf{V}\dot{\mathbf{z} }=\mathbf{A}\mathbf{V}\mathbf{z}+ \mathbf{B}\mathbf{u}. \]

Premultiplying by \( \mathbf{V}^{-1} \) yields

\[ \dot{\mathbf{z} }=\underbrace{\mathbf{V}^{-1}\mathbf{A}\mathbf{V} }_{\boldsymbol{\Lambda} } \mathbf{z}+\underbrace{\mathbf{V}^{-1}\mathbf{B} }_{\mathbf{B}_m} \mathbf{u}. \]

The output equation becomes

\[ \mathbf{y}=\underbrace{\mathbf{C}\mathbf{V} }_{\mathbf{C}_m} \mathbf{z}+\mathbf{D}\mathbf{u}. \]

Therefore, the modal realization is

\[ \boxed{\dot{\mathbf{z} }=\boldsymbol{\Lambda}\mathbf{z}+ \mathbf{B}_m\mathbf{u},\qquad \mathbf{y}=\mathbf{C}_m\mathbf{z}+\mathbf{D}\mathbf{u}.} \]

Unlike controllable and observable companion forms, the diagonal modal form exposes the eigenvalues directly on the diagonal. The cost is that \( \mathbf{B}_m \) and \( \mathbf{C}_m \) may become dense and may be numerically sensitive if the eigenvector matrix is ill-conditioned.

Since \( \boldsymbol{\Lambda} \) is diagonal, the modal state equation separates into \( n \) scalar first-order differential equations. If \( \mathbf{b}_{m,i} \) denotes the \( i \)-th row of \( \mathbf{B}_m \), then

\[ \dot{z}_i(t)=\lambda_i z_i(t)+\mathbf{b}_{m,i}\mathbf{u}(t), \qquad i=1,\dots,n. \]

For the homogeneous response \( \mathbf{u}(t)=\mathbf{0} \), each modal coordinate evolves independently:

\[ z_i(t)=e^{\lambda_i(t-t_0)}z_i(t_0),\qquad \mathbf{z}(t)=e^{\boldsymbol{\Lambda}(t-t_0)}\mathbf{z}(t_0). \]

\[ e^{\boldsymbol{\Lambda}t}=\operatorname{diag}(e^{\lambda_1t}, e^{\lambda_2t},\dots,e^{\lambda_nt}),\qquad e^{\mathbf{A}t}=\mathbf{V}e^{\boldsymbol{\Lambda}t}\mathbf{V}^{-1}. \]

For a general input, variation of constants gives

\[ z_i(t)=e^{\lambda_i(t-t_0)}z_i(t_0)+ \int_{t_0}^t e^{\lambda_i(t-\tau)}\mathbf{b}_{m,i} \mathbf{u}(\tau)\,d\tau. \]

Thus modal form converts an \( n \)-state coupled vector equation into \( n \) forced scalar equations. This is the main conceptual reason modal coordinates are so useful for interpretation.

flowchart LR
  U["Input u"] --> Z1["mode 1: zdot1 = lambda1 z1 + row1(Bm) u"]
  U --> Z2["mode 2: zdot2 = lambda2 z2 + row2(Bm) u"]
  U --> ZN["mode n: zdotn = lambdan zn + rown(Bm) u"]
  Z1 --> Y["Output y = Cm z + D u"]
  Z2 --> Y
  ZN --> Y
        

Let \( \mathbf{p}_i^T \) be the \( i \)-th row of \( \mathbf{V}^{-1} \). Then \( \mathbf{p}_i^T \) is a left eigenvector of \( \mathbf{A} \) and satisfies the biorthogonality relation with the right eigenvectors:

\[ \mathbf{p}_i^T\mathbf{A}=\lambda_i\mathbf{p}_i^T, \qquad \mathbf{p}_i^T\mathbf{v}_j=\delta_{ij}. \]

This gives an explicit formula for modal coordinates:

\[ z_i=\mathbf{p}_i^T\mathbf{x},\qquad \mathbf{x}=\sum_{i=1}^n \mathbf{v}_i z_i. \]

The modal input and output directions can also be written mode by mode:

\[ \mathbf{B}_m=\begin{bmatrix}\mathbf{p}_1^T\mathbf{B}\\ \mathbf{p}_2^T\mathbf{B}\\ \vdots\\ \mathbf{p}_n^T\mathbf{B} \end{bmatrix},\qquad \mathbf{C}_m=\begin{bmatrix}\mathbf{C}\mathbf{v}_1& \mathbf{C}\mathbf{v}_2&\cdots&\mathbf{C}\mathbf{v}_n \end{bmatrix}. \]

Hence \( \mathbf{p}_i^T\mathbf{B} \) measures how strongly the inputs excite mode \( i \), while \( \mathbf{C}\mathbf{v}_i \) measures how strongly mode \( i \) appears in the measured output.

Because similarity transformations preserve controllability and observability, the pairs \( (\mathbf{A},\mathbf{B}) \) and \( (\boldsymbol{\Lambda},\mathbf{B}_m) \) have the same controllability status. Likewise, \( (\mathbf{A},\mathbf{C}) \) and \( (\boldsymbol{\Lambda},\mathbf{C}_m) \) have the same observability status.

For distinct eigenvalues, the PBH tests take an especially transparent modal form. Mode \( i \) is controllable precisely when its modal input row is nonzero:

\[ \operatorname{rank}\begin{bmatrix}\lambda_i\mathbf{I}- \mathbf{A}&\mathbf{B}\end{bmatrix}=n \quad\Longleftrightarrow\quad \mathbf{p}_i^T\mathbf{B}\neq\mathbf{0}. \]

Similarly, mode \( i \) is observable precisely when its modal output column is nonzero:

\[ \operatorname{rank}\begin{bmatrix}\lambda_i\mathbf{I}- \mathbf{A}\\ \mathbf{C}\end{bmatrix}=n \quad\Longleftrightarrow\quad \mathbf{C}\mathbf{v}_i\neq\mathbf{0}. \]

Thus, in diagonal modal form, an uncontrollable mode is one whose row in \( \mathbf{B}_m \) is zero, and an unobservable mode is one whose column in \( \mathbf{C}_m \) is zero. This gives a direct modal reading of the rank tests developed in earlier chapters.

The transfer matrix is invariant under similarity transformation:

\[ \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+ \mathbf{D}=\mathbf{C}_m(s\mathbf{I}-\boldsymbol{\Lambda})^{-1} \mathbf{B}_m+\mathbf{D}. \]

Since \( \boldsymbol{\Lambda} \) is diagonal,

\[ (s\mathbf{I}-\boldsymbol{\Lambda})^{-1}= \operatorname{diag}\left(\frac{1}{s-\lambda_1},\frac{1}{s-\lambda_2}, \dots,\frac{1}{s-\lambda_n}\right). \]

Therefore the transfer matrix decomposes into modal residues:

\[ \boxed{\mathbf{G}(s)=\mathbf{D}+ \sum_{i=1}^n \frac{(\mathbf{C}\mathbf{v}_i)(\mathbf{p}_i^T\mathbf{B})}{s-\lambda_i}.} \]

For a SISO system, the residue of pole \( \lambda_i \) is the scalar

\[ r_i=(\mathbf{C}\mathbf{v}_i)(\mathbf{p}_i^T\mathbf{B}). \]

A zero residue may indicate that the pole is hidden from the input-output transfer function because the mode is not excited by the input, not seen by the output, or canceled by the product of input and output directions. This point connects modal form to the minimal-realization ideas that will be developed later.

If \( \mathbf{A} \) is real and has a nonreal eigenvalue \( \lambda=\alpha+j\omega \), its conjugate \( \bar{\lambda}=\alpha-j\omega \) is also an eigenvalue. Over \( \mathbb{C} \), the modal block is diagonal. Over \( \mathbb{R} \), the corresponding real modal block is the two-dimensional rotation-decay block

\[ \mathbf{A}_{\alpha,\omega}=\begin{bmatrix}\alpha&\omega\\ -\omega&\alpha\end{bmatrix}. \]

This real block gives the response components

\[ e^{\mathbf{A}_{\alpha,\omega}t}=e^{\alpha t} \begin{bmatrix}\cos(\omega t)&\sin(\omega t)\\ -\sin(\omega t)&\cos(\omega t)\end{bmatrix}. \]

In this lesson the diagonal form is developed over the complex field for mathematical clarity. In numerical real-valued software, modal routines often return real block-diagonal forms rather than a fully complex diagonal form.

Although diagonal modal form is mathematically elegant, it may be a poor computational coordinate system when eigenvectors are nearly linearly dependent. The relevant diagnostic is the condition number of \( \mathbf{V} \):

\[ \kappa(\mathbf{V})=\|\mathbf{V}\|\,\|\mathbf{V}^{-1}\|. \]

A small perturbation in \( \mathbf{A} \) may cause a large perturbation in modal vectors when \( \kappa(\mathbf{V}) \) is large. At the level of the transformation,

\[ \frac{\|\Delta \mathbf{z}\|}{\|\mathbf{z}\|} \lesssim \kappa(\mathbf{V})\, \frac{\|\Delta \mathbf{x}\|}{\|\mathbf{x}\|},\qquad \mathbf{z}=\mathbf{V}^{-1}\mathbf{x}. \]

Therefore modal form is excellent for theoretical interpretation and low-order hand analysis, but Schur-based or real block modal forms are often preferred in robust numerical computations.

The Python implementation uses numpy and scipy to compute the eigenvector matrix, transform \( \mathbf{B} \) and \( \mathbf{C} \), and verify that physical-coordinate and modal-coordinate zero-input simulations agree.

Chapter17_Lesson3.py

"""
Chapter17_Lesson3.py
Diagonal modal form for a continuous-time LTI system with distinct eigenvalues.

Libraries:
    numpy      : matrix and vector operations
    scipy      : matrix exponential and linear solves
Optional:
    python-control: higher-level state-space workflows, not required here
"""

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


def modal_form(A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray):
    """Return modal matrices Lambda, Bm, Cm, D, and eigenvector matrix V.

    Coordinate convention:
        x = V z,    z = V^{-1} x
        zdot = Lambda z + Bm u
        y    = Cm z + D u
    """
    eigenvalues, V = eig(A)
    Vinv = inv(V)
    Lambda = Vinv @ A @ V
    Bm = Vinv @ B
    Cm = C @ V
    return eigenvalues, V, Vinv, Lambda, Bm, Cm, D


def simulate_zero_input(A: np.ndarray, C: np.ndarray, x0: np.ndarray, t_grid):
    """Compute x(t) and y(t) for u(t)=0 using the matrix exponential."""
    xs, ys = [], []
    for t in t_grid:
        x_t = expm(A * t) @ x0
        y_t = C @ x_t
        xs.append(x_t)
        ys.append(y_t)
    return np.array(xs), np.array(ys)


def modal_zero_input(eigenvalues, V, Vinv, C, x0, t_grid):
    """Compute zero-input response from modal coordinates."""
    z0 = Vinv @ x0
    xs, ys = [], []
    for t in t_grid:
        exp_diag = np.diag(np.exp(eigenvalues * t))
        z_t = exp_diag @ z0
        x_t = V @ z_t
        y_t = C @ x_t
        xs.append(x_t)
        ys.append(y_t)
    return np.array(xs), np.array(ys)


def main():
    # A has distinct eigenvalues -1, -2, -4 and is diagonalizable.
    A = np.array([[0.0, 1.0, 0.0],
                  [-2.0, -3.0, 0.0],
                  [0.5, 0.0, -4.0]])
    B = np.array([[0.0],
                  [1.0],
                  [0.0]])
    C = np.array([[1.0, 0.0, 1.0]])
    D = np.array([[0.0]])
    x0 = np.array([1.0, 0.0, -0.5])

    lam, V, Vinv, Lambda, Bm, Cm, Dm = modal_form(A, B, C, D)

    print("Eigenvalues:")
    print(lam)
    print("\nEigenvector matrix V:")
    print(V)
    print("\nModal A matrix Lambda = V^{-1} A V:")
    print(np.real_if_close(Lambda))
    print("\nModal input matrix Bm = V^{-1} B:")
    print(np.real_if_close(Bm))
    print("\nModal output matrix Cm = C V:")
    print(np.real_if_close(Cm))
    print("\nCondition number of V:", np.linalg.cond(V))

    t_grid = np.linspace(0.0, 5.0, 6)
    x_direct, y_direct = simulate_zero_input(A, C, x0, t_grid)
    x_modal, y_modal = modal_zero_input(lam, V, Vinv, C, x0, t_grid)

    print("\nAgreement check ||x_direct - x_modal||:", norm(x_direct - x_modal))
    print("Agreement check ||y_direct - y_modal||:", norm(y_direct - y_modal))
    print("\nSample output response y(t):")
    for t, y in zip(t_grid, y_modal):
        print(f"t={t:4.2f}, y={np.real_if_close(y[0]): .6f}")


if __name__ == "__main__":
    main()

The C++ implementation uses the Eigen library. It computes complex eigenvectors even when the example happens to have real eigenvalues, making the structure closer to the general modal algorithm.

Chapter17_Lesson3.cpp

/*
Chapter17_Lesson3.cpp
Diagonal modal form for a continuous-time LTI system with distinct real eigenvalues.

Required library:
    Eigen 3  (https://eigen.tuxfamily.org)
Compile example:
    g++ Chapter17_Lesson3.cpp -std=c++17 -I /path/to/eigen -O2 -o modal_form
*/

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

int main() {
    using Matrix = Eigen::MatrixXd;
    using CMatrix = Eigen::MatrixXcd;
    using CVector = Eigen::VectorXcd;

    Matrix A(3, 3);
    A << 0.0, 1.0, 0.0,
        -2.0, -3.0, 0.0,
         0.5, 0.0, -4.0;

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

    Matrix C(1, 3);
    C << 1.0, 0.0, 1.0;

    Eigen::EigenSolver<Matrix> solver(A);
    CVector lambda = solver.eigenvalues();
    CMatrix V = solver.eigenvectors();
    CMatrix Vinv = V.inverse();
    CMatrix Lambda = Vinv * A.cast<std::complex<double>>() * V;
    CMatrix Bm = Vinv * B.cast<std::complex<double>>();
    CMatrix Cm = C.cast<std::complex<double>>() * V;

    std::cout << "Eigenvalues:\n" << lambda << "\n\n";
    std::cout << "Modal A matrix Lambda = V^{-1} A V:\n" << Lambda << "\n\n";
    std::cout << "Modal input matrix Bm = V^{-1} B:\n" << Bm << "\n\n";
    std::cout << "Modal output matrix Cm = C V:\n" << Cm << "\n\n";

    CVector x0(3);
    x0 << 1.0, 0.0, -0.5;
    CVector z0 = Vinv * x0;

    std::cout << "Zero-input modal response y(t):\n";
    for (double t = 0.0; t <= 5.0; t += 1.0) {
        CMatrix expLambda = CMatrix::Zero(3, 3);
        for (int i = 0; i < 3; ++i) {
            expLambda(i, i) = std::exp(lambda(i) * t);
        }
        CVector zt = expLambda * z0;
        CVector xt = V * zt;
        std::complex<double> yt = (C.cast<std::complex<double>>() * xt)(0);
        std::cout << "t=" << t << ", y=" << yt << "\n";
    }

    return 0;
}

The Java implementation uses Apache Commons Math. The chosen example has distinct real eigenvalues, so the implementation can stay in real arithmetic. For complex modes in production code, use a numerical linear algebra library that exposes complex Schur or complex eigenvector workflows.

Chapter17_Lesson3.java

/*
Chapter17_Lesson3.java
Diagonal modal form for a continuous-time LTI system with distinct real eigenvalues.

Required library:
    Apache Commons Math 3.x
Compile example:
    javac -cp commons-math3-3.6.1.jar Chapter17_Lesson3.java
Run example:
    java -cp .:commons-math3-3.6.1.jar Chapter17_Lesson3
*/

import org.apache.commons.math3.linear.Array2DRowRealMatrix;
import org.apache.commons.math3.linear.EigenDecomposition;
import org.apache.commons.math3.linear.LUDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;

public class Chapter17_Lesson3 {
    public static void main(String[] args) {
        double[][] aData = {
            {0.0, 1.0, 0.0},
            {-2.0, -3.0, 0.0},
            {0.5, 0.0, -4.0}
        };
        double[][] bData = {
            {0.0},
            {1.0},
            {0.0}
        };
        double[][] cData = {
            {1.0, 0.0, 1.0}
        };

        RealMatrix A = new Array2DRowRealMatrix(aData);
        RealMatrix B = new Array2DRowRealMatrix(bData);
        RealMatrix C = new Array2DRowRealMatrix(cData);

        EigenDecomposition ed = new EigenDecomposition(A);
        RealMatrix V = ed.getV();
        RealMatrix Vinv = new LUDecomposition(V).getSolver().getInverse();
        RealMatrix Lambda = Vinv.multiply(A).multiply(V);
        RealMatrix Bm = Vinv.multiply(B);
        RealMatrix Cm = C.multiply(V);

        System.out.println("Eigenvalues:");
        for (int i = 0; i < A.getRowDimension(); i++) {
            System.out.printf("lambda_%d = %.8f%n", i + 1, ed.getRealEigenvalue(i));
        }

        printMatrix("Modal A matrix Lambda = V^{-1} A V", Lambda);
        printMatrix("Modal input matrix Bm = V^{-1} B", Bm);
        printMatrix("Modal output matrix Cm = C V", Cm);

        double[] x0Data = {1.0, 0.0, -0.5};
        RealVector x0 = new Array2DRowRealMatrix(new double[][]{ {1.0}, {0.0}, {-0.5} }).getColumnVector(0);
        RealVector z0 = Vinv.operate(x0);

        System.out.println("Zero-input modal response y(t):");
        for (int step = 0; step <= 5; step++) {
            double t = step;
            double[][] expData = new double[3][3];
            for (int i = 0; i < 3; i++) {
                expData[i][i] = Math.exp(ed.getRealEigenvalue(i) * t);
            }
            RealMatrix expLambda = new Array2DRowRealMatrix(expData);
            RealVector zt = expLambda.operate(z0);
            RealVector xt = V.operate(zt);
            double y = C.operate(xt).getEntry(0);
            System.out.printf("t=%.1f, y=%.8f%n", t, y);
        }
    }

    private static void printMatrix(String title, RealMatrix M) {
        System.out.println("\n" + title + ":");
        for (int i = 0; i < M.getRowDimension(); i++) {
            for (int j = 0; j < M.getColumnDimension(); j++) {
                System.out.printf("%12.6f ", M.getEntry(i, j));
            }
            System.out.println();
        }
    }
}

MATLAB can compute the diagonal modal coordinates directly using eig. With Control System Toolbox, modalreal can compute modal realizations and may return real block forms for complex conjugate modes. In Simulink, each modal equation can be represented by an integrator with feedback gain \( \lambda_i \) and input gain from the corresponding row of \( \mathbf{B}_m \).

Chapter17_Lesson3.m

% Chapter17_Lesson3.m
% Diagonal modal form for a continuous-time LTI system with distinct eigenvalues.
%
% MATLAB libraries/toolboxes:
%   eig, inv, expm are available in base MATLAB.
%   ss, initial, modalreal require Control System Toolbox.
%   Simulink implementation can use a State-Space block with A,B,C,D or
%   separate Integrator/Gain/Sum blocks for each modal state.

clear; clc;

A = [ 0.0  1.0  0.0;
     -2.0 -3.0  0.0;
      0.5  0.0 -4.0];
B = [0.0; 1.0; 0.0];
C = [1.0 0.0 1.0];
D = 0.0;
x0 = [1.0; 0.0; -0.5];

% Modal transformation x = V z, z = inv(V) x
[V,Lambda] = eig(A);
Vinv = inv(V);
Bm = Vinv * B;
Cm = C * V;
z0 = Vinv * x0;

disp('Eigenvalues:');
disp(diag(Lambda));
disp('Modal A matrix Lambda = inv(V)*A*V:');
disp(Vinv * A * V);
disp('Modal input matrix Bm = inv(V)*B:');
disp(Bm);
disp('Modal output matrix Cm = C*V:');
disp(Cm);
disp('Condition number of V:');
disp(cond(V));

% Zero-input response from physical and modal coordinates
t = linspace(0,5,101);
yPhysical = zeros(size(t));
yModal = zeros(size(t));
for k = 1:length(t)
    xPhysical = expm(A*t(k)) * x0;
    zModal = expm(Lambda*t(k)) * z0;
    xFromModal = V * zModal;
    yPhysical(k) = C * xPhysical;
    yModal(k) = C * xFromModal;
end

fprintf('Agreement norm between physical and modal outputs: %.3e\n', norm(yPhysical-yModal));

% Control System Toolbox workflow
sys = ss(A,B,C,D);
try
    [sysModal, blocks] = modalreal(sys);
    disp('modalreal block information:');
    disp(blocks);
    disp('Modal realization returned by modalreal:');
    sysModal
catch ME
    disp('modalreal is unavailable or Control System Toolbox is missing.');
    disp(ME.message);
end

% Simulink note:
% To build the modal model manually, create one Integrator block for each z_i:
%   dz_i/dt = lambda_i z_i + row_i(Bm) u
% Then compute y = Cm z + D u using Gain and Sum blocks.
% For complex conjugate eigenvalues, use real 2-by-2 modal blocks instead.

Mathematica is convenient for symbolic verification of \( \mathbf{V}^{-1}\mathbf{A}\mathbf{V}=\boldsymbol{\Lambda} \) and for expressing modal output responses using MatrixExp.

Chapter17_Lesson3.nb

(* Chapter17_Lesson3.nb *)
A = { {0, 1, 0}, {-2, -3, 0}, {0.5, 0, -4} };
B = { {0}, {1}, {0} };
Cmat = { {1, 0, 1} };
Dmat = { {0} };

{eVals, eVecs} = Eigensystem[A];
V = Transpose[eVecs];
Vinv = Inverse[V];
Lambda = Simplify[Vinv . A . V];
Bm = Simplify[Vinv . B];
Cm = Simplify[Cmat . V];

x0 = {1, 0, -0.5};
z0 = Vinv . x0;
y[t_] = Simplify[Cmat . V . MatrixExp[Lambda t] . z0];

Print["Eigenvalues = ", eVals];
Print["Lambda = ", MatrixForm[Lambda]];
Print["Bm = ", MatrixForm[Bm]];
Print["Cm = ", MatrixForm[Cm]];
Table[{tt, N[y[tt]]}, {tt, 0, 5}]

Problem 1 (Constructing Modal Form): Consider \( \mathbf{A}=\begin{bmatrix}0&1\\-2&-3\end{bmatrix} \), \( \mathbf{B}=\begin{bmatrix}0\\1\end{bmatrix} \), \( \mathbf{C}=\begin{bmatrix}1&0\end{bmatrix} \), and \( D=0 \). Find a diagonal modal realization.

Solution: The characteristic polynomial is

\[ \det(s\mathbf{I}-\mathbf{A})=s^2+3s+2=(s+1)(s+2). \]

Hence \( \lambda_1=-1 \) and \( \lambda_2=-2 \). Eigenvectors may be chosen as

\[ \mathbf{v}_1=\begin{bmatrix}1\\-1\end{bmatrix},\qquad \mathbf{v}_2=\begin{bmatrix}1\\-2\end{bmatrix},\qquad \mathbf{V}=\begin{bmatrix}1&1\\-1&-2\end{bmatrix}. \]

\[ \mathbf{V}^{-1}=\begin{bmatrix}2&1\\-1&-1\end{bmatrix}, \qquad \boldsymbol{\Lambda}=\mathbf{V}^{-1}\mathbf{A}\mathbf{V} =\begin{bmatrix}-1&0\\0&-2\end{bmatrix}. \]

The transformed input and output matrices are

\[ \mathbf{B}_m=\mathbf{V}^{-1}\mathbf{B} =\begin{bmatrix}1\\-1\end{bmatrix},\qquad \mathbf{C}_m=\mathbf{C}\mathbf{V}=\begin{bmatrix}1&1\end{bmatrix}. \]

Problem 2 (Modal Solution): For the system in Problem 1 with zero input and \( \mathbf{x}(0)=\begin{bmatrix}3\\-4\end{bmatrix} \), compute \( \mathbf{x}(t) \) using modal coordinates.

Solution: First compute \( \mathbf{z}(0)=\mathbf{V}^{-1}\mathbf{x}(0) \):

\[ \mathbf{z}(0)=\begin{bmatrix}2&1\\-1&-1\end{bmatrix} \begin{bmatrix}3\\-4\end{bmatrix} =\begin{bmatrix}2\\1\end{bmatrix}. \]

\[ \mathbf{z}(t)=\begin{bmatrix}2e^{-t}\\e^{-2t}\end{bmatrix}, \qquad \mathbf{x}(t)=\mathbf{V}\mathbf{z}(t) =\begin{bmatrix}2e^{-t}+e^{-2t}\\-2e^{-t}-2e^{-2t}\end{bmatrix}. \]

Problem 3 (Modal Observability): For a system with distinct eigenvalues, prove that mode \( i \) is unobservable if \( \mathbf{C}\mathbf{v}_i=\mathbf{0} \).

Solution: In modal coordinates, the zero-input output is

\[ \mathbf{y}(t)=\sum_{k=1}^n (\mathbf{C}\mathbf{v}_k)e^{\lambda_k t}z_k(0). \]

If \( \mathbf{C}\mathbf{v}_i=\mathbf{0} \), then the term involving \( z_i(0) \) is absent from the output for all time. Therefore a change in the initial amplitude along \( \mathbf{v}_i \) cannot be inferred from the output, so the corresponding mode is unobservable.

Problem 4 (Transfer-Function Residues): Show that for a SISO diagonal modal system, the transfer function can be written as a sum of simple fractions.

Solution: With \( \boldsymbol{\Lambda}=\operatorname{diag}(\lambda_1, \dots,\lambda_n) \), \( \mathbf{B}_m=\begin{bmatrix}b_1&\cdots&b_n\end{bmatrix}^T \), and \( \mathbf{C}_m=\begin{bmatrix}c_1&\cdots&c_n\end{bmatrix} \),

\[ G(s)=D+\mathbf{C}_m(s\mathbf{I}-\boldsymbol{\Lambda})^{-1} \mathbf{B}_m =D+\sum_{i=1}^n \frac{c_i b_i}{s-\lambda_i}. \]

The scalar \( c_i b_i \) is the residue associated with pole \( \lambda_i \).

Problem 5 (Conditioning): Explain why a modal transformation can be numerically unreliable even when all eigenvalues are distinct.

Solution: Distinct eigenvalues guarantee diagonalizability but do not guarantee that the eigenvector matrix is well conditioned. If \( \kappa(\mathbf{V}) \) is large, small perturbations in the original state or in the computed eigenvectors can produce large perturbations in modal coordinates because \( \mathbf{z}=\mathbf{V}^{-1}\mathbf{x} \). This is why numerical packages often prefer orthogonal Schur forms or real block modal forms for computation.

For distinct eigenvalues, the state matrix admits an eigenvector basis and can be transformed to diagonal modal form. The transformation \( \mathbf{x}=\mathbf{V}\mathbf{z} \) produces \( \boldsymbol{\Lambda}=\mathbf{V}^{-1}\mathbf{A}\mathbf{V} \), \( \mathbf{B}_m=\mathbf{V}^{-1}\mathbf{B} \), and \( \mathbf{C}_m=\mathbf{C}\mathbf{V} \). Each modal coordinate evolves as a scalar first-order system. Modal form also makes the PBH tests visually transparent: input rows reveal controllable modes and output columns reveal observable modes. The main limitation is numerical conditioning of the eigenvector matrix.