Consider the continuous-time LTI state equation (single or multi-input): \( \dot{\mathbf{x} }(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t) \), with initial condition \( \mathbf{x}(0)=\mathbf{x}_0 \). From Lesson 2, the variation-of-constants formula gives:

\[ \mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0 + \int_{0}^{t} e^{\mathbf{A}(t - \tau)}\,\mathbf{B}\,\mathbf{u}(\tau)\, d\tau. \]

The integral above is a convolution-like operator with kernel \( e^{\mathbf{A}(t-\tau)}\mathbf{B} \). For special classes of inputs—constant, step, and polynomials—we can evaluate this integral in closed form using matrix identities. A key theme is replacing “integrate numerically” with “evaluate a matrix function”.

flowchart TD
  S["Start: xdot = A x + B u(t)"] --> F["Use x(t) = exp(A t) x0 + \nintegral exp(A(t-s)) B u(s) ds"]
  F --> T{"Input type?"}
  T --> C["Constant: \nu(t)=u0"]
  T --> H["Step: \nu(t)=u0 * 1(t-ts)"]
  T --> P["Polynomial: \nu(t)=sum u_k t^k"]
  C --> C1["Compute integral via\n A^{-1}(exp(A t)-I) if possible"]
  C --> C2["Or via phi_1(A t) \nseries when A singular"]
  H --> H1["Piecewise: t < ts \nand t >= ts"]
  P --> P1["Use phi_{k+1}(A t) \nfor each monomial"]
  C1 --> A1["Assemble x(t)"]
  C2 --> A1
  H1 --> A1
  P1 --> A1
  A1 --> E["Validate: differentiate to \nrecover dynamics"]
        

Let \( \mathbf{u}(t)=\mathbf{u}_0 \) (constant vector). Then

\[ \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 + \left(\int_{0}^{t} e^{\mathbf{A}(t-\tau)}\,\mathbf{B}\, d\tau\right)\mathbf{u}_0. \]

Substitute \( s=t-\tau \) (so \( d\tau = -ds \) and limits swap), yielding:

\[ \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\,\mathbf{B}\, d\tau = \int_{0}^{t} e^{\mathbf{A}s}\,\mathbf{B}\, ds. \]

Case 1: \( \mathbf{A} \) is invertible. Define \( \mathbf{F}(t)=\mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I}) \). Then:

\[ \frac{d}{dt}\mathbf{F}(t) = \mathbf{A}^{-1}(\mathbf{A}e^{\mathbf{A}t}) = e^{\mathbf{A}t}, \qquad \mathbf{F}(0)=\mathbf{0}. \]

Since \( \mathbf{F}'(t)=e^{\mathbf{A}t} \) and \( \mathbf{F}(0)=\mathbf{0} \), the fundamental theorem of calculus implies:

\[ \int_{0}^{t} e^{\mathbf{A}s}\, ds = \mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I}). \]

Therefore, the constant-input solution becomes:

\[ \mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0 + \mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I})\,\mathbf{B}\mathbf{u}_0. \]

Verification (direct proof by differentiation). Let \( \mathbf{x}(t) \) be as above. Using \( \frac{d}{dt}e^{\mathbf{A}t}=\mathbf{A}e^{\mathbf{A}t} \) and the identity \( \frac{d}{dt}\left[\mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I})\right]=e^{\mathbf{A}t} \), we obtain:

\[ \dot{\mathbf{x} }(t) = \mathbf{A}e^{\mathbf{A}t}\mathbf{x}_0 + e^{\mathbf{A}t}\mathbf{B}\mathbf{u}_0 = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}_0, \]

and \( \mathbf{x}(0)=\mathbf{x}_0 \) holds because \( e^{\mathbf{A}0}=\mathbf{I} \).

Case 2: \( \mathbf{A} \) may be singular. The expression \( \mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I}) \) may not exist. Instead define the matrix function \( \boldsymbol{\varphi}_1(\mathbf{Z}) \) by the absolutely convergent series:

\[ \boldsymbol{\varphi}_1(\mathbf{Z}) \triangleq \sum_{j=0}^{\infty}\frac{\mathbf{Z}^j}{(j+1)!}. \]

Then:

\[ \int_{0}^{t} e^{\mathbf{A}s}\,ds = t\,\boldsymbol{\varphi}_1(\mathbf{A}t), \qquad \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 + t\,\boldsymbol{\varphi}_1(\mathbf{A}t)\,\mathbf{B}\mathbf{u}_0. \]

When \( \mathbf{A} \) is invertible, \( \boldsymbol{\varphi}_1(\mathbf{A}t)=(e^{\mathbf{A}t}-\mathbf{I})(\mathbf{A}t)^{-1} \) and the formulas coincide. The series definition guarantees correctness for singular \( \mathbf{A} \) without introducing generalized inverses.

Let a step of amplitude \( \mathbf{u}_0 \) occur at time \( t_s \ge 0 \), with the indicator \( \mathbf{1}(t-t_s)=0 \) for \( t < t_s \) and \( 1 \) for \( t \ge t_s \). The forced term becomes:

\[ \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\,\mathbf{u}_0\,\mathbf{1}(\tau-t_s)\, d\tau = \int_{t_s}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\,\mathbf{u}_0\, d\tau, \quad (t \ge t_s). \]

Therefore the solution is naturally piecewise:

\[ \mathbf{x}(t)= \begin{cases} e^{\mathbf{A}t}\mathbf{x}_0, & t < t_s,\\[6pt] e^{\mathbf{A}t}\mathbf{x}_0 + \displaystyle\int_{t_s}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\mathbf{u}_0\, d\tau, & t \ge t_s. \end{cases} \]

Evaluate the integral by shifting time. Let \( \delta=t-t_s \) and set \( s=t-\tau \); when \( \tau=t_s \), we have \( s=\delta \), and when \( \tau=t \), we have \( s=0 \). Thus:

\[ \int_{t_s}^{t} e^{\mathbf{A}(t-\tau)}\, d\tau = \int_{0}^{\delta} e^{\mathbf{A}s}\, ds = \delta\,\boldsymbol{\varphi}_1(\mathbf{A}\delta). \]

So for \( t \ge t_s \),

\[ \mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0 + (t-t_s)\,\boldsymbol{\varphi}_1(\mathbf{A}(t-t_s))\,\mathbf{B}\mathbf{u}_0. \]

Continuity at the switching time. From the integral form, the forced term is zero when the upper and lower limits coincide: at \( t=t_s \), the integral is \( \int_{t_s}^{t_s}(\cdot)=0 \), hence \( \lim_{t ↓ t_s}\mathbf{x}(t)=\lim_{t ↑ t_s}\mathbf{x}(t)=e^{\mathbf{A}t_s}\mathbf{x}_0 \). The state is continuous for bounded inputs.

Let the input be a vector-valued polynomial: \( \mathbf{u}(t)=\sum_{k=0}^{p}\mathbf{u}_k t^k \). By linearity of the integral, it suffices to solve for a monomial input \( \mathbf{u}(t)=\mathbf{u}_k t^k \) and then sum.

For a monomial, the forced term is:

\[ \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\,(\tau^k \mathbf{u}_k)\, d\tau. \]

Perform the normalization \( \tau=t s \), so \( d\tau=t\,ds \) and \( s \in [0,1] \). Then:

\[ \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\,(\tau^k \mathbf{u}_k)\, d\tau = t^{k+1}\left(\int_{0}^{1} e^{\mathbf{A}t(1-s)} s^k\, ds\right)\mathbf{B}\mathbf{u}_k. \]

Define the matrix “phi” functions by the integral representation (valid for all square matrices):

\[ \boldsymbol{\varphi}_{m}(\mathbf{Z}) \triangleq \frac{1}{(m-1)!}\int_{0}^{1} e^{(1-s)\mathbf{Z} }\, s^{m-1}\, ds, \qquad m=1,2,\dots \]

With \( m=k+1 \) we obtain:

\[ \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\,(\tau^k \mathbf{u}_k)\, d\tau = t^{k+1}\,k!\,\boldsymbol{\varphi}_{k+1}(\mathbf{A}t)\,\mathbf{B}\mathbf{u}_k. \]

Summing over \( k \) yields the polynomial-input closed form:

\[ \mathbf{x}(t) = e^{\mathbf{A}t}\mathbf{x}_0 + \sum_{k=0}^{p} t^{k+1}\,k!\,\boldsymbol{\varphi}_{k+1}(\mathbf{A}t)\,\mathbf{B}\mathbf{u}_k. \]

Series form (useful for computation and for singular \( \mathbf{A} \)). Expanding \( e^{(1-s)\mathbf{Z} }=\sum_{j=0}^{\infty}\frac{(1-s)^j\mathbf{Z}^j}{j!} \) and integrating term-by-term gives the standard power series:

\[ \boldsymbol{\varphi}_{m}(\mathbf{Z})=\sum_{j=0}^{\infty}\frac{\mathbf{Z}^j}{(j+m)!}, \qquad m=1,2,\dots \]

This is the exact matrix analogue of scalar special functions used in exponential integrators. In particular, \( \boldsymbol{\varphi}_1(\mathbf{Z})=\sum_{j=0}^{\infty}\frac{\mathbf{Z}^j}{(j+1)!} \) recovers the constant-input formula in Section 2.

Useful recurrence. From the series definition one can show:

\[ \mathbf{Z}\,\boldsymbol{\varphi}_{m+1}(\mathbf{Z}) = \boldsymbol{\varphi}_{m}(\mathbf{Z}) - \frac{1}{m!}\mathbf{I}, \qquad m \ge 1, \]

which is helpful for deriving identities and for numerical evaluation when combined with stable algorithms.

A practical and theoretically clean way to compute integrals like \( \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\,d\tau \) is to embed them into a single matrix exponential. Consider the block matrix \( \mathbf{M}=\begin{bmatrix}\mathbf{A} & \mathbf{B}\\ \mathbf{0} & \mathbf{0}\end{bmatrix} \). Then the exponential has the structure:

\[ e^{\mathbf{M}t} = \begin{bmatrix} e^{\mathbf{A}t} & \displaystyle\int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\, d\tau\\[10pt] \mathbf{0} & \mathbf{I} \end{bmatrix}. \]

Proof sketch (matrix ODE argument). Let \( \mathbf{Z}(t)=e^{\mathbf{M}t} \). Then \( \dot{\mathbf{Z} }(t)=\mathbf{M}\mathbf{Z}(t) \), \( \mathbf{Z}(0)=\mathbf{I} \). Writing blocks \( \mathbf{Z}(t)=\begin{bmatrix}\mathbf{Z}_{11}(t)&\mathbf{Z}_{12}(t)\\ \mathbf{0}&\mathbf{Z}_{22}(t)\end{bmatrix} \) yields: \( \dot{\mathbf{Z} }_{11}=\mathbf{A}\mathbf{Z}_{11} \), \( \dot{\mathbf{Z} }_{12}=\mathbf{A}\mathbf{Z}_{12}+\mathbf{B}\mathbf{Z}_{22} \), and \( \dot{\mathbf{Z} }_{22}=\mathbf{0} \). Hence \( \mathbf{Z}_{11}=e^{\mathbf{A}t} \), \( \mathbf{Z}_{22}=\mathbf{I} \), and variation of constants gives \( \mathbf{Z}_{12}(t)=\int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\, d\tau \).

flowchart TD
  A1["Goal: compute integral K(t) = int_0^t exp(A(t-s)) B ds"] --> B1["Form block matrix M = [[A, B],[0, 0]]"]
  B1 --> C1["Compute E = exp(M t) using a matrix exp routine"]
  C1 --> D1["Extract blocks: E11 = exp(A t), E12 = K(t)"]
  D1 --> E1["Constant input: x(t)=E11 x0 + E12 u0"]
  D1 --> F1["Reuse for step: use t-ts and piecewise logic"]
        

This block approach avoids explicit inversion of \( \mathbf{A} \) and remains valid if \( \mathbf{A} \) is singular. For polynomial inputs, one may either (i) compute \( \boldsymbol{\varphi}_{m}(\mathbf{A}t) \) via series/recurrences, or (ii) use larger block embeddings (conceptually similar to the above) to extract higher-order moments.

The code below computes closed-form solutions using the series definition of \( \boldsymbol{\varphi}_m(\cdot) \), and verifies them using numerical integration. Relevant libraries in Python’s control ecosystem include: numpy, scipy.linalg, scipy.integrate, and (later chapters) control (python-control).

File: Chapter7_Lesson3.py

"""
Chapter7_Lesson3.py
Modern Control — Chapter 7 (Solutions of LTI State Equations), Lesson 3
Solutions for Constant, Step, and Polynomial Inputs
"""

import numpy as np

try:
    from scipy.linalg import expm
    from scipy.integrate import solve_ivp
except ImportError as e:
    raise SystemExit("This script requires SciPy (scipy.linalg.expm, scipy.integrate.solve_ivp).") from e


def phi_series(Z: np.ndarray, m: int, terms: int = 40) -> np.ndarray:
    """
    phi_m(Z) = sum_{j=0}^inf Z^j/(j+m)!
    """
    import math
    n = Z.shape[0]
    I = np.eye(n)
    S = I / math.factorial(m)
    Zpow = I.copy()
    for j in range(1, terms):
        Zpow = Zpow @ Z
        S = S + Zpow / math.factorial(j + m)
    return S


def x_constant_input(A, B, x0, u0, t):
    At = A * t
    E = expm(At)
    Phi1 = phi_series(At, m=1)
    return E @ x0 + (t * (Phi1 @ (B @ u0)))


def x_step_input(A, B, x0, u0, t, ts):
    if t < ts:
        return expm(A * t) @ x0
    dt = t - ts
    E = expm(A * t)
    Phi1 = phi_series(A * dt, m=1)
    return E @ x0 + (dt * (Phi1 @ (B @ u0)))


def x_poly_input(A, B, x0, coeffs, t):
    """
    u(t)=sum_{k=0}^p u_k t^k
    x(t)=exp(A t)x0 + sum t^{k+1} k! phi_{k+1}(A t) B u_k
    """
    import math
    At = A * t
    E = expm(At)
    x = E @ x0
    for k, uk in enumerate(coeffs):
        Phik1 = phi_series(At, m=k + 1)
        x = x + (t ** (k + 1)) * math.factorial(k) * (Phik1 @ (B @ uk))
    return x


def simulate_ivp(A, B, x0, u_fun, t_span, t_eval):
    A = np.asarray(A)
    B = np.asarray(B)
    x0 = np.asarray(x0)

    def f(t, x):
        return (A @ x) + (B @ u_fun(t))

    sol = solve_ivp(f, t_span=t_span, y0=x0, t_eval=t_eval, rtol=1e-10, atol=1e-12)
    return sol.t, sol.y.T


def main():
    import math

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

    u0 = np.array([2.0])
    t_eval = np.linspace(0.0, 5.0, 251)

    # Constant
    x_cf = np.vstack([x_constant_input(A, B, x0, u0, t) for t in t_eval])
    _, x_num = simulate_ivp(A, B, x0, lambda t: u0, (t_eval[0], t_eval[-1]), t_eval)
    print("Constant input max error:", np.max(np.linalg.norm(x_cf - x_num, axis=1)))

    # Step at ts
    ts = 1.5
    x_step_cf = np.vstack([x_step_input(A, B, x0, u0, t, ts) for t in t_eval])
    def u_step(t):
        return u0 if t >= ts else np.zeros_like(u0)
    _, x_num2 = simulate_ivp(A, B, x0, u_step, (t_eval[0], t_eval[-1]), t_eval)
    print("Step input max error:", np.max(np.linalg.norm(x_step_cf - x_num2, axis=1)))

    # Polynomial: u(t)=u0 + u1 t
    u1 = np.array([0.5])
    coeffs = [u0, u1]
    x_poly_cf = np.vstack([x_poly_input(A, B, x0, coeffs, t) for t in t_eval])
    def u_poly(t):
        return u0 + u1 * t
    _, x_num3 = simulate_ivp(A, B, x0, u_poly, (t_eval[0], t_eval[-1]), t_eval)
    print("Polynomial input max error:", np.max(np.linalg.norm(x_poly_cf - x_num3, axis=1)))


if __name__ == "__main__":
    main()
      

Implementation note: the truncated series for \( \boldsymbol{\varphi}_m(\cdot) \) is adequate for small-to-moderate \( \|\mathbf{A}t\| \). For larger problems, robust methods include scaling-and-squaring, Krylov subspace methods, and block-matrix exponential extraction (Section 5).

In C++, a common scientific stack for modern control implementations uses Eigen for linear algebra. The unsupported Eigen module MatrixFunctions provides a matrix exponential routine suitable for demonstrations.

File: Chapter7_Lesson3.cpp

/*
Chapter7_Lesson3.cpp
Dependencies:
  - Eigen
  - unsupported/Eigen/MatrixFunctions (for matrix exponential)
*/

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

using Eigen::MatrixXd;
using Eigen::VectorXd;

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

static MatrixXd phi_series(const MatrixXd& Z, int m, int terms = 40) {
    const int n = static_cast<int>(Z.rows());
    MatrixXd I = MatrixXd::Identity(n, n);
    MatrixXd S = I / factorial(m);
    MatrixXd Zpow = I;
    for (int j = 1; j < terms; ++j) {
        Zpow = Zpow * Z;
        S += Zpow / factorial(j + m);
    }
    return S;
}

static VectorXd x_constant_input(const MatrixXd& A,
                                 const MatrixXd& B,
                                 const VectorXd& x0,
                                 const VectorXd& u0,
                                 double t) {
    MatrixXd At = A * t;
    MatrixXd E = At.exp();
    MatrixXd Phi1 = phi_series(At, 1);
    return (E * x0) + (t * (Phi1 * (B * u0)));
}

static VectorXd x_step_input(const MatrixXd& A,
                             const MatrixXd& B,
                             const VectorXd& x0,
                             const VectorXd& u0,
                             double t,
                             double ts) {
    if (t < ts) return (A * t).exp() * x0;
    MatrixXd E = (A * t).exp();
    double dt = t - ts;
    MatrixXd Phi1 = phi_series(A * dt, 1);
    return (E * x0) + (dt * (Phi1 * (B * u0)));
}

static VectorXd x_poly_input(const MatrixXd& A,
                             const MatrixXd& B,
                             const VectorXd& x0,
                             const std::vector<VectorXd>& coeffs,
                             double t) {
    MatrixXd At = A * t;
    MatrixXd E = At.exp();
    VectorXd x = E * x0;
    for (int k = 0; k < static_cast<int>(coeffs.size()); ++k) {
        MatrixXd Phik1 = phi_series(At, k + 1);
        x += std::pow(t, k + 1) * factorial(k) * (Phik1 * (B * coeffs[k]));
    }
    return x;
}

int main() {
    MatrixXd A(2,2);
    A << 0.0, 1.0,
        -2.0, -3.0;

    MatrixXd B(2,1);
    B << 0.0,
         1.0;

    VectorXd x0(2);
    x0 << 1.0, 0.0;

    VectorXd u0(1); u0 << 2.0;
    VectorXd u1(1); u1 << 0.5;
    std::vector<VectorXd> coeffs = {u0, u1};

    double ts = 1.5;

    for (double t : {0.0, 1.0, 2.0, 5.0}) {
        VectorXd xc = x_constant_input(A, B, x0, u0, t);
        VectorXd xs = x_step_input(A, B, x0, u0, t, ts);
        VectorXd xp = x_poly_input(A, B, x0, coeffs, t);

        std::cout << "t=" << t << "\n";
        std::cout << "  x_const=" << xc.transpose() << "\n";
        std::cout << "  x_step =" << xs.transpose() << "\n";
        std::cout << "  x_poly =" << xp.transpose() << "\n\n";
    }
    return 0;
}
      

Java control stacks often rely on EJML or Apache Commons Math for linear algebra. For pedagogical transparency (and to satisfy “from scratch” requirements), the code below implements small dense-matrix routines and evaluates \( e^{\mathbf{A}t} \) by scaling-and-squaring with a truncated Taylor series, plus \( \boldsymbol{\varphi}_m \) by series.

File: Chapter7_Lesson3.java

/*
Chapter7_Lesson3.java
Self-contained demonstration for small matrices.
*/

import java.util.Arrays;

public class Chapter7_Lesson3 {

    static double[][] eye(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    static double[][] add(double[][] A, double[][] B) {
        int n = A.length, m = A[0].length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                C[i][j] = A[i][j] + B[i][j];
        return C;
    }

    static double[][] scale(double[][] A, double s) {
        int n = A.length, m = A[0].length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++)
            for (int j = 0; j < m; j++)
                C[i][j] = s * A[i][j];
        return C;
    }

    static double[][] mul(double[][] A, double[][] B) {
        int n = A.length, p = A[0].length, m = B[0].length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++) {
            for (int k = 0; k < p; k++) {
                double aik = A[i][k];
                for (int j = 0; j < m; j++) C[i][j] += aik * B[k][j];
            }
        }
        return C;
    }

    static double[] mul(double[][] A, double[] x) {
        int n = A.length, m = A[0].length;
        double[] y = new double[n];
        for (int i = 0; i < n; i++) {
            double s = 0.0;
            for (int j = 0; j < m; j++) s += A[i][j] * x[j];
            y[i] = s;
        }
        return y;
    }

    static double norm1(double[][] A) {
        int n = A.length, m = A[0].length;
        double max = 0.0;
        for (int j = 0; j < m; j++) {
            double col = 0.0;
            for (int i = 0; i < n; i++) col += Math.abs(A[i][j]);
            max = Math.max(max, col);
        }
        return max;
    }

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

    static double[][] expm(double[][] A) {
        int n = A.length;
        double aNorm = norm1(A);
        int s = 0;
        if (aNorm > 0.5) {
            s = (int) Math.ceil(Math.log(aNorm / 0.5) / Math.log(2.0));
        }
        double[][] As = scale(A, 1.0 / Math.pow(2.0, s));

        int K = 30;
        double[][] X = eye(n);
        double[][] term = eye(n);
        for (int k = 1; k <= K; k++) {
            term = mul(term, As);
            X = add(X, scale(term, 1.0 / factorial(k)));
        }
        for (int i = 0; i < s; i++) X = mul(X, X);
        return X;
    }

    static double[][] phiSeries(double[][] Z, int m, int terms) {
        int n = Z.length;
        double[][] S = scale(eye(n), 1.0 / factorial(m));
        double[][] Zpow = eye(n);
        for (int j = 1; j < terms; j++) {
            Zpow = mul(Zpow, Z);
            S = add(S, scale(Zpow, 1.0 / factorial(j + m)));
        }
        return S;
    }

    static double[] xConstant(double[][] A, double[][] B, double[] x0, double[] u0, double t) {
        double[][] At = scale(A, t);
        double[][] E = expm(At);
        double[][] Phi1 = phiSeries(At, 1, 40);

        double[] Bu0 = mul(B, u0);
        double[] forced = mul(scale(Phi1, t), Bu0);

        double[] x = mul(E, x0);
        for (int i = 0; i < x.length; i++) x[i] += forced[i];
        return x;
    }

    static double[] xStep(double[][] A, double[][] B, double[] x0, double[] u0, double t, double ts) {
        if (t < ts) return mul(expm(scale(A, t)), x0);

        double[][] E = expm(scale(A, t));
        double dt = t - ts;
        double[][] Phi1 = phiSeries(scale(A, dt), 1, 40);

        double[] Bu0 = mul(B, u0);
        double[] forced = mul(scale(Phi1, dt), Bu0);

        double[] x = mul(E, x0);
        for (int i = 0; i < x.length; i++) x[i] += forced[i];
        return x;
    }

    static double[] xPoly(double[][] A, double[][] B, double[] x0, double[][] coeffs, double t) {
        double[][] At = scale(A, t);
        double[][] E = expm(At);
        double[] x = mul(E, x0);

        for (int k = 0; k < coeffs.length; k++) {
            double[][] Phik1 = phiSeries(At, k + 1, 45);
            double[] Buk = mul(B, coeffs[k]);
            double scaleFac = Math.pow(t, k + 1) * factorial(k);
            double[] add = mul(scale(Phik1, scaleFac), Buk);
            for (int i = 0; i < x.length; i++) x[i] += add[i];
        }
        return x;
    }

    static String fmt(double[] x) {
        StringBuilder sb = new StringBuilder();
        sb.append("[");
        for (int i = 0; i < x.length; i++) {
            sb.append(String.format("% .6f", x[i]));
            if (i + 1 < x.length) sb.append(", ");
        }
        sb.append("]");
        return sb.toString();
    }

    public static void main(String[] args) {
        double[][] A = new double[][] { {0.0, 1.0}, {-2.0, -3.0} };
        double[][] B = new double[][] { {0.0}, {1.0} };
        double[] x0 = new double[] {1.0, 0.0};

        double[] u0 = new double[] {2.0};
        double[] u1 = new double[] {0.5};
        double[][] coeffs = new double[][] {u0, u1};

        double ts = 1.5;

        for (double t : new double[] {0.0, 1.0, 2.0, 5.0}) {
            System.out.println("t=" + t);
            System.out.println("  x_const=" + fmt(xConstant(A, B, x0, u0, t)));
            System.out.println("  x_step =" + fmt(xStep(A, B, x0, u0, t, ts)));
            System.out.println("  x_poly =" + fmt(xPoly(A, B, x0, coeffs, t)));
            System.out.println();
        }
    }
}
      

MATLAB’s expm provides reliable matrix exponentials, and ode45 can verify results. For Simulink, a State-Space block driven by Step or Ramp matches the inputs studied here.

File: Chapter7_Lesson3.m

% Chapter7_Lesson3.m
% Solutions for Constant, Step, and Polynomial Inputs

clear; clc;

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

u0 = 2;
u1 = 0.5;
ts = 1.5;

tgrid = linspace(0, 5, 251);

phi_series = @(Z,m,terms) local_phi_series(Z,m,terms);

x_const = zeros(2, numel(tgrid));
x_step  = zeros(2, numel(tgrid));
x_poly  = zeros(2, numel(tgrid));

for i = 1:numel(tgrid)
    t = tgrid(i);
    At = A*t;
    E = expm(At);

    Phi1 = phi_series(At, 1, 40);

    % constant input
    x_const(:,i) = E*x0 + t * (Phi1*(B*u0));

    % step input
    if t < ts
        x_step(:,i) = expm(A*t)*x0;
    else
        dt = t - ts;
        Phi1_dt = phi_series(A*dt, 1, 40);
        x_step(:,i) = expm(A*t)*x0 + dt * (Phi1_dt*(B*u0));
    end

    % polynomial input u(t)=u0 + u1 t
    Phi2 = phi_series(At, 2, 45);
    x_poly(:,i) = E*x0 + t * (Phi1*(B*u0)) + (t^2) * (Phi2*(B*u1));
end

opts = odeset('RelTol',1e-10,'AbsTol',1e-12);

% verify constant
ode_const = @(t,x) A*x + B*u0;
[tn1, xn1] = ode45(ode_const, [tgrid(1) tgrid(end)], x0, opts);
xn1i = interp1(tn1, xn1, tgrid);
err1 = max(vecnorm((x_const' - xn1i), 2, 2));
fprintf('Constant input max error: %.3e\n', err1);

% verify step
ode_step = @(t,x) A*x + B*(u0*(t>=ts));
[tn2, xn2] = ode45(ode_step, [tgrid(1) tgrid(end)], x0, opts);
xn2i = interp1(tn2, xn2, tgrid);
err2 = max(vecnorm((x_step' - xn2i), 2, 2));
fprintf('Step input max error: %.3e\n', err2);

% verify polynomial
ode_poly = @(t,x) A*x + B*(u0 + u1*t);
[tn3, xn3] = ode45(ode_poly, [tgrid(1) tgrid(end)], x0, opts);
xn3i = interp1(tn3, xn3, tgrid);
err3 = max(vecnorm((x_poly' - xn3i), 2, 2));
fprintf('Polynomial input max error: %.3e\n', err3);

figure;
plot(tgrid, x_const(1,:), tgrid, x_step(1,:), tgrid, x_poly(1,:));
xlabel('t'); ylabel('x_1(t)'); legend('constant','step','polynomial'); grid on;

function Phi = local_phi_series(Z, m, terms)
    n = size(Z,1);
    Phi = eye(n)/factorial(m);
    Zpow = eye(n);
    for j = 1:terms-1
        Zpow = Zpow*Z;
        Phi = Phi + Zpow/factorial(j+m);
    end
end
      

Mathematica can evaluate the convolution integral symbolically for many low-dimensional systems, providing an independent verification of the analytic expressions.

File: Chapter7_Lesson3.nb

(* Chapter7_Lesson3.nb (Wolfram Language notebook expression)
   Symbolic verification for constant, step, and polynomial inputs. *)

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

u0 = 2;
xConst[t_] := MatrixExp[A t].x0 +
  Integrate[MatrixExp[A (t - \[Tau])].B u0, {\[Tau], 0, t}] // Simplify;

ts = 3/2;
uStep[t_] := u0*UnitStep[t - ts];
xStep[t_] := MatrixExp[A t].x0 +
  Integrate[MatrixExp[A (t - \[Tau])].B uStep[\[Tau]], {\[Tau], 0, t}] //
    PiecewiseExpand // Simplify;

u1 = 1/2;
uPoly[t_] := u0 + u1 t;
xPoly[t_] := MatrixExp[A t].x0 +
  Integrate[MatrixExp[A (t - \[Tau])].B uPoly[\[Tau]], {\[Tau], 0, t}] // Simplify;

{xConst[5], xStep[2], xPoly[5]} // N
      

Problem 1 (Constant input, invertible \( \mathbf{A} \)): Assume \( \mathbf{A} \) is invertible and \( \mathbf{u}(t)=\mathbf{u}_0 \). Starting from \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 + \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\mathbf{u}_0\, d\tau \), derive the formula \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 + \mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I})\mathbf{B}\mathbf{u}_0 \).

Solution: Substitute \( s=t-\tau \) to rewrite the integral as \( \int_{0}^{t}e^{\mathbf{A}s}\,ds \). Define \( \mathbf{F}(t)=\mathbf{A}^{-1}(e^{\mathbf{A}t}-\mathbf{I}) \). Then \( \mathbf{F}'(t)=\mathbf{A}^{-1}(\mathbf{A}e^{\mathbf{A}t})=e^{\mathbf{A}t \!} \) and \( \mathbf{F}(0)=\mathbf{0} \), so \( \int_{0}^{t}e^{\mathbf{A}s}\,ds=\mathbf{F}(t) \). Substituting yields the result.

Problem 2 (Step input time shift): Let \( \mathbf{u}(t)=\mathbf{u}_0\,\mathbf{1}(t-t_s) \). Show that for \( t \ge t_s \), \( \int_{t_s}^{t} e^{\mathbf{A}(t-\tau)}\,d\tau = (t-t_s)\boldsymbol{\varphi}_1(\mathbf{A}(t-t_s)) \).

Solution: Let \( \delta=t-t_s \) and substitute \( s=t-\tau \). The integral becomes \( \int_{0}^{\delta} e^{\mathbf{A}s}\,ds \). By the definition \( \int_{0}^{\delta} e^{\mathbf{A}s}\,ds=\delta\,\boldsymbol{\varphi}_1(\mathbf{A}\delta) \) (series or integral representation), concluding the proof.

Problem 3 (Polynomial input, affine case): For \( \mathbf{u}(t)=\mathbf{u}_0+\mathbf{u}_1 t \), derive: \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 + t\boldsymbol{\varphi}_1(\mathbf{A}t)\mathbf{B}\mathbf{u}_0 + t^2 \boldsymbol{\varphi}_2(\mathbf{A}t)\mathbf{B}\mathbf{u}_1 \).

Solution: Use linearity to split the input into constant and monomial terms. The constant term contributes \( t\boldsymbol{\varphi}_1(\mathbf{A}t)\mathbf{B}\mathbf{u}_0 \). For the monomial term \( \mathbf{u}_1 t \), set \( k=1 \) in \( \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}(\tau^k\mathbf{u}_k)\,d\tau = t^{k+1}k!\boldsymbol{\varphi}_{k+1}(\mathbf{A}t)\mathbf{B}\mathbf{u}_k \), giving \( t^2\boldsymbol{\varphi}_2(\mathbf{A}t)\mathbf{B}\mathbf{u}_1 \).

Problem 4 (Recurrence of \( \boldsymbol{\varphi}_m \)): Starting from the series definition \( \boldsymbol{\varphi}_{m}(\mathbf{Z})=\sum_{j=0}^{\infty}\frac{\mathbf{Z}^j}{(j+m)!} \), show: \( \mathbf{Z}\boldsymbol{\varphi}_{m+1}(\mathbf{Z})=\boldsymbol{\varphi}_{m}(\mathbf{Z})-\frac{1}{m!}\mathbf{I} \).

Solution: Compute: \( \mathbf{Z}\boldsymbol{\varphi}_{m+1}(\mathbf{Z})=\sum_{j=0}^{\infty}\frac{\mathbf{Z}^{j+1} }{(j+m+1)!} =\sum_{j=1}^{\infty}\frac{\mathbf{Z}^{j} }{(j+m)!} \). Meanwhile, \( \boldsymbol{\varphi}_{m}(\mathbf{Z})=\frac{1}{m!}\mathbf{I}+\sum_{j=1}^{\infty}\frac{\mathbf{Z}^{j} }{(j+m)!} \). Subtracting gives the recurrence.

We converted the convolution integral solution of LTI state equations into closed-form expressions for constant, step, and polynomial inputs. The main tool is the family of matrix functions \( \boldsymbol{\varphi}_m(\cdot) \), which yields formulas valid even when \( \mathbf{A} \) is singular. We also presented a block-matrix exponential embedding that extracts forced-response integrals from a single exponential evaluation—an idea that later supports stable numerical simulation methods.