Consider the continuous-time LTI state equation introduced earlier: \( \dot{\mathbf{x} }(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t) \), with initial condition \( \mathbf{x}(0) = \mathbf{x}_0 \). Here \( \mathbf{x}(t)\in\mathbb{R}^n \), \( \mathbf{u}(t)\in\mathbb{R}^m \), \( \mathbf{A}\in\mathbb{R}^{n\times n} \), \( \mathbf{B}\in\mathbb{R}^{n\times m} \).

From Lesson 1, the homogeneous solution of \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x} \) is \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}(0) \). We now derive the full solution with input, and isolate the forced response (also called the zero-state response) that depends only on \( \mathbf{u} \).

We will frequently use the matrix exponential properties already established: \( \frac{d}{dt}e^{\mathbf{A}t}=\mathbf{A}e^{\mathbf{A}t}=e^{\mathbf{A}t}\mathbf{A} \) and \( e^{\mathbf{A}0}=\mathbf{I} \).

The central idea is to seek a solution of the form \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{c}(t) \), where \( \mathbf{c}(t)\in\mathbb{R}^n \) is an unknown time-varying vector.

2.1 Theorem (General solution of LTI state equation)

Assume \( \mathbf{u}(\cdot) \) is piecewise continuous on \( [0,T] \). Then the unique solution on \( [0,T] \) is

\[ \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0+\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds,\quad 0\le t\le T. \]

2.2 Proof

Let \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{c}(t) \). Differentiate:

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

Substitute into the state equation \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \):

\[ \mathbf{A}e^{\mathbf{A}t}\mathbf{c}(t)+e^{\mathbf{A}t}\dot{\mathbf{c} }(t) =\mathbf{A}\left(e^{\mathbf{A}t}\mathbf{c}(t)\right)+\mathbf{B}\mathbf{u}(t). \]

Cancel the common term \( \mathbf{A}e^{\mathbf{A}t}\mathbf{c}(t) \) to obtain:

\[ e^{\mathbf{A}t}\dot{\mathbf{c} }(t)=\mathbf{B}\mathbf{u}(t) \quad\Rightarrow\quad \dot{\mathbf{c} }(t)=e^{-\mathbf{A}t}\mathbf{B}\mathbf{u}(t). \]

Integrate from \( 0 \) to \( t \):

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

Since \( \mathbf{x}(0)=e^{\mathbf{A}0}\mathbf{c}(0)=\mathbf{c}(0)=\mathbf{x}_0 \), we have \( \mathbf{c}(0)=\mathbf{x}_0 \), hence:

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

Multiply by \( e^{\mathbf{A}t} \):

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

This is the desired formula. Uniqueness follows from standard existence/uniqueness for linear ODEs (already discussed in Chapter 3), since the difference of any two solutions satisfies the homogeneous equation.

Decompose the state as the sum of zero-input and zero-state parts:

\[ \mathbf{x}(t)=\underbrace{e^{\mathbf{A}t}\mathbf{x}_0}_{\mathbf{x}_{zi}(t)} +\underbrace{\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds}_{\mathbf{x}_{zs}(t)}. \]

The forced response is the mapping \( \mathbf{u}(\cdot)\mapsto\mathbf{x}_{zs}(\cdot) \). Define the matrix-valued kernel (state impulse response kernel):

\[ \mathbf{H}(t)\triangleq e^{\mathbf{A}t}\mathbf{B},\quad t\ge 0, \qquad \mathbf{H}(t)=\mathbf{0},\quad t<0. \]

Then the forced response can be written in the standard convolution form (matrix convolution):

\[ \mathbf{x}_{zs}(t)=\int_{-\infty}^{\infty} \mathbf{H}(t-s)\mathbf{u}(s)\,ds =\int_{0}^{t} \mathbf{H}(t-s)\mathbf{u}(s)\,ds, \quad t\ge 0. \]

The second equality uses causality: \( \mathbf{H}(t-s)=\mathbf{0} \) when \( s>t \).

flowchart TD
  U["input u(t)"] --> K["kernel H(t) = exp(A t) B"]
  K --> I["integrate: x_zs(t) = integral_0^t H(t-s) u(s) ds"]
  X["initial state x0"] --> ZI["zero-input: x_zi(t) = exp(A t) x0"]
  I --> XS["forced part x_zs(t)"]
  ZI --> SUM["sum"]
  XS --> SUM
  SUM --> XT["state x(t)"]
        

Define the forced-response operator \( \mathcal{T} \) by \( (\mathcal{T}\mathbf{u})(t)=\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds \). We now show the core system-theoretic properties.

4.1 Linearity

For scalars \( \alpha,\beta\in\mathbb{R} \) and inputs \( \mathbf{u}_1,\mathbf{u}_2 \):

\[ (\mathcal{T}(\alpha\mathbf{u}_1+\beta\mathbf{u}_2))(t) =\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}(\alpha\mathbf{u}_1(s)+\beta\mathbf{u}_2(s))\,ds =\alpha(\mathcal{T}\mathbf{u}_1)(t)+\beta(\mathcal{T}\mathbf{u}_2)(t). \]

4.2 Causality

If two inputs satisfy \( \mathbf{u}_1(s)=\mathbf{u}_2(s) \) for all \( s\in[0,t] \), then their forced responses at time \( t \) are equal, because the integral uses only values on \( [0,t] \).

4.3 Time invariance (with zero initial state)

Consider the forced response with \( \mathbf{x}_0=\mathbf{0} \). Let the time-shifted input be \( \mathbf{u}_\Delta(t)=\mathbf{u}(t-\Delta) \) with \( \Delta\ge 0 \), and define \( \mathbf{u}(t)=\mathbf{0} \) for \( t<0 \). Then:

\[ (\mathcal{T}\mathbf{u}_\Delta)(t) =\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s-\Delta)\,ds =\int_{-\Delta}^{t-\Delta} e^{\mathbf{A}(t-(r+\Delta))}\mathbf{B}\mathbf{u}(r)\,dr. \]

Since \( \mathbf{u}(r)=\mathbf{0} \) for \( r<0 \), the lower limit can be replaced by \( 0 \):

\[ (\mathcal{T}\mathbf{u}_\Delta)(t)=\int_{0}^{t-\Delta} e^{\mathbf{A}((t-\Delta)-r)}\mathbf{B}\mathbf{u}(r)\,dr =(\mathcal{T}\mathbf{u})(t-\Delta). \]

Thus, the forced-response mapping is time-invariant (in the standard causal sense) when the initial state is zero.

It is often useful to verify directly that \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0+\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds \) satisfies the differential equation. Let \( \mathbf{x}_{zs}(t)=\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds \).

Using Leibniz’s rule for differentiating an integral with variable upper limit:

\[ \frac{d}{dt}\mathbf{x}_{zs}(t) =e^{\mathbf{A}(t-t)}\mathbf{B}\mathbf{u}(t)+\int_{0}^{t}\frac{\partial}{\partial t}\left(e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\right)\,ds. \]

Since \( e^{\mathbf{A}0}=\mathbf{I} \) and \( \frac{\partial}{\partial t}e^{\mathbf{A}(t-s)}=\mathbf{A}e^{\mathbf{A}(t-s)} \), we get:

\[ \dot{\mathbf{x} }_{zs}(t)=\mathbf{B}\mathbf{u}(t)+\int_{0}^{t}\mathbf{A}e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds =\mathbf{B}\mathbf{u}(t)+\mathbf{A}\mathbf{x}_{zs}(t). \]

Similarly, \( \frac{d}{dt}(e^{\mathbf{A}t}\mathbf{x}_0)=\mathbf{A}e^{\mathbf{A}t}\mathbf{x}_0 \). Summing the derivatives yields:

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t), \qquad \mathbf{x}(0)=\mathbf{x}_0. \]

This identity is also the starting point for several numerical schemes (quadrature, time-stepping), which we implement next.

For most nontrivial inputs \( \mathbf{u}(t) \), the forced response integral is not closed-form and must be computed. A standard approach is to discretize time \( t_k = k\Delta t \) and approximate:

\[ \mathbf{x}(t_k)\approx e^{\mathbf{A}t_k}\mathbf{x}_0+\sum_{j=0}^{k-1} e^{\mathbf{A}(t_k-t_j)}\mathbf{B}\mathbf{u}(t_j)\,\Delta t \quad \text{(left-rectangle quadrature)}. \]

Higher accuracy can be obtained with trapezoidal or Simpson rules, but the conceptual structure remains the same: a weighted sum of \( e^{\mathbf{A}(t_k-t_j)}\mathbf{B} \) applied to sampled inputs.

flowchart TD
  S["given A, B, x0, input u(t), step dt, horizon T"] --> G["build time grid t_k = k dt"]
  G --> P["precompute Phi_k = exp(A * (k*dt)) for needed k"]
  P --> Q["for each k: sum_j exp(A * (t_k - t_j)) * B * u(t_j) * w_j"]
  Q --> X["x(t_k) = exp(A*t_k) * x0 + sum"]
  X --> V["validate: check xdot approx A x + B u (optional)"]
        

Python offers two complementary workflows: (i) compute the convolution integral directly using scipy.linalg.expm and numerical quadrature; and (ii) validate against a state-space simulator (python-control or scipy.signal.StateSpace).

7.1 Direct convolution integral (quadrature)

import numpy as np
from scipy.linalg import expm

def forced_response_convolution(A, B, u_fun, t_grid):
    """
    Compute x_zs(t_k) = ∫_0^{t_k} exp(A (t_k - s)) B u(s) ds
    using trapezoidal rule on the provided time grid.
    """
    n = A.shape[0]
    xzs = np.zeros((len(t_grid), n))
    for k, tk in enumerate(t_grid):
        if k == 0:
            xzs[k] = 0.0
            continue
        # trapezoidal weights
        s = t_grid[:k+1]
        us = np.array([u_fun(si) for si in s])  # shape (k+1, m)
        # accumulate integral
        acc = np.zeros((n,))
        for j in range(k):
            sj, sj1 = s[j], s[j+1]
            dt = sj1 - sj
            K0 = expm(A * (tk - sj))  @ B @ us[j]
            K1 = expm(A * (tk - sj1)) @ B @ us[j+1]
            acc += 0.5 * (K0 + K1) * dt
        xzs[k] = acc
    return xzs

def total_state(A, x0, xzs, t_grid):
    X = np.zeros_like(xzs)
    for k, tk in enumerate(t_grid):
        X[k] = expm(A * tk) @ x0 + xzs[k]
    return X

# Example (2-state, 1-input)
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])

# Define a piecewise-continuous input u(t)
def u_fun(t):
    # u(t) = 1 for t in [0,2], else 0 (unit pulse)
    return np.array([1.0 if 0.0 <= t <= 2.0 else 0.0])

t_grid = np.linspace(0.0, 6.0, 601)  # dt = 0.01
xzs = forced_response_convolution(A, B, u_fun, t_grid)
x = total_state(A, x0, xzs, t_grid)

print("x(0) =", x[0])
print("x(6) =", x[-1])

7.2 Cross-check with a state-space simulator

import numpy as np
from scipy.signal import StateSpace, lsim

# Build an equivalent state-space model for simulation check
C = np.eye(2)
D = np.zeros((2,1))
sys = StateSpace(A, B, C, D)

# Create input samples on t_grid
u_samples = np.array([u_fun(t)[0] for t in t_grid])

# lsim returns state trajectory if requested (SciPy version dependent)
tout, yout, xout = lsim(sys, U=u_samples, T=t_grid, X0=x0)

# Compare xout (simulation) to x (convolution)
err = np.max(np.linalg.norm(xout - x, axis=1))
print("max ||x_sim - x_conv|| =", err)

Recommended libraries for modern control workflows in Python: \( \)

• numpy, scipy (matrix exponential, numerical integration, simulation)
• control (python-control) for state-space models, forced response utilities, and LTI analysis
• sympy (optional) for symbolic checks in small examples

In C++, a standard approach uses Eigen for linear algebra and its matrix exponential (unsupported module). The implementation below computes the trapezoidal approximation of the convolution integral.

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

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

using InputFun = std::function<VectorXd(double)>;

std::vector<VectorXd> forcedResponseConvolution(
    const MatrixXd& A,
    const MatrixXd& B,
    const InputFun& u_fun,
    const std::vector<double>& t_grid)
{
    const int n = A.rows();
    std::vector<VectorXd> xzs(t_grid.size(), VectorXd::Zero(n));

    for (size_t k = 1; k < t_grid.size(); ++k)
    {
        double tk = t_grid[k];
        VectorXd acc = VectorXd::Zero(n);

        for (size_t j = 0; j < k; ++j)
        {
            double sj  = t_grid[j];
            double sj1 = t_grid[j+1];
            double dt  = sj1 - sj;

            VectorXd uj  = u_fun(sj);
            VectorXd uj1 = u_fun(sj1);

            MatrixXd K0 = (A * (tk - sj)).exp()  * B;
            MatrixXd K1 = (A * (tk - sj1)).exp() * B;

            acc += 0.5 * (K0 * uj + K1 * uj1) * dt;
        }
        xzs[k] = acc;
    }
    return xzs;
}

std::vector<VectorXd> totalState(
    const MatrixXd& A,
    const VectorXd& x0,
    const std::vector<VectorXd>& xzs,
    const std::vector<double>& t_grid)
{
    std::vector<VectorXd> x(t_grid.size(), x0);
    for (size_t k = 0; k < t_grid.size(); ++k)
    {
        double tk = t_grid[k];
        x[k] = (A * tk).exp() * x0 + xzs[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;

    auto u_fun = [](double t) -> VectorXd {
        VectorXd u(1);
        u(0) = (t >= 0.0 && t <= 2.0) ? 1.0 : 0.0; // unit pulse
        return u;
    };

    std::vector<double> t_grid;
    double T = 6.0, dt = 0.01;
    for (int k = 0; k <= static_cast<int>(T/dt); ++k)
        t_grid.push_back(k * dt);

    auto xzs = forcedResponseConvolution(A, B, u_fun, t_grid);
    auto x   = totalState(A, x0, xzs, t_grid);

    std::cout << "x(0) = " << x.front().transpose() << std::endl;
    std::cout << "x(6) = " << x.back().transpose() << std::endl;
    return 0;
}

Common C++ libraries used in modern control computations include Eigen (linear algebra), Armadillo, and Boost (numerics). For high reliability, pay attention to matrix exponential accuracy and time-step selection.

Java workflows often use EJML (Efficient Java Matrix Library). Since a robust matrix exponential is not always included by default, a practical approach is (a) implement a scaling-and-squaring Padé exponential, or (b) use a numerical ODE solver to validate. Below is a compact educational implementation using a truncated series for \( e^{\mathbf{A}t} \) (suitable only for small \( \|\mathbf{A}t\| \)).

For production-grade work, replace expmSeries with a Padé/scaling-and-squaring exponential.

import java.util.function.DoubleFunction;

public class ConvolutionForcedResponse {

    // Minimal matrix utilities (educational only)
    static double[][] matMul(double[][] A, double[][] B) {
        int n = A.length, m = B[0].length, p = B.length;
        double[][] C = new double[n][m];
        for (int i = 0; i < n; i++)
            for (int k = 0; k < p; k++)
                for (int j = 0; j < m; j++)
                    C[i][j] += A[i][k] * B[k][j];
        return C;
    }

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

    static double[][] matAdd(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[][] matScale(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[][] eye(int n) {
        double[][] I = new double[n][n];
        for (int i = 0; i < n; i++) I[i][i] = 1.0;
        return I;
    }

    // Truncated series expm(A*t) ≈ Σ_{k=0}^{K} (A*t)^k / k!
    // Educational only; replace with Padé + scaling/squaring for robustness.
    static double[][] expmSeries(double[][] A, double t, int K) {
        int n = A.length;
        double[][] At = matScale(A, t);
        double[][] E = eye(n);
        double[][] term = eye(n);

        for (int k = 1; k <= K; k++) {
            term = matMul(term, At);
            double invFact = 1.0;
            for (int j = 2; j <= k; j++) invFact /= j;
            E = matAdd(E, matScale(term, invFact));
        }
        return E;
    }

    static double[] addVec(double[] a, double[] b) {
        double[] c = new double[a.length];
        for (int i = 0; i < a.length; i++) c[i] = a[i] + b[i];
        return c;
    }

    static double[] scaleVec(double[] a, double s) {
        double[] b = new double[a.length];
        for (int i = 0; i < a.length; i++) b[i] = s * a[i];
        return b;
    }

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

        DoubleFunction<double[]> uFun = (double tt) -> {
            return new double[] { (tt >= 0.0 && tt <= 2.0) ? 1.0 : 0.0 };
        };

        double T = 6.0, dt = 0.01;
        int N = (int)Math.round(T/dt);

        double[][] x = new double[N+1][2];

        // Precompute expm(A * (k*dt)) as needed (simple, but expensive if done naively)
        for (int k = 0; k <= N; k++) {
            double tk = k * dt;

            // forced part via trapezoidal sum
            double[] acc = new double[] {0.0, 0.0};
            for (int j = 0; j < k; j++) {
                double sj  = j * dt;
                double sj1 = (j+1) * dt;

                double[] uj  = uFun.apply(sj);
                double[] uj1 = uFun.apply(sj1);

                double[][] K0 = matMul(expmSeries(A, tk - sj, 20), B);
                double[][] K1 = matMul(expmSeries(A, tk - sj1, 20), B);

                double[] v0 = matVec(K0, uj);
                double[] v1 = matVec(K1, uj1);
                acc = addVec(acc, scaleVec(addVec(v0, v1), 0.5 * dt));
            }

            double[] zi = matVec(expmSeries(A, tk, 20), x0);
            x[k] = addVec(zi, acc);
        }

        System.out.println("x(0) = [" + x[0][0] + ", " + x[0][1] + "]");
        System.out.println("x(6) = [" + x[N][0] + ", " + x[N][1] + "]");
    }
}

Practical note: the truncated series is included to keep the lesson self-contained. In real modern-control code, use a robust matrix exponential (Padé + scaling/squaring) or integrate the ODE directly.

10.1 MATLAB: direct convolution integral

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

% Input u(t): unit pulse on [0,2]
u = @(t) (t >= 0 && t <= 2);

% Forced response x_zs(t) = ∫_0^t exp(A*(t-s))*B*u(s) ds
x_zs = @(t) integral(@(s) expm(A*(t - s)) * B * u(s), 0, t, 'ArrayValued', true);

% Total state x(t) = expm(A*t)*x0 + x_zs(t)
t = linspace(0, 6, 601);
X = zeros(2, numel(t));
for k = 1:numel(t)
    tk = t(k);
    X(:,k) = expm(A*tk)*x0 + x_zs(tk);
end

disp(X(:,1));
disp(X(:,end));

10.2 MATLAB: validation using built-in simulation

A = [0 1; -2 -3];
B = [0; 1];
C = eye(2);
D = zeros(2,1);
sys = ss(A,B,C,D);

t = linspace(0, 6, 601);
u = double(t >= 0 & t <= 2);  % sampled pulse input
x0 = [1;0];

[y, t_out, x_out] = lsim(sys, u, t, x0);

% Compare x_out' (2xN) with your convolution result X
% (Assuming X computed as above)
err = max(vecnorm((x_out' - X), 2, 1));
disp(err);

10.3 Simulink workflow (conceptual)

In Simulink, implement the forced response without explicitly coding the integral by using the continuous State-Space block (which internally solves \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{u} \)). For educational alignment with this lesson, you can also build a “convolution evaluator” using:

• A MATLAB Function block that computes \( e^{\mathbf{A}(t-s)} \) and a numerical quadrature over a buffer of past inputs.
• A Transport Delay + Integrator structure to approximate the convolution (less precise, but conceptually illustrative).

Mathematica can compute both symbolic and numerical variants of the forced response. For most multi-state systems, the numerical approach is robust and concise.

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

u[t_] := If[0 <= t <= 2, {1}, {0}]; (* 1-input as a 1-vector *)

xZS[t_?NumericQ] := NIntegrate[
  MatrixExp[A (t - s)].B.u[s],
  {s, 0, t}
] // Flatten;

xTotal[t_?NumericQ] := MatrixExp[A t].x0 + xZS[t];

xTotal[0]
xTotal[6]

For symbolic experiments (small dimension, simple inputs), you may try Integrate in place of NIntegrate, but expect expressions to grow quickly.

The following problems are designed to reinforce the convolution integral structure and its role as the forced response in LTI state equations. Unless otherwise stated, assume \( \mathbf{u}(\cdot) \) is piecewise continuous.

Problem 1 (Verification by Differentiation): Let \( \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0+\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds \). Prove that \( \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t) \) and \( \mathbf{x}(0)=\mathbf{x}_0 \).

Solution: Evaluate \( \mathbf{x}(0) \):

\[ \mathbf{x}(0)=e^{\mathbf{A}0}\mathbf{x}_0+\int_{0}^{0}(\cdots)\,ds=\mathbf{I}\mathbf{x}_0+\mathbf{0}=\mathbf{x}_0. \]

Differentiate using Leibniz’s rule exactly as in Section 5 to obtain \( \dot{\mathbf{x} }_{zs}(t)=\mathbf{A}\mathbf{x}_{zs}(t)+\mathbf{B}\mathbf{u}(t) \) and \( \frac{d}{dt}(e^{\mathbf{A}t}\mathbf{x}_0)=\mathbf{A}e^{\mathbf{A}t}\mathbf{x}_0 \). Summing yields the result.

Problem 2 (Scalar Convolution Form): Consider the scalar system \( \dot{x}(t)=a x(t)+b u(t) \), \( x(0)=x_0 \). Derive the forced response explicitly and identify the impulse-response kernel \( h(t) \).

Solution: Here \( e^{At} \) reduces to \( e^{at} \):

\[ x(t)=e^{at}x_0+\int_{0}^{t} e^{a(t-s)}\,b\,u(s)\,ds. \]

With \( x_0=0 \), define \( h(t)=b e^{at} \) for \( t\ge 0 \) and \( h(t)=0 \) for \( t<0 \), then

\[ x_{zs}(t)=\int_{0}^{t} h(t-s)u(s)\,ds=(h*u)(t), \]

which matches the classical convolution form of a causal LTI system.

Problem 3 (Diagonalizable A: Kernel Representation): Assume \( \mathbf{A}=\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1} \) with \( \boldsymbol{\Lambda}=\mathrm{diag}(\lambda_1,\dots,\lambda_n) \). Show that the forced response can be written as \( \mathbf{x}_{zs}(t)=\mathbf{V}\int_{0}^{t} e^{\boldsymbol{\Lambda}(t-s)}\mathbf{V}^{-1}\mathbf{B}\mathbf{u}(s)\,ds \).

Solution: Use the similarity invariance of the exponential:

\[ e^{\mathbf{A}t}=e^{\mathbf{V}\boldsymbol{\Lambda}\mathbf{V}^{-1}t} =\mathbf{V}e^{\boldsymbol{\Lambda}t}\mathbf{V}^{-1}. \]

Substitute into the forced response:

\[ \mathbf{x}_{zs}(t)=\int_{0}^{t}\mathbf{V}e^{\boldsymbol{\Lambda}(t-s)}\mathbf{V}^{-1}\mathbf{B}\mathbf{u}(s)\,ds =\mathbf{V}\int_{0}^{t} e^{\boldsymbol{\Lambda}(t-s)}\mathbf{V}^{-1}\mathbf{B}\mathbf{u}(s)\,ds. \]

This representation makes clear that the kernel is a superposition of modal exponentials \( e^{\lambda_i(t-s)} \) (basis-dependent through \( \mathbf{V} \)).

Problem 4 (Pulse Input at a Final Time): Let \( \mathbf{x}_0=\mathbf{0} \) and suppose \( \mathbf{u}(t)=\mathbf{u}_0 \) for \( 0\le t\le T \), and \( \mathbf{u}(t)=\mathbf{0} \) for \( t>T \). Derive a compact expression for \( \mathbf{x}(T) \).

Solution: For \( t=T \):

\[ \mathbf{x}(T)=\int_{0}^{T} e^{\mathbf{A}(T-s)}\mathbf{B}\mathbf{u}_0\,ds =\left(\int_{0}^{T} e^{\mathbf{A}(T-s)}\,ds\right)\mathbf{B}\mathbf{u}_0. \]

Substitute \( r=T-s \), so \( s=0 \Rightarrow r=T \) and \( s=T \Rightarrow r=0 \):

\[ \mathbf{x}(T)=\left(\int_{T}^{0} e^{\mathbf{A}r}(-dr)\right)\mathbf{B}\mathbf{u}_0 =\left(\int_{0}^{T} e^{\mathbf{A}r}\,dr\right)\mathbf{B}\mathbf{u}_0. \]

If \( \mathbf{A} \) is invertible, one may further simplify via \( \int_{0}^{T} e^{\mathbf{A}r}\,dr=\mathbf{A}^{-1}(e^{\mathbf{A}T}-\mathbf{I}) \):

\[ \mathbf{x}(T)=\mathbf{A}^{-1}\left(e^{\mathbf{A}T}-\mathbf{I}\right)\mathbf{B}\mathbf{u}_0. \]

Problem 5 (Causality of the Kernel): Define \( \mathbf{H}(t)=e^{\mathbf{A}t}\mathbf{B} \) for \( t\ge 0 \) and \( \mathbf{H}(t)=\mathbf{0} \) for \( t<0 \). Prove that the convolution \( \mathbf{x}_{zs}(t)=\int_{-\infty}^{\infty}\mathbf{H}(t-s)\mathbf{u}(s)\,ds \) reduces to \( \int_{0}^{t} \mathbf{H}(t-s)\mathbf{u}(s)\,ds \) for \( t\ge 0 \).

Solution: For \( t\ge 0 \):

\[ \int_{-\infty}^{\infty}\mathbf{H}(t-s)\mathbf{u}(s)\,ds =\int_{-\infty}^{t}\mathbf{H}(t-s)\mathbf{u}(s)\,ds+\int_{t}^{\infty}\mathbf{H}(t-s)\mathbf{u}(s)\,ds. \]

In the second integral, \( s>t \Rightarrow t-s<0 \Rightarrow \mathbf{H}(t-s)=\mathbf{0} \), so it vanishes. If additionally \( \mathbf{u}(s)=\mathbf{0} \) for \( s<0 \) (causal inputs), the first integral becomes \( \int_{0}^{t}\mathbf{H}(t-s)\mathbf{u}(s)\,ds \).

We derived the full solution of the LTI state equation using variation of constants and identified the forced response as a convolution-type integral: \( \mathbf{x}_{zs}(t)=\int_{0}^{t} e^{\mathbf{A}(t-s)}\mathbf{B}\mathbf{u}(s)\,ds \). This formulation establishes a rigorous input-to-state map, clarifies causality and time invariance (for zero initial state), and provides a computational blueprint for simulation via quadrature. In the next lesson, we specialize this integral to common input classes (constant, step, polynomial) to obtain closed-form expressions and efficient evaluation strategies.