In Chapters 6–11, students studied deterministic linear systems through transfer functions and state-space models. In Chapter 17 (Lessons 1–2), random processes, stationarity, autocorrelation, and power spectral density (PSD) were introduced. This lesson combines these ideas: when the input is random, the system output is also random, but its mean, covariance, and variance can often be propagated using deterministic equations.

The core engineering questions are:

• How does the \( \mu_u(t) \) (input mean) affect the \( \mu_y(t) \) (output mean)?
• How does input uncertainty (variance/covariance or PSD) affect output uncertainty?
• Under what conditions does the variance converge to a finite steady-state value?

flowchart TD
  A["Linear model (ODE / state-space / transfer function)"] --> B["Split input into mean + zero-mean fluctuation"]
  B --> C["Propagate mean with \ndeterministic dynamics"]
  B --> D["Propagate covariance / variance"]
  D --> E["Transient uncertainty vs steady-state uncertainty"]
  E --> F["Time-domain view (covariance)"]
  E --> G["Frequency-domain view (PSD)"]
  F --> H["Validate using Monte Carlo simulation"]
  G --> H
        

Consider a causal continuous-time LTI system with impulse response \( g(t) \) and random input \( u(t) \). The output is

\[ y(t) = \int_{0}^{t} g(t-\tau)\,u(\tau)\,d\tau \]

when initial conditions are zero (or when the homogeneous part has been treated separately). Write the input as

\[ u(t) = \mu_u(t) + \tilde{u}(t), \qquad E[\tilde{u}(t)] = 0 \]

By linearity of expectation, the output mean is

\[ \mu_y(t) = E[y(t)] = \int_{0}^{t} g(t-\tau)\,\mu_u(\tau)\,d\tau \]

so the mean propagates through the same LTI operator as a deterministic input.

For second-order statistics, define the zero-mean output fluctuation \( \tilde{y}(t)=y(t)-\mu_y(t) \). Then

\[ \tilde{y}(t)=\int_{0}^{t} g(t-\tau)\,\tilde{u}(\tau)\,d\tau \]

and the (possibly non-stationary) covariance is

\[ R_{\tilde{y} }(t_1,t_2)=E[\tilde{y}(t_1)\tilde{y}(t_2)] = \int_{0}^{t_1}\int_{0}^{t_2} g(t_1-\tau_1)g(t_2-\tau_2) R_{\tilde{u} }(\tau_1,\tau_2)\,d\tau_1 d\tau_2 \]

This formula is exact and is the time-domain basis for mean/variance propagation.

Proof sketch (mean propagation):

Since integration and expectation can be interchanged for finite-energy kernels under standard regularity assumptions,

\[ E[y(t)] = E\!\left[\int_0^t g(t-\tau)u(\tau)\,d\tau\right] = \int_0^t g(t-\tau)E[u(\tau)]\,d\tau = \int_0^t g(t-\tau)\mu_u(\tau)\,d\tau \]

For control engineering, a state-space description is usually more convenient than direct convolution. Let

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t), \qquad \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t) \]

Decompose the random input as \( \mathbf{u}(t)=\boldsymbol{\mu}_u(t)+\tilde{\mathbf{u} }(t) \) with \( E[\tilde{\mathbf{u} }(t)]=\mathbf{0} \). Taking expectation gives

\[ \dot{\boldsymbol{\mu} }_x(t)=\mathbf{A}\boldsymbol{\mu}_x(t)+\mathbf{B}\boldsymbol{\mu}_u(t), \qquad \boldsymbol{\mu}_y(t)=\mathbf{C}\boldsymbol{\mu}_x(t)+\mathbf{D}\boldsymbol{\mu}_u(t) \]

where \( \boldsymbol{\mu}_x(t)=E[\mathbf{x}(t)] \). Thus the mean state evolves exactly like a deterministic state driven by the input mean.

Closed-form mean solution:

\[ \boldsymbol{\mu}_x(t)=e^{\mathbf{A}t}\boldsymbol{\mu}_x(0) + \int_{0}^{t} e^{\mathbf{A}(t-\tau)}\mathbf{B}\boldsymbol{\mu}_u(\tau)\,d\tau \]

If \( \boldsymbol{\mu}_u(t)=\bar{\mathbf{u} } \) is constant and \( \mathbf{A} \) is nonsingular, then

\[ \boldsymbol{\mu}_x(t)=e^{\mathbf{A}t}\boldsymbol{\mu}_x(0) + \mathbf{A}^{-1}\!\left(e^{\mathbf{A}t}-\mathbf{I}\right)\mathbf{B}\bar{\mathbf{u} } \]

and if the system is asymptotically stable (all eigenvalues of \( \mathbf{A} \) have negative real parts), the steady-state mean is

\[ \boldsymbol{\mu}_{x,\infty} = -\mathbf{A}^{-1}\mathbf{B}\bar{\mathbf{u} } \]

A very important model in random vibration and control is a linear system driven by an additive zero-mean white process:

\[ \dot{\mathbf{x} }(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\boldsymbol{\mu}_u(t)+\mathbf{\Gamma}\mathbf{w}(t) \]

\[ E[\mathbf{w}(t)] = \mathbf{0}, \qquad E[\mathbf{w}(t)\mathbf{w}(\tau)^{T}] = \mathbf{Q}\,\delta(t-\tau) \]

Here \( \mathbf{\Gamma} \) maps the white input into the state equation and \( \mathbf{Q} \) is the input intensity matrix.

Define the mean-removed state \( \tilde{\mathbf{x} }(t)=\mathbf{x}(t)-\boldsymbol{\mu}_x(t) \), so that

\[ \dot{\tilde{\mathbf{x} } }(t)=\mathbf{A}\tilde{\mathbf{x} }(t)+\mathbf{\Gamma}\mathbf{w}(t) \]

and define the covariance matrix \( \mathbf{P}(t)=E[\tilde{\mathbf{x} }(t)\tilde{\mathbf{x} }(t)^T] \). Then the covariance satisfies the matrix differential equation

\[ \dot{\mathbf{P} }(t)=\mathbf{A}\mathbf{P}(t)+\mathbf{P}(t)\mathbf{A}^{T}+\mathbf{\Gamma}\mathbf{Q}\mathbf{\Gamma}^{T} \]

This is the continuous-time covariance propagation equation (also called the differential Lyapunov equation).

Derivation (product rule):

\[ \frac{d}{dt}\mathbf{P}(t) = \frac{d}{dt}E[\tilde{\mathbf{x} }\tilde{\mathbf{x} }^T] = E[\dot{\tilde{\mathbf{x} } }\tilde{\mathbf{x} }^T] + E[\tilde{\mathbf{x} }\dot{\tilde{\mathbf{x} } }^{T}] \]

\[ = E[(\mathbf{A}\tilde{\mathbf{x} }+\mathbf{\Gamma}\mathbf{w})\tilde{\mathbf{x} }^T] + E[\tilde{\mathbf{x} }(\mathbf{A}\tilde{\mathbf{x} }+\mathbf{\Gamma}\mathbf{w})^T] \]

\[ = \mathbf{A}E[\tilde{\mathbf{x} }\tilde{\mathbf{x} }^T] + E[\tilde{\mathbf{x} }\tilde{\mathbf{x} }^T]\mathbf{A}^{T} + E[\mathbf{\Gamma}\mathbf{w}\tilde{\mathbf{x} }^T] + E[\tilde{\mathbf{x} }\mathbf{w}^{T}\mathbf{\Gamma}^{T}] \]

For ideal white driving terms, the quadratic variation contributes \( \mathbf{\Gamma}\mathbf{Q}\mathbf{\Gamma}^{T} \), yielding the stated result. (In a rigorous stochastic-calculus derivation, this term appears from the second-order differential product.)

The output covariance follows from \( \mathbf{y}=\mathbf{C}\mathbf{x}+\mathbf{D}\mathbf{u} \). If the direct term is deterministic or absent (common in plant outputs), then

\[ \mathbf{P}_y(t)=E[(\mathbf{y}-\boldsymbol{\mu}_y)(\mathbf{y}-\boldsymbol{\mu}_y)^T] = \mathbf{C}\mathbf{P}(t)\mathbf{C}^{T} \]

If \( \mathbf{A} \) is asymptotically stable and the white input intensity \( \mathbf{Q} \) is constant, then the covariance converges (under standard conditions) to a constant matrix \( \mathbf{P}_{\infty} \) satisfying

\[ \mathbf{A}\mathbf{P}_{\infty}+\mathbf{P}_{\infty}\mathbf{A}^{T} + \mathbf{\Gamma}\mathbf{Q}\mathbf{\Gamma}^{T}= \mathbf{0} \]

This is the continuous-time algebraic Lyapunov equation.

The finite steady-state variance condition is directly tied to stability: if any mode of \( \mathbf{A} \) has real part \( \ge 0 \), the covariance may fail to converge (or may diverge).

Scalar special case: for

\[ \dot{x}(t) = -a x(t) + b\mu_u + b w(t), \qquad a > 0,\qquad E[w(t)w(\tau)] = q\,\delta(t-\tau) \]

the mean and variance obey

\[ \dot{\mu}_x = -a\mu_x + b\mu_u, \qquad \dot{P} = -2aP + b^2 q \]

with steady-state solutions

\[ \mu_{x,\infty} = \frac{b\mu_u}{a}, \qquad P_{\infty} = \frac{b^2 q}{2a} \]

For a stable LTI system and a wide-sense stationary (WSS) zero-mean input, Lesson 2 showed that the PSD describes how power is distributed across frequency. If the transfer function is \( G(s) \), then the output PSD is

\[ S_y(\omega) = \left|G(j\omega)\right|^2 S_u(\omega) \]

By the Wiener–Khintchine relation, the output variance is the zero-lag autocorrelation:

\[ \sigma_y^2 = R_y(0) = \frac{1}{2\pi}\int_{-\infty}^{\infty} S_y(\omega)\,d\omega = \frac{1}{2\pi}\int_{-\infty}^{\infty}\left|G(j\omega)\right|^2 S_u(\omega)\,d\omega \]

This formula is the frequency-domain counterpart of covariance propagation.

Consistency check for the first-order low-pass system:

\[ G(s)=\frac{b}{s+a}, \qquad a > 0 \]

For white input with PSD \( S_u(\omega)=q \),

\[ \sigma_y^2 = \frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{b^2 q}{a^2+\omega^2}\,d\omega = \frac{b^2 q}{2\pi}\cdot\frac{\pi}{a} = \frac{b^2 q}{2a} \]

exactly matching the time-domain Lyapunov result in Section 5.

Because Chapter 16 already introduced discrete-time dynamics, it is useful to state the stochastic propagation rules in discrete time as well:

\[ \mathbf{x}_{k+1}=\mathbf{A}_d\mathbf{x}_k+\mathbf{B}_d\mathbf{u}_k+\mathbf{G}_d\mathbf{w}_k \]

\[ E[\mathbf{w}_k]=\mathbf{0}, \qquad E[\mathbf{w}_k\mathbf{w}_j^T]= \begin{cases} \mathbf{Q}_d, & k=j \\ \mathbf{0}, & k\neq j \end{cases} \]

If \( \mathbf{w}_k \) is independent of \( \mathbf{x}_k \), then

\[ \boldsymbol{\mu}_{x,k+1} = \mathbf{A}_d\boldsymbol{\mu}_{x,k} + \mathbf{B}_d\boldsymbol{\mu}_{u,k} \]

\[ \mathbf{P}_{k+1} = \mathbf{A}_d\mathbf{P}_k\mathbf{A}_d^T + \mathbf{G}_d\mathbf{Q}_d\mathbf{G}_d^T \]

For scalar \( x_{k+1}=a x_k + g w_k \) with variance \( \sigma_w^2 \), the variance recursion is

\[ P_{k+1}=a^2 P_k + g^2 \sigma_w^2 \]

and if \( |a| < 1 \), the steady-state variance is

\[ P_{\infty}=\frac{g^2 \sigma_w^2}{1-a^2} \]

flowchart TD
  A["Random input model"] --> B["Mean recursion"]
  A --> C["Covariance recursion"]
  B --> D["Transient response of expected trajectory"]
  C --> E["Transient/steady uncertainty"]
  D --> F["Output mean and confidence bands"]
  E --> F
        

Consider the second-order mechanical model (from earlier dynamics chapters)

\[ m\ddot{x} + c\dot{x} + kx = f(t) \]

with random force \( f(t)=\bar{f}+w(t) \), where \( E[w(t)]=0 \) and \( E[w(t)w(\tau)] = q\delta(t-\tau) \). Defining state \( \mathbf{x}=[x,\dot{x}]^T \), we obtain

\[ \dot{\mathbf{x} } = \begin{bmatrix} 0 & 1 \\ -k/m & -c/m \end{bmatrix}\mathbf{x} + \begin{bmatrix} 0 \\ 1/m \end{bmatrix}\bar{f} + \begin{bmatrix} 0 \\ 1/m \end{bmatrix}w(t) \]

Therefore, the mean and covariance equations are

\[ \dot{\boldsymbol{\mu} }_x = \begin{bmatrix} 0 & 1 \\ -k/m & -c/m \end{bmatrix}\boldsymbol{\mu}_x + \begin{bmatrix} 0 \\ 1/m \end{bmatrix}\bar{f} \]

\[ \dot{\mathbf{P} } = \mathbf{A}\mathbf{P}+\mathbf{P}\mathbf{A}^{T} + \begin{bmatrix} 0 \\ 1/m \end{bmatrix} q \begin{bmatrix} 0 & 1/m \end{bmatrix} \]

The displacement variance at time \( t \) is the \( (1,1) \)-entry of \( \mathbf{P}(t) \). This is exactly what the Python and MATLAB codes compute and compare against Monte Carlo simulation.

Engineering interpretation: the damping \( c \) does not only shape the mean transient; it also dissipates stochastic energy, reducing steady-state covariance. Lower damping typically implies higher displacement variance under the same white-noise force intensity.

The following code examples implement mean/variance propagation and Monte Carlo verification in multiple languages used in control engineering workflows.

9.1 Python — Chapter17_Lesson3.py

# Chapter17_Lesson3.py
# Response of Linear Systems to Random Inputs: Mean and Variance Propagation
# Python implementation: continuous-time covariance propagation + Monte Carlo verification

import numpy as np

np.random.seed(7)

# Mass-spring-damper (state-space)
# x1 = displacement, x2 = velocity
m = 1.0
c = 0.6
k = 4.0

A = np.array([[0.0, 1.0],
              [-k/m, -c/m]])
Gamma = np.array([[0.0],
                  [1.0/m]])   # white-noise force channel
C = np.array([[1.0, 0.0]])    # observe displacement

# Deterministic mean input enters through B*u_mean
B = np.array([[0.0],
              [1.0/m]])
u_mean = 1.0  # constant mean force

# White noise force spectral intensity (continuous-time)
q = 0.8
Q = np.array([[q]])

# Simulation horizon
T = 20.0
dt = 1e-3
N = int(T / dt)

# Mean and covariance propagation (Euler integration of moment ODEs)
mx = np.zeros((2, 1))
P = np.zeros((2, 2))

mean_hist = np.zeros((N, 2))
var_x1_hist = np.zeros(N)

for n in range(N):
    mdot = A @ mx + B * u_mean
    Pdot = A @ P + P @ A.T + Gamma @ Q @ Gamma.T
    mx = mx + dt * mdot
    P = P + dt * Pdot

    mean_hist[n, :] = mx[:, 0]
    var_x1_hist[n] = P[0, 0]

# Monte Carlo (Euler-Maruyama) to verify the covariance result
M = 4000  # trajectories
x = np.zeros((M, 2))
x1_mean_mc = np.zeros(N)
x1_var_mc = np.zeros(N)

sqrt_qdt = np.sqrt(q * dt)

for n in range(N):
    # white-noise increment: w(t) dt ~ sqrt(q dt) * N(0,1)
    noise = np.random.randn(M, 1)
    drift = (x @ A.T) + (u_mean * np.ones((M, 1))) @ B.T
    diffusion = noise @ (sqrt_qdt * Gamma.T)
    x = x + dt * drift + diffusion

    x1 = x[:, 0]
    x1_mean_mc[n] = np.mean(x1)
    x1_var_mc[n] = np.var(x1, ddof=1)

# Steady-state covariance by solving the continuous Lyapunov equation numerically
# A P + P A^T + Gamma Q Gamma^T = 0
# We vectorize for a 2x2 system: vec(AP + PA^T) = (I kron A + A kron I) vec(P)
I2 = np.eye(2)
L = np.kron(I2, A) + np.kron(A, I2)
rhs = -(Gamma @ Q @ Gamma.T).reshape(-1, 1)
vecP_ss = np.linalg.solve(L, rhs)
P_ss = vecP_ss.reshape(2, 2)

print("Final propagated mean (theory):", mean_hist[-1, :])
print("Final propagated displacement variance (theory):", var_x1_hist[-1])
print("Final displacement mean (Monte Carlo):", x1_mean_mc[-1])
print("Final displacement variance (Monte Carlo):", x1_var_mc[-1])
print("Steady-state covariance from Lyapunov solve:")
print(P_ss)

# A compact discrete-time example (scalar) for lecture continuity
# x[k+1] = a_d x[k] + b_d * eta[k], eta[k] ~ N(0, sigma_eta^2)
a_d = 0.92
b_d = 1.0
sigma_eta2 = 0.5
P_k = 0.0
for kstep in range(50):
    P_k = a_d**2 * P_k + b_d**2 * sigma_eta2
print("Scalar discrete-time variance after 50 steps:", P_k)
print("Scalar steady-state variance:", (b_d**2 * sigma_eta2) / (1 - a_d**2))

9.2 C++ — Chapter17_Lesson3.cpp

// Chapter17_Lesson3.cpp
// Response of Linear Systems to Random Inputs: Mean and Variance Propagation
// C++ implementation: discrete-time scalar and 2-state examples (Monte Carlo + theory)

#include <iostream>
#include <vector>
#include <random>
#include <iomanip>
#include <cmath>

int main() {
    std::cout << std::fixed << std::setprecision(6);

    // Example 1: Scalar discrete-time system
    // x_{k+1} = a x_k + b mu_u + g w_k,  w_k ~ N(0, q)
    const double a = 0.95;
    const double b = 0.40;
    const double g = 1.0;
    const double mu_u = 1.2;
    const double q = 0.30;
    const int K = 120;
    const int M = 20000;

    // Theoretical mean and variance propagation
    double m = 0.0;
    double P = 0.0;
    for (int k = 0; k < K; ++k) {
        m = a * m + b * mu_u;
        P = a * a * P + g * g * q;
    }

    // Monte Carlo verification
    std::mt19937 rng(42);
    std::normal_distribution<double> N01(0.0, 1.0);

    std::vector<double> x(M, 0.0);
    for (int k = 0; k < K; ++k) {
        for (int i = 0; i < M; ++i) {
            double wk = std::sqrt(q) * N01(rng);
            x[i] = a * x[i] + b * mu_u + g * wk;
        }
    }

    double mean_mc = 0.0;
    for (double xi : x) mean_mc += xi;
    mean_mc /= M;

    double var_mc = 0.0;
    for (double xi : x) {
        double d = xi - mean_mc;
        var_mc += d * d;
    }
    var_mc /= (M - 1);

    double P_inf = (g * g * q) / (1.0 - a * a);
    double m_inf = (b * mu_u) / (1.0 - a);

    std::cout << "Scalar example (after " << K << " steps)\n";
    std::cout << "Theory mean      = " << m << "\n";
    std::cout << "Monte Carlo mean = " << mean_mc << "\n";
    std::cout << "Theory variance  = " << P << "\n";
    std::cout << "Monte Carlo var  = " << var_mc << "\n";
    std::cout << "Steady mean      = " << m_inf << "\n";
    std::cout << "Steady variance  = " << P_inf << "\n\n";

    // Example 2: 2-state stable linear system with white process noise
    // x_{k+1} = A x_k + G w_k, w_k ~ N(0, Q), mean zero input
    double A[2][2] = { {0.90, 0.10},
                      {-0.20, 0.80} };
    double G[2][1] = { {0.0},
                      {1.0} };
    double Q = 0.20;

    // Propagate covariance P_{k+1} = A P A^T + G Q G^T
    double Pm[2][2] = { {0.0, 0.0},
                       {0.0, 0.0} };

    for (int k = 0; k < 200; ++k) {
        // T = A P
        double T[2][2];
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                T[i][j] = A[i][0] * Pm[0][j] + A[i][1] * Pm[1][j];
            }
        }

        // APAT = T A^T
        double APAT[2][2];
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                APAT[i][j] = T[i][0] * A[j][0] + T[i][1] * A[j][1];
            }
        }

        // G Q G^T
        double GQGT[2][2] = {
            {G[0][0] * Q * G[0][0], G[0][0] * Q * G[1][0]},
            {G[1][0] * Q * G[0][0], G[1][0] * Q * G[1][0]}
        };

        for (int i = 0; i < 2; ++i)
            for (int j = 0; j < 2; ++j)
                Pm[i][j] = APAT[i][j] + GQGT[i][j];
    }

    std::cout << "2-state covariance after 200 steps (theory recursion):\n";
    std::cout << "[" << Pm[0][0] << ", " << Pm[0][1] << "]\n";
    std::cout << "[" << Pm[1][0] << ", " << Pm[1][1] << "]\n";

    return 0;
}

9.3 Java — Chapter17_Lesson3.java

// Chapter17_Lesson3.java
// Response of Linear Systems to Random Inputs: Mean and Variance Propagation
// Java implementation: discrete-time mean/variance recursion + Monte Carlo

import java.util.Random;

public class Chapter17_Lesson3 {
    public static void main(String[] args) {
        final double a = 0.93;
        final double b = 0.50;
        final double g = 1.0;
        final double muU = 0.8;
        final double q = 0.25;      // variance of white sequence w_k
        final int K = 150;
        final int M = 30000;

        // Theoretical propagation
        double m = 0.0;
        double P = 0.0;
        for (int k = 0; k < K; k++) {
            m = a * m + b * muU;
            P = a * a * P + g * g * q;
        }

        // Monte Carlo
        Random rng = new Random(1234);
        double[] x = new double[M];
        for (int k = 0; k < K; k++) {
            for (int i = 0; i < M; i++) {
                double wk = Math.sqrt(q) * rng.nextGaussian();
                x[i] = a * x[i] + b * muU + g * wk;
            }
        }

        double meanMC = 0.0;
        for (double xi : x) meanMC += xi;
        meanMC /= M;

        double varMC = 0.0;
        for (double xi : x) {
            double d = xi - meanMC;
            varMC += d * d;
        }
        varMC /= (M - 1);

        double steadyMean = (b * muU) / (1.0 - a);
        double steadyVar = (g * g * q) / (1.0 - a * a);

        System.out.printf("After %d steps:%n", K);
        System.out.printf("Theory mean      = %.6f%n", m);
        System.out.printf("Monte Carlo mean = %.6f%n", meanMC);
        System.out.printf("Theory variance  = %.6f%n", P);
        System.out.printf("Monte Carlo var  = %.6f%n", varMC);
        System.out.printf("Steady mean      = %.6f%n", steadyMean);
        System.out.printf("Steady variance  = %.6f%n", steadyVar);

        // Optional extension: covariance recursion for 2x2 systems
        double[][] A = { {0.88, 0.12}, {-0.15, 0.84} };
        double[][] G = { {0.0}, {1.0} };
        double Q = 0.10;
        double[][] P2 = { {0.0, 0.0}, {0.0, 0.0} };

        for (int step = 0; step < 250; step++) {
            double[][] AP = matMul(A, P2);
            double[][] APAT = matMul(AP, transpose(A));
            double[][] GQGT = {
                {G[0][0] * Q * G[0][0], G[0][0] * Q * G[1][0]},
                {G[1][0] * Q * G[0][0], G[1][0] * Q * G[1][0]}
            };
            P2 = matAdd(APAT, GQGT);
        }

        System.out.println("\n2x2 covariance after 250 steps:");
        System.out.printf("[%.6f, %.6f]%n", P2[0][0], P2[0][1]);
        System.out.printf("[%.6f, %.6f]%n", P2[1][0], P2[1][1]);
    }

    static double[][] matMul(double[][] A, double[][] B) {
        int r = A.length;
        int c = B[0].length;
        int n = B.length;
        double[][] C = new double[r][c];
        for (int i = 0; i < r; i++) {
            for (int j = 0; j < c; j++) {
                double s = 0.0;
                for (int k = 0; k < n; k++) s += A[i][k] * B[k][j];
                C[i][j] = s;
            }
        }
        return C;
    }

    static double[][] transpose(double[][] A) {
        double[][] T = new double[A[0].length][A.length];
        for (int i = 0; i < A.length; i++)
            for (int j = 0; j < A[0].length; j++)
                T[j][i] = A[i][j];
        return T;
    }

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

9.4 MATLAB/Simulink — Chapter17_Lesson3.m

% Chapter17_Lesson3.m
% Response of Linear Systems to Random Inputs: Mean and Variance Propagation
% MATLAB/Simulink-oriented implementation (moment equations + Monte Carlo)

clear; clc; rng(10);

% Continuous-time mass-spring-damper
m = 1.0; c = 0.5; k = 3.0;
A = [0 1; -k/m -c/m];
B = [0; 1/m];          % deterministic mean-force channel
Gamma = [0; 1/m];      % white-noise channel
C = [1 0];

u_mean = 1.0;
q = 0.7;               % white-noise spectral intensity

% Time integration for mean/covariance ODEs
dt = 1e-3; T = 15; N = round(T/dt);
mx = [0;0];
P  = zeros(2);

mx_hist = zeros(2,N);
var_hist = zeros(1,N);

for n = 1:N
    mdot = A*mx + B*u_mean;
    Pdot = A*P + P*A' + Gamma*q*Gamma';
    mx = mx + dt*mdot;
    P  = P  + dt*Pdot;

    mx_hist(:,n) = mx;
    var_hist(n) = P(1,1);
end

% Steady-state covariance via Lyapunov equation
% A*Pss + Pss*A' + Gamma*q*Gamma' = 0
Pss = lyap(A, Gamma*q*Gamma');

fprintf('Final propagated mean (theory): [%g, %g]\n', mx(1), mx(2));
fprintf('Final displacement variance (theory): %g\n', P(1,1));
disp('Steady-state covariance Pss (Lyapunov solution):');
disp(Pss);

% Monte Carlo verification (Euler-Maruyama)
M = 3000;
X = zeros(2,M);
sqrt_qdt = sqrt(q*dt);

for n = 1:N
    W = randn(1,M);
    drift = A*X + B*u_mean;
    diff  = Gamma * (sqrt_qdt * W);
    X = X + dt*drift + diff;
end

x1 = X(1,:);
fprintf('Monte Carlo displacement mean: %g\n', mean(x1));
fprintf('Monte Carlo displacement variance: %g\n', var(x1,1));

% -------------------------------------------------------------
% Simulink recipe (programmatic outline)
% -------------------------------------------------------------
% The following commands create a simple Simulink model with:
% Constant mean force + Band-Limited White Noise -> Sum -> State-Space block
%
% Uncomment and run in MATLAB with Simulink installed.
%
% modelName = 'Chapter17_Lesson3_Simulink';
% new_system(modelName); open_system(modelName);
% add_block('simulink/Sources/Constant', [modelName '/MeanForce'], ...
%           'Value', num2str(u_mean), 'Position', [30 40 80 70]);
% add_block('simulink/Sources/Band-Limited White Noise', [modelName '/Noise'], ...
%           'NoisePower', 'q', 'SampleTime', num2str(dt), 'Position', [30 100 160 130]);
% add_block('simulink/Math Operations/Sum', [modelName '/Sum'], ...
%           'Inputs', '++', 'Position', [210 55 240 115]);
% add_block('simulink/Continuous/State-Space', [modelName '/Plant'], ...
%           'A', mat2str(A), 'B', mat2str(B), 'C', mat2str(C), 'D', '0', ...
%           'Position', [300 50 450 120]);
% add_block('simulink/Sinks/Scope', [modelName '/Scope'], ...
%           'Position', [500 60 560 110]);
% add_line(modelName, 'MeanForce/1', 'Sum/1');
% add_line(modelName, 'Noise/1', 'Sum/2');
% add_line(modelName, 'Sum/1', 'Plant/1');
% add_line(modelName, 'Plant/1', 'Scope/1');
% set_param(modelName, 'StopTime', num2str(T));
% save_system(modelName);

9.5 Wolfram Mathematica — Chapter17_Lesson3.nb

(* Chapter17_Lesson3.nb *)
(* Response of Linear Systems to Random Inputs: Mean and Variance Propagation *)
(* Wolfram Language code (can be pasted directly into a Mathematica notebook) *)

ClearAll["Global`*"];

(* Continuous-time Lyapunov example *)
A = { {0, 1}, {-4, -0.6} };
Gamma = { {0}, {1} };
q = 0.8;

(* Solve A.P + P.Transpose[A] + Gamma q Gamma^T == 0 *)
Pss = LyapunovSolve[A, -Gamma.q.Transpose[Gamma]];
Print["Steady-state covariance Pss = "];
Print[MatrixForm[Pss]];

(* Deterministic mean propagation for constant mean input *)
B = { {0}, {1} };
uMean = 1.0;
mx0 = {0, 0};

(* Closed-form mean: m(t) = exp(A t) m0 + integral_0^t exp(A (t-tau)) B uMean d tau *)
m[t_] := MatrixExp[A t].mx0 + NIntegrate[MatrixExp[A (t - tau)].B[[All, 1]] uMean, {tau, 0, t}];
Print["Mean state at t=10: ", N[m[10], 8]];

(* Discrete-time scalar recursion and Monte Carlo check *)
a = 0.94; b = 0.4; g = 1.0; mu = 1.0; qk = 0.3;
K = 120; M = 20000;

mth = 0.0; Pth = 0.0;
Do[
  mth = a mth + b mu;
  Pth = a^2 Pth + g^2 qk;
, {k, 1, K}];

SeedRandom[1234];
x = ConstantArray[0.0, M];
Do[
  w = Sqrt[qk] RandomVariate[NormalDistribution[0, 1], M];
  x = a x + b mu + g w;
, {kk, 1, K}];

Print["Theory mean/variance after K steps: ", {mth, Pth}];
Print["Monte Carlo mean/variance after K steps: ", N[{Mean[x], Variance[x]}, 8]];
Print["Steady-state mean/variance: ", N[{(b mu)/(1 - a), (g^2 qk)/(1 - a^2)}, 8]];

Problem 1 (Mean propagation in a scalar linear system): Consider \( \dot{x}(t)=-a x(t)+b u(t) \) with \( a > 0 \), and \( u(t)=\mu_u + \tilde{u}(t) \), where \( E[\tilde{u}(t)]=0 \). Derive an ODE for \( \mu_x(t)=E[x(t)] \) and the steady-state mean.

Solution: Taking expectation gives

\[ \dot{\mu}_x(t) = -a\mu_x(t) + b\mu_u \]

because \( E[\tilde{u}(t)] = 0 \). The steady-state value solves \( 0 = -a\mu_{x,\infty} + b\mu_u \), hence

\[ \mu_{x,\infty} = \frac{b\mu_u}{a} \]

Problem 2 (Variance propagation for white-noise forcing): For \( \dot{x}(t)=-a x(t)+b w(t) \) with \( a > 0 \), \( E[w(t)]=0 \), and \( E[w(t)w(\tau)] = q\delta(t-\tau) \), derive the variance ODE for \( P(t)=E[x(t)^2] \).

Solution: The scalar covariance equation is a special case of Section 4:

\[ \dot{P}(t)=(-a)P + P(-a) + b^2 q = -2aP(t)+b^2 q \]

Therefore the variance converges exponentially to \( P_{\infty}=b^2 q/(2a) \).

Problem 3 (Frequency-domain verification of variance): For the same system in Problem 2, verify the steady-state variance using the PSD method.

Solution: The transfer function from \( w \) to \( x \) is

\[ G(s)=\frac{b}{s+a} \]

and for white input \( S_w(\omega)=q \):

\[ S_x(\omega)=\frac{b^2 q}{a^2+\omega^2} \]

Thus

\[ \sigma_x^2=\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{b^2 q}{a^2+\omega^2}\,d\omega = \frac{b^2 q}{2a} \]

which matches the time-domain covariance result exactly.

Problem 4 (Discrete-time variance recursion): Let \( x_{k+1}=a x_k + g w_k \), \( E[w_k]=0 \), \( E[w_k^2]=\sigma_w^2 \), and assume \( w_k \) is independent of \( x_k \). Derive the variance recursion and steady-state variance.

Solution: Since \( E[x_k]=0 \) (if initialized at zero mean),

\[ P_{k+1}=E[x_{k+1}^2]=E[(a x_k + g w_k)^2] = a^2 E[x_k^2] + 2agE[x_k w_k] + g^2 E[w_k^2] \]

Independence implies \( E[x_k w_k]=0 \), so

\[ P_{k+1}=a^2 P_k + g^2 \sigma_w^2 \]

If \( |a| < 1 \), then

\[ P_{\infty}=\frac{g^2 \sigma_w^2}{1-a^2} \]

Problem 5 (Matrix Lyapunov equation for a 2-state system): Let

\[ \mathbf{A}=\begin{bmatrix}0 & 1\\ -4 & -0.6\end{bmatrix},\qquad \mathbf{\Gamma}=\begin{bmatrix}0\\1\end{bmatrix},\qquad \mathbf{Q}=[q] \]

with \( q > 0 \). Write the steady-state covariance equation and explain how to solve it.

Solution: The steady-state covariance satisfies

\[ \mathbf{A}\mathbf{P}_{\infty}+\mathbf{P}_{\infty}\mathbf{A}^{T} + \mathbf{\Gamma}q\mathbf{\Gamma}^{T} = \mathbf{0} \]

Let \( \mathbf{P}_{\infty}=\begin{bmatrix}p_{11}&p_{12}\\p_{12}&p_{22}\end{bmatrix} \). Substituting into the matrix equation produces three linear equations in \( p_{11}, p_{12}, p_{22} \). Solve this linear system directly or use a numerical Lyapunov solver (as done in the Python, MATLAB, and Mathematica codes).

In this lesson, we established that linear systems propagate the mean of a random input through the same deterministic dynamics used in standard system analysis. We then derived covariance propagation equations for state-space models, including the differential and algebraic Lyapunov equations for white-process forcing. Finally, we connected the time-domain covariance view to the frequency-domain PSD formula and validated the theory with multi-language Monte Carlo simulations.

These results are foundational for noise analysis, random vibration, sensor/actuator uncertainty modeling (next lesson), and later topics such as estimation and filtering.