Deterministic system dynamics models assume that inputs, parameters, and initial conditions are known exactly. In practice, excitation forces, sensor noise, manufacturing tolerances, and environmental disturbances are uncertain. To model these uncertainties rigorously, we use \( \Omega \) (sample space), random variables \( X:\Omega \to \mathbb{R} \), and random processes \( X(t,\omega) \).

In random vibration and stochastic dynamics, the key object is a time-indexed family of random variables. At each fixed time \( t \), \( X(t) \) is a random variable; for each fixed outcome \( \omega \), \( t \mapsto X(t,\omega) \) is a sample path.

\[ X : T \times \Omega \to \mathbb{R}, \qquad (t,\omega) \mapsto X(t,\omega) \]

We will distinguish:

• Random variable level: distributions, means, variances, moments.
• Random process level: finite-dimensional distributions, mean function, covariance function.
• Stationarity level: invariance of statistics under time shift.

flowchart TD
  A["Uncertain source in dynamics (force, noise, parameter)"] --> B["Model as random variable or random process"]
  B --> C["Define ensemble statistics: mean, variance, covariance"]
  C --> D["Check time-shift invariance"]
  D --> E["Classify: nonstationary / WSS / strict stationary"]
  E --> F["Use in later lessons: PSD and random response"]
        

A random variable is a measurable mapping from the probability space \( (\Omega,\mathcal{F},\mathbb{P}) \) to \( \mathbb{R} \). For a random variable \( X \), the cumulative distribution function (CDF) is \( F_X(x)=\mathbb{P}(X \leq x) \). If \( X \) is continuous, its probability density function (PDF) \( f_X(x) \) satisfies

\[ F_X(x) = \int_{-\infty}^{x} f_X(\xi)\,d\xi, \qquad \int_{-\infty}^{\infty} f_X(x)\,dx = 1 \]

The expectation operator is the probabilistic analog of averaging over infinitely many independent repetitions. For any integrable function \( g \),

\[ \mathbb{E}[g(X)] = \int_{-\infty}^{\infty} g(x)\,f_X(x)\,dx \]

The first two central moments are particularly important for engineering uncertainty quantification:

\[ \mu_X = \mathbb{E}[X], \qquad \sigma_X^2 = \operatorname{Var}(X) = \mathbb{E}\!\left[(X-\mu_X)^2\right] \]

For two random variables \( X \) and \( Y \), covariance is

\[ \operatorname{Cov}(X,Y) = \mathbb{E}\!\left[(X-\mu_X)(Y-\mu_Y)\right] = \mathbb{E}[XY]-\mu_X\mu_Y \]

The covariance matrix of a random vector \( \mathbf{X}\in\mathbb{R}^n \) is

\[ \mathbf{C}_{\mathbf{X}} = \mathbb{E}\!\left[(\mathbf{X}-\boldsymbol{\mu})(\mathbf{X}-\boldsymbol{\mu})^T\right] \]

This matrix is always symmetric positive semidefinite: \( \mathbf{a}^T \mathbf{C}_{\mathbf{X}} \mathbf{a} \geq 0 \) for every \( \mathbf{a}\in\mathbb{R}^n \). This positivity property will reappear later as a crucial constraint on valid covariance functions of random processes.

A random process is a family of random variables indexed by time (continuous or discrete). We write \( \{X(t): t \in T\} \), where \( T=\mathbb{R} \) for continuous-time processes and \( T=\mathbb{Z} \) (or \( \mathbb{N} \)) for discrete-time processes.

The most complete probabilistic description is given by the finite-dimensional distributions: for any \( n \in \mathbb{N} \) and times \( t_1,\dots,t_n \in T \),

\[ F_{X(t_1),\dots,X(t_n)}(x_1,\dots,x_n) = \mathbb{P}\!\left(X(t_1)\leq x_1,\dots,X(t_n)\leq x_n\right) \]

In engineering practice, we often work with lower-order summaries:

• Mean function: \( m_X(t)=\mathbb{E}[X(t)] \)
• Second moment: \( M_X(t_1,t_2)=\mathbb{E}[X(t_1)X(t_2)] \)
• Covariance function: \( C_X(t_1,t_2)=\operatorname{Cov}(X(t_1),X(t_2)) \)

\[ C_X(t_1,t_2)=M_X(t_1,t_2)-m_X(t_1)m_X(t_2) \]

A process is second-order if \( \mathbb{E}[X(t)^2] < \infty \) for all \( t \). Most random vibration models studied in control and mechanical systems are second-order, because mean-square quantities correspond directly to physical energy-like measures.

In experiments, we often observe only one realization (one sample path) of a process. The theory, however, is built on ensemble averages over repeated experiments. For a sample path \( x(t) \) observed on \( [0,T] \), the time average is

\[ \overline{x}_T = \frac{1}{T}\int_{0}^{T} x(t)\,dt \]

The corresponding ensemble mean at time \( t \) is \( m_X(t)=\mathbb{E}[X(t)] \). Equality of time and ensemble averages is an ergodicity question, which is stronger than stationarity and will be used carefully later. In this lesson, the key point is that stationarity is defined by ensemble invariance, not by visual appearance of a single record.

For numerical work, ensemble estimates from \( M \) independent simulations are

\[ \widehat{m}_X(t) = \frac{1}{M}\sum_{r=1}^{M} X^{(r)}(t), \qquad \widehat{C}_X(t_1,t_2)=\frac{1}{M}\sum_{r=1}^{M}\left(X^{(r)}(t_1)-\widehat{m}_X(t_1)\right)\left(X^{(r)}(t_2)-\widehat{m}_X(t_2)\right) \]

These estimators appear in the code examples below and form the basis of Monte Carlo analysis in Chapter 17, Lesson 5.

Stationarity is a time-shift invariance property. The strongest common form is strict-sense stationarity (SSS), which requires invariance of all finite-dimensional distributions under any common time shift.

Definition (SSS): The process \( X(t) \) is strictly stationary if for every \( n \in \mathbb{N} \), every \( t_1,\dots,t_n \in T \), every shift \( \tau \in T \), and every \( x_1,\dots,x_n \in \mathbb{R} \),

\[ \mathbb{P}\!\left(X(t_1)\leq x_1,\dots,X(t_n)\leq x_n\right) = \mathbb{P}\!\left(X(t_1+\tau)\leq x_1,\dots,X(t_n+\tau)\leq x_n\right) \]

In many engineering applications, second-order statistics are sufficient. Then we use wide-sense stationarity (WSS), also called weak stationarity.

Definition (WSS): A second-order process \( X(t) \) is WSS if

1. \( m_X(t)=\mu \) is constant, independent of \( t \), and
2. \( C_X(t_1,t_2)=C_X(t_1+\tau,t_2+\tau) \) for every shift \( \tau \).

Under WSS, covariance depends only on lag \( \eta=t_2-t_1 \), so we write

\[ C_X(t_1,t_2)=C_X(\eta), \qquad R_X(t_1,t_2)=R_X(\eta) \]

where \( R_X \) is the autocorrelation function \( R_X(t_1,t_2)=\mathbb{E}[X(t_1)X(t_2)] \). If the mean is constant \( \mu \),

\[ R_X(\eta)=C_X(\eta)+\mu^2 \]

Theorem 1 (SSS implies WSS under finite second moments). Suppose \( X(t) \) is strictly stationary and \( \mathbb{E}[X(t)^2] < \infty \) for all \( t \). Then \( X(t) \) is wide-sense stationary.

Proof: Since strict stationarity implies equality of all one-dimensional distributions under shifts, \( X(t) \) and \( X(t+\tau) \) have the same distribution. Therefore, whenever the expectation exists,

\[ \mathbb{E}[X(t+\tau)] = \mathbb{E}[X(t)] = \mu \]

so the mean is constant. Next, strict stationarity also implies equality of all two-dimensional distributions:

\[ (X(t_1),X(t_2)) \overset{d}{=} (X(t_1+\tau),X(t_2+\tau)) \]

Since second moments exist, equal distributions imply equality of the mixed second moment:

\[ \mathbb{E}[X(t_1)X(t_2)] = \mathbb{E}[X(t_1+\tau)X(t_2+\tau)] \]

Subtracting the constant mean product \( \mu^2 \) from both sides gives \( C_X(t_1,t_2)=C_X(t_1+\tau,t_2+\tau) \). Hence the process is WSS. \(\square\)

Theorem 2 (WSS Gaussian process is SSS). If \( X(t) \) is a Gaussian process and is WSS, then it is strictly stationary.

Proof sketch: A Gaussian process is fully determined by its finite-dimensional mean vectors and covariance matrices. Under WSS, the mean vector of \( (X(t_1),\dots,X(t_n)) \) is \( \mu \mathbf{1} \), which is invariant under a common shift. The covariance entries are

\[ \operatorname{Cov}(X(t_i),X(t_j)) = C_X(t_j-t_i) \]

and thus are unchanged when each time is shifted by the same \( \tau \). Therefore the shifted and unshifted Gaussian vectors have the same mean and covariance, hence the same distribution. This holds for all \( n \), so the process is strictly stationary. \(\square\)

Theorem 3 (Affine transformation preserves WSS). If \( X(t) \) is WSS and \( Y(t)=aX(t)+b \) with constants \( a,b \in \mathbb{R} \), then \( Y(t) \) is WSS.

Proof: The mean is constant:

\[ m_Y(t)=\mathbb{E}[aX(t)+b]=a\,m_X(t)+b=a\mu+b \]

For covariance,

\[ C_Y(t_1,t_2)=\operatorname{Cov}(aX(t_1)+b,aX(t_2)+b)=a^2 C_X(t_1,t_2) \]

Since \( C_X(t_1,t_2) \) depends only on time lag, so does \( C_Y(t_1,t_2) \). Therefore \( Y(t) \) is WSS. \(\square\)

Example 1 (Random-phase harmonic process). Let \( X(t)=A\cos(\omega_0 t + \Theta) \) where \( A \) and \( \omega_0 \) are constants and \( \Theta \sim \operatorname{Unif}(0,2\pi) \). Then the process is strictly stationary.

The mean is

\[ m_X(t)=\mathbb{E}[A\cos(\omega_0 t+\Theta)] = \frac{A}{2\pi}\int_{0}^{2\pi}\cos(\omega_0 t+\theta)\,d\theta = 0 \]

The autocorrelation is

\[ R_X(t_1,t_2)=\mathbb{E}[X(t_1)X(t_2)] = \frac{A^2}{2}\cos\!\left(\omega_0 (t_2-t_1)\right) \]

which depends only on lag. The strict-stationarity proof follows from the fact that \( \Theta+\omega_0\tau \pmod{2\pi} \) is also uniform on \( [0,2\pi) \).

Example 2 (Nonstationary random slope/intercept process). Let \( Y(t)=\alpha t+\beta \) with random constants \( \alpha,\beta \). Then

\[ m_Y(t)=\mathbb{E}[\alpha]\,t+\mathbb{E}[\beta] \]

\[ \operatorname{Var}(Y(t)) = t^2 \operatorname{Var}(\alpha) + 2t\,\operatorname{Cov}(\alpha,\beta) + \operatorname{Var}(\beta) \]

Unless very special parameter relations hold, these moments vary with \( t \), so the process is not WSS.

Example 3 (AR(1) discrete-time process). Consider \( x[k]=a\,x[k-1]+w[k] \) where \( w[k] \) is zero-mean white-noise innovation with variance \( \sigma_w^2 \). If \( |a| < 1 \) and the initial state is drawn from the invariant distribution, then the process is WSS with

\[ \mathbb{E}[x[k]]=0, \qquad \operatorname{Var}(x[k])=\frac{\sigma_w^2}{1-a^2} \]

This is a standard stochastic state update used in system identification, random disturbances, and sampled-data models.

In control engineering and vibration analysis, stochastic forcing enters differential equations as a random process. For a single-degree-of-freedom oscillator,

\[ m\ddot{q}(t)+c\dot{q}(t)+kq(t)=F(t) \]

If \( F(t) \) is deterministic, the response \( q(t) \) is deterministic. If \( F(t) \) is random, then \( q(t) \) is also a random process. Before computing response variance or spectral content, we must first characterize the input process:

• Is the mean constant or drifting?
• Does covariance depend only on lag?
• Is a Gaussian model reasonable?
• Is the process continuous-time or sampled (discrete-time)?

These questions determine whether later tools (autocorrelation, PSD, transfer-based variance propagation) can be applied directly.

The following programs generate three illustrative processes: (i) a random-phase harmonic process (stationary), (ii) a nonstationary affine-in-time process, and (iii) a discrete AR(1) process initialized in its stationary regime. They estimate ensemble mean and variance to numerically verify stationarity concepts.

Code: Chapter17_Lesson1.py

# Chapter17_Lesson1.py
# Random Variables, Random Processes, and Stationarity Concepts
# Python implementation for System Dynamics (Chapter 17, Lesson 1)

import numpy as np
import matplotlib.pyplot as plt

np.random.seed(17)

# ------------------------------------------------------------
# Example 1: Random-phase harmonic process (wide-sense stationary)
# X(t) = A cos(omega t + Theta), Theta ~ Uniform(0, 2*pi)
# ------------------------------------------------------------
A = 2.0
omega = 3.0
M = 2000                  # number of realizations
t = np.linspace(0.0, 4.0, 401)

theta = np.random.uniform(0.0, 2.0 * np.pi, size=M)
X = np.array([A * np.cos(omega * t + th) for th in theta])   # shape: (M, Nt)

# Ensemble mean estimate (function of time)
mean_est = X.mean(axis=0)

# Ensemble variance estimate (function of time)
var_est = X.var(axis=0, ddof=0)

# Autocorrelation estimate at a few lags (ensemble-based)
def ensemble_autocorr(X_ensemble, lag_index):
    x1 = X_ensemble[:, :-lag_index] if lag_index > 0 else X_ensemble
    x2 = X_ensemble[:, lag_index:]  if lag_index > 0 else X_ensemble
    return np.mean(x1 * x2)

lags = [0, 10, 30, 60]
dt = t[1] - t[0]
print("Random-phase harmonic process estimates:")
print("Mean at selected times:", mean_est[::100])
print("Variance at selected times:", var_est[::100])
for L in lags:
    tau = L * dt
    R_hat = ensemble_autocorr(X, L)
    R_theory = 0.5 * A**2 * np.cos(omega * tau)
    print(f"lag={L:3d}, tau={tau:.3f}, R_hat={R_hat:.5f}, R_theory={R_theory:.5f}")

# ------------------------------------------------------------
# Example 2: Nonstationary process
# Y(t) = alpha * t + beta, alpha and beta random (constant per realization)
# mean and variance depend on t
# ------------------------------------------------------------
alpha = np.random.normal(loc=0.5, scale=0.2, size=M)
beta = np.random.normal(loc=0.0, scale=1.0, size=M)
Y = np.array([a * t + b for a, b in zip(alpha, beta)])

Y_mean = Y.mean(axis=0)
Y_var = Y.var(axis=0, ddof=0)

print("\nNonstationary affine-in-time process estimates:")
for idx in [0, 100, 200, 400]:
    print(f"t={t[idx]:.2f}, mean={Y_mean[idx]:.4f}, var={Y_var[idx]:.4f}")

# ------------------------------------------------------------
# Example 3: Discrete-time AR(1) process (stationary regime if |a| < 1)
# x[k] = a*x[k-1] + w[k], w[k] ~ N(0, sigma_w^2)
# ------------------------------------------------------------
a = 0.8
sigma_w = 1.0
K = 400
M_ar = 500

# Initialize from invariant Gaussian distribution for stationarity:
# Var[x] = sigma_w^2 / (1 - a^2)
sigma_x0 = sigma_w / np.sqrt(1.0 - a**2)
x = np.zeros((M_ar, K))
x[:, 0] = np.random.normal(0.0, sigma_x0, size=M_ar)
for k in range(1, K):
    w = np.random.normal(0.0, sigma_w, size=M_ar)
    x[:, k] = a * x[:, k - 1] + w

mean_ar = x.mean(axis=0)
var_ar = x.var(axis=0, ddof=0)

print("\nAR(1) stationary regime estimates:")
print("Mean at k = [0, 100, 200, 399]:", [float(mean_ar[i]) for i in [0, 100, 200, 399]])
print("Var  at k = [0, 100, 200, 399]:", [float(var_ar[i]) for i in [0, 100, 200, 399]])
print("Theoretical stationary variance:", sigma_w**2 / (1.0 - a**2))

# ------------------------------------------------------------
# Visualization
# ------------------------------------------------------------
fig, ax = plt.subplots(2, 1, figsize=(9, 8), constrained_layout=True)

# Random-phase process statistics
ax[0].plot(t, mean_est, label="Estimated ensemble mean")
ax[0].plot(t, np.zeros_like(t), "--", label="Theory mean = 0")
ax[0].set_title("Random-phase harmonic process: ensemble mean")
ax[0].set_xlabel("t [s]")
ax[0].set_ylabel("Mean")
ax[0].grid(True)
ax[0].legend()

# Nonstationary process variance
ax[1].plot(t, Y_var, label="Estimated Var[Y(t)]")
ax[1].set_title("Nonstationary process Y(t)=alpha t + beta: variance changes with time")
ax[1].set_xlabel("t [s]")
ax[1].set_ylabel("Variance")
ax[1].grid(True)
ax[1].legend()

plt.show()

      

Code: Chapter17_Lesson1.cpp

/* Chapter17_Lesson1.cpp
   Random Variables, Random Processes, and Stationarity Concepts
   C++ implementation for System Dynamics (Chapter 17, Lesson 1)
*/

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

static double mean_of_vector(const std::vector<double>& v) {
    double s = std::accumulate(v.begin(), v.end(), 0.0);
    return s / static_cast<double>(v.size());
}

static double variance_of_vector(const std::vector<double>& v) {
    double m = mean_of_vector(v);
    double s2 = 0.0;
    for (double x : v) {
        double d = x - m;
        s2 += d * d;
    }
    return s2 / static_cast<double>(v.size());
}

int main() {
    std::mt19937 rng(17);
    std::uniform_real_distribution<double> unif(0.0, 2.0 * M_PI);
    std::normal_distribution<double> n_alpha(0.5, 0.2);
    std::normal_distribution<double> n_beta(0.0, 1.0);

    const double A = 2.0;
    const double omega = 3.0;
    const int Nt = 401;
    const int M = 2000;
    const double t0 = 0.0;
    const double t1 = 4.0;
    const double dt = (t1 - t0) / (Nt - 1);

    std::vector<double> t(Nt);
    for (int i = 0; i < Nt; ++i) {
        t[i] = t0 + i * dt;
    }

    // Ensemble means and variances for X(t) = A cos(omega t + Theta)
    std::vector<double> meanX(Nt, 0.0), meanX2(Nt, 0.0);

    for (int m = 0; m < M; ++m) {
        double theta = unif(rng);
        for (int i = 0; i < Nt; ++i) {
            double x = A * std::cos(omega * t[i] + theta);
            meanX[i] += x;
            meanX2[i] += x * x;
        }
    }
    for (int i = 0; i < Nt; ++i) {
        meanX[i] /= static_cast<double>(M);
        meanX2[i] /= static_cast<double>(M);
    }

    std::vector<double> varX(Nt, 0.0);
    for (int i = 0; i < Nt; ++i) {
        varX[i] = meanX2[i] - meanX[i] * meanX[i];
    }

    // Ensemble autocorrelation estimates at selected lags
    std::vector<int> lags = {0, 10, 30, 60};
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Random-phase harmonic process estimates\n";
    for (int L : lags) {
        double sum = 0.0;
        int count = 0;
        for (int m = 0; m < M; ++m) {
            double theta = unif(rng); // independent realizations for correlation estimate
            for (int i = 0; i < Nt - L; ++i) {
                double x1 = A * std::cos(omega * t[i] + theta);
                double x2 = A * std::cos(omega * t[i + L] + theta);
                sum += x1 * x2;
                ++count;
            }
        }
        double Rhat = sum / static_cast<double>(count);
        double tau = L * dt;
        double Rth = 0.5 * A * A * std::cos(omega * tau);
        std::cout << "lag=" << std::setw(3) << L
                  << ", tau=" << tau
                  << ", R_hat=" << Rhat
                  << ", R_theory=" << Rth << "\n";
    }

    std::cout << "\nSelected mean/variance samples for X(t):\n";
    for (int idx : {0, 100, 200, 400}) {
        std::cout << "t=" << t[idx]
                  << ", mean=" << meanX[idx]
                  << ", var=" << varX[idx] << "\n";
    }

    // Nonstationary process Y(t) = alpha t + beta
    std::vector<double> meanY(Nt, 0.0), meanY2(Nt, 0.0);
    for (int m = 0; m < M; ++m) {
        double alpha = n_alpha(rng);
        double beta = n_beta(rng);
        for (int i = 0; i < Nt; ++i) {
            double y = alpha * t[i] + beta;
            meanY[i] += y;
            meanY2[i] += y * y;
        }
    }
    for (int i = 0; i < Nt; ++i) {
        meanY[i] /= static_cast<double>(M);
        meanY2[i] /= static_cast<double>(M);
    }

    std::cout << "\nNonstationary affine-in-time process estimates:\n";
    for (int idx : {0, 100, 200, 400}) {
        double varY = meanY2[idx] - meanY[idx] * meanY[idx];
        std::cout << "t=" << t[idx]
                  << ", mean=" << meanY[idx]
                  << ", var=" << varY << "\n";
    }

    // AR(1) process in stationary regime
    const double a = 0.8;
    const double sigma_w = 1.0;
    const int K = 400;
    const int Mar = 500;
    std::normal_distribution<double> n_w(0.0, sigma_w);
    std::normal_distribution<double> n_x0(0.0, sigma_w / std::sqrt(1.0 - a * a));

    std::vector<double> meanAR(K, 0.0), meanAR2(K, 0.0);
    for (int m = 0; m < Mar; ++m) {
        std::vector<double> x(K, 0.0);
        x[0] = n_x0(rng);
        for (int k = 1; k < K; ++k) {
            x[k] = a * x[k - 1] + n_w(rng);
        }
        for (int k = 0; k < K; ++k) {
            meanAR[k] += x[k];
            meanAR2[k] += x[k] * x[k];
        }
    }
    for (int k = 0; k < K; ++k) {
        meanAR[k] /= static_cast<double>(Mar);
        meanAR2[k] /= static_cast<double>(Mar);
    }

    std::cout << "\nAR(1) stationary regime estimates:\n";
    for (int k : {0, 100, 200, 399}) {
        double varAR = meanAR2[k] - meanAR[k] * meanAR[k];
        std::cout << "k=" << k
                  << ", mean=" << meanAR[k]
                  << ", var=" << varAR << "\n";
    }
    std::cout << "Theoretical stationary variance = " << (sigma_w * sigma_w) / (1.0 - a * a) << "\n";

    return 0;
}

      

Code: Chapter17_Lesson1.java

/*
 * Chapter17_Lesson1.java
 * Random Variables, Random Processes, and Stationarity Concepts
 * Java implementation for System Dynamics (Chapter 17, Lesson 1)
 */

import java.util.Random;

public class Chapter17_Lesson1 {

    static double mean(double[] x) {
        double s = 0.0;
        for (double v : x) s += v;
        return s / x.length;
    }

    static double variance(double[] x) {
        double m = mean(x);
        double s2 = 0.0;
        for (double v : x) {
            double d = v - m;
            s2 += d * d;
        }
        return s2 / x.length;
    }

    static double gaussian(Random rng, double mu, double sigma) {
        return mu + sigma * rng.nextGaussian();
    }

    public static void main(String[] args) {
        Random rng = new Random(17);

        final double A = 2.0;
        final double omega = 3.0;
        final int Nt = 401;
        final int M = 2000;
        final double t0 = 0.0;
        final double t1 = 4.0;
        final double dt = (t1 - t0) / (Nt - 1);

        double[] t = new double[Nt];
        for (int i = 0; i < Nt; i++) {
            t[i] = t0 + i * dt;
        }

        // Random-phase harmonic process X(t)=A cos(omega t + Theta)
        double[] meanX = new double[Nt];
        double[] meanX2 = new double[Nt];

        for (int m = 0; m < M; m++) {
            double theta = 2.0 * Math.PI * rng.nextDouble();
            for (int i = 0; i < Nt; i++) {
                double x = A * Math.cos(omega * t[i] + theta);
                meanX[i] += x;
                meanX2[i] += x * x;
            }
        }
        for (int i = 0; i < Nt; i++) {
            meanX[i] /= M;
            meanX2[i] /= M;
        }

        System.out.println("Random-phase harmonic process estimates:");
        int[] sel = {0, 100, 200, 400};
        for (int idx : sel) {
            double varX = meanX2[idx] - meanX[idx] * meanX[idx];
            System.out.printf("t=%.2f, mean=%.6f, var=%.6f%n", t[idx], meanX[idx], varX);
        }

        // Correlation estimate at selected lags
        int[] lags = {0, 10, 30, 60};
        for (int L : lags) {
            double sum = 0.0;
            long count = 0L;
            for (int m = 0; m < M; m++) {
                double theta = 2.0 * Math.PI * rng.nextDouble();
                for (int i = 0; i < Nt - L; i++) {
                    double x1 = A * Math.cos(omega * t[i] + theta);
                    double x2 = A * Math.cos(omega * t[i + L] + theta);
                    sum += x1 * x2;
                    count++;
                }
            }
            double tau = L * dt;
            double Rhat = sum / count;
            double Rth = 0.5 * A * A * Math.cos(omega * tau);
            System.out.printf("lag=%d, tau=%.3f, Rhat=%.6f, Rtheory=%.6f%n", L, tau, Rhat, Rth);
        }

        // Nonstationary process Y(t)=alpha t + beta
        double[] meanY = new double[Nt];
        double[] meanY2 = new double[Nt];
        for (int m = 0; m < M; m++) {
            double alpha = gaussian(rng, 0.5, 0.2);
            double beta = gaussian(rng, 0.0, 1.0);
            for (int i = 0; i < Nt; i++) {
                double y = alpha * t[i] + beta;
                meanY[i] += y;
                meanY2[i] += y * y;
            }
        }
        for (int i = 0; i < Nt; i++) {
            meanY[i] /= M;
            meanY2[i] /= M;
        }

        System.out.println("\nNonstationary affine-in-time process estimates:");
        for (int idx : sel) {
            double varY = meanY2[idx] - meanY[idx] * meanY[idx];
            System.out.printf("t=%.2f, mean=%.6f, var=%.6f%n", t[idx], meanY[idx], varY);
        }

        // AR(1) process x[k] = a x[k-1] + w[k], stationary if |a| < 1
        final double a = 0.8;
        final double sigmaW = 1.0;
        final int K = 400;
        final int Mar = 500;

        double[] meanAR = new double[K];
        double[] meanAR2 = new double[K];
        double sigmaX0 = sigmaW / Math.sqrt(1.0 - a * a);

        for (int m = 0; m < Mar; m++) {
            double[] x = new double[K];
            x[0] = gaussian(rng, 0.0, sigmaX0);
            for (int k = 1; k < K; k++) {
                x[k] = a * x[k - 1] + gaussian(rng, 0.0, sigmaW);
            }
            for (int k = 0; k < K; k++) {
                meanAR[k] += x[k];
                meanAR2[k] += x[k] * x[k];
            }
        }
        for (int k = 0; k < K; k++) {
            meanAR[k] /= Mar;
            meanAR2[k] /= Mar;
        }

        System.out.println("\nAR(1) stationary regime estimates:");
        int[] ksel = {0, 100, 200, 399};
        for (int k : ksel) {
            double varAR = meanAR2[k] - meanAR[k] * meanAR[k];
            System.out.printf("k=%d, mean=%.6f, var=%.6f%n", k, meanAR[k], varAR);
        }
        System.out.printf("Theoretical stationary variance = %.6f%n", sigmaW * sigmaW / (1.0 - a * a));
    }
}

      

Code: Chapter17_Lesson1.m (MATLAB + optional Simulink model generation)

% Chapter17_Lesson1.m
% Random Variables, Random Processes, and Stationarity Concepts
% MATLAB / Simulink implementation for System Dynamics (Chapter 17, Lesson 1)

clear; clc; close all;
rng(17);

%% Example 1: Random-phase harmonic process (WSS)
% X(t) = A*cos(omega*t + Theta), Theta ~ U(0, 2*pi)
A = 2.0;
omega = 3.0;
M = 2000;
t = linspace(0, 4, 401);
Nt = numel(t);

theta = 2*pi*rand(M,1);
X = zeros(M, Nt);
for m = 1:M
    X(m,:) = A*cos(omega*t + theta(m));
end

meanX = mean(X, 1);
varX  = var(X, 1, 1);   % population variance (normalize by N)

lags = [0, 10, 30, 60];
dt = t(2)-t(1);
fprintf('Random-phase harmonic process estimates:\n');
for L = lags
    if L == 0
        x1 = X;
        x2 = X;
    else
        x1 = X(:,1:end-L);
        x2 = X(:,1+L:end);
    end
    Rhat = mean(x1(:).*x2(:));
    tau = L*dt;
    Rth = 0.5*A^2*cos(omega*tau);
    fprintf('lag=%3d, tau=%7.4f, Rhat=%9.5f, Rtheory=%9.5f\n', L, tau, Rhat, Rth);
end

%% Example 2: Nonstationary process Y(t) = alpha*t + beta
alpha = 0.5 + 0.2*randn(M,1);
beta  = randn(M,1);
Y = alpha*t + beta;

meanY = mean(Y,1);
varY  = var(Y,1,1);

fprintf('\nNonstationary process estimates (time-varying moments):\n');
sel = [1, 101, 201, 401];
for idx = sel
    fprintf('t=%6.3f, mean=%9.5f, var=%9.5f\n', t(idx), meanY(idx), varY(idx));
end

%% Example 3: Discrete AR(1) process in stationary regime
% x[k] = a*x[k-1] + w[k], stationarity requires |a| < 1
a = 0.8;
sigma_w = 1.0;
K = 400;
Mar = 500;

sigma_x0 = sigma_w/sqrt(1-a^2); % invariant std
Xar = zeros(Mar, K);
Xar(:,1) = sigma_x0*randn(Mar,1);

for k = 2:K
    Xar(:,k) = a*Xar(:,k-1) + sigma_w*randn(Mar,1);
end

meanAR = mean(Xar,1);
varAR  = var(Xar,1,1);

fprintf('\nAR(1) stationary regime estimates:\n');
kSel = [1, 101, 201, 400];
for idx = kSel
    fprintf('k=%3d, mean=%9.5f, var=%9.5f\n', idx-1, meanAR(idx), varAR(idx));
end
fprintf('Theoretical stationary variance = %9.5f\n', sigma_w^2/(1-a^2));

%% Plot results
figure('Name','Chapter17 Lesson1 - Stationarity Concepts');
subplot(2,1,1);
plot(t, meanX, 'LineWidth', 1.2); hold on;
plot(t, zeros(size(t)), '--', 'LineWidth', 1.0);
grid on;
xlabel('t [s]'); ylabel('Ensemble Mean');
title('Random-phase harmonic process: estimated ensemble mean');
legend('Estimated','Theory mean = 0','Location','best');

subplot(2,1,2);
plot(t, varY, 'LineWidth', 1.2);
grid on;
xlabel('t [s]'); ylabel('Variance');
title('Nonstationary process Y(t)=alpha t + beta: variance depends on time');

%% Optional Simulink model creation (AR(1)-like recursion)
% This block programmatically creates a simple Simulink diagram illustrating
% a stochastic recursion x[k] = a*x[k-1] + w[k]. It requires Simulink.
try
    mdl = 'Chapter17_Lesson1_Simulink';
    if bdIsLoaded(mdl), close_system(mdl, 0); end
    new_system(mdl); open_system(mdl);

    add_block('simulink/Sources/Random Number', [mdl '/w'], ...
        'Position', [40 80 120 110]);
    add_block('simulink/Math Operations/Gain', [mdl '/sigma_w'], ...
        'Gain', '1.0', 'Position', [150 75 220 115]);
    add_block('simulink/Math Operations/Sum', [mdl '/Sum'], ...
        'Inputs', '++', 'Position', [280 70 310 120]);
    add_block('simulink/Math Operations/Gain', [mdl '/a'], ...
        'Gain', '0.8', 'Position', [350 130 420 170]);
    add_block('simulink/Discrete/Unit Delay', [mdl '/z^-1'], ...
        'Position', [460 125 520 175], 'SampleTime', '0.01');
    add_block('simulink/Sinks/Scope', [mdl '/Scope'], ...
        'Position', [560 70 620 110]);

    add_line(mdl, 'w/1', 'sigma_w/1');
    add_line(mdl, 'sigma_w/1', 'Sum/1');
    add_line(mdl, 'Sum/1', 'Scope/1');
    add_line(mdl, 'Sum/1', 'a/1');
    add_line(mdl, 'a/1', 'z^-1/1');
    add_line(mdl, 'z^-1/1', 'Sum/2');

    set_param(mdl, 'StopTime', '10');
    fprintf('\nSimulink model "%s" created successfully.\n', mdl);
catch ME
    fprintf('\nSimulink model creation skipped: %s\n', ME.message);
end

      

Code: Chapter17_Lesson1.nb (Wolfram Mathematica input cells)

\
SeedRandom[17];

(* Example 1: Random-phase harmonic process *)
A = 2.0; omega = 3.0; m = 2000;
tGrid = N@Subdivide[0, 4, 400];
theta = RandomReal[{0, 2 Pi}, m];
xEnsemble = Table[A Cos[omega tGrid + theta[[i]]], {i, 1, m}];

meanX = Mean[xEnsemble];
varX = Variance /@ Transpose[xEnsemble];

rHat[L_Integer] := Module[{x1, x2},
  If[L == 0,
    Mean[Flatten[xEnsemble*xEnsemble]],
    x1 = xEnsemble[[All, 1 ;; -L - 1]];
    x2 = xEnsemble[[All, 1 + L ;; -1]];
    Mean[Flatten[x1*x2]]
  ]
];

lags = {0, 10, 30, 60};
dt = tGrid[[2]] - tGrid[[1]];
Print["Random-phase harmonic process estimates:"];
Do[
  tau = L*dt;
  rTheory = 0.5*A^2*Cos[omega*tau];
  Print[
    "lag=", L, ", tau=", NumberForm[tau, {5, 3}],
    ", Rhat=", NumberForm[rHat[L], {8, 5}],
    ", Rtheory=", NumberForm[rTheory, {8, 5}]
  ],
  {L, lags}
];

(* Example 2: Nonstationary process Y(t) = alpha t + beta *)
alpha = RandomVariate[NormalDistribution[0.5, 0.2], m];
beta = RandomVariate[NormalDistribution[0, 1.0], m];
yEnsemble = Table[alpha[[i]]*tGrid + beta[[i]], {i, 1, m}];
meanY = Mean[yEnsemble];
varY = Variance /@ Transpose[yEnsemble];

Print["Nonstationary affine-in-time process estimates:"];
Do[
  idx = k;
  Print[
    "t=", NumberForm[tGrid[[idx]], {5, 2}],
    ", mean=", NumberForm[meanY[[idx]], {8, 5}],
    ", var=", NumberForm[varY[[idx]], {8, 5}]
  ],
  {k, {1, 101, 201, 401}}
];

(* Example 3: AR(1) process in stationary regime *)
a = 0.8; sigmaW = 1.0; K = 400; mAR = 500;
sigmaX0 = sigmaW/Sqrt[1 - a^2];

arPaths = Table[
   Module[{x = ConstantArray[0., K]},
     x[[1]] = RandomVariate[NormalDistribution[0, sigmaX0]];
     Do[
       x[[k]] = a*x[[k - 1]] + RandomVariate[NormalDistribution[0, sigmaW]],
       {k, 2, K}
     ];
     x
   ],
   {mAR}
];

meanAR = Mean[arPaths];
varAR = Variance /@ Transpose[arPaths];

Print["AR(1) stationary regime estimates:"];
Do[
  idx = k;
  Print[
    "k=", idx - 1,
    ", mean=", NumberForm[meanAR[[idx]], {8, 5}],
    ", var=", NumberForm[varAR[[idx]], {8, 5}]
  ],
  {k, {1, 101, 201, 400}}
];
Print["Theoretical stationary variance = ", NumberForm[sigmaW^2/(1 - a^2), {8, 5}]];

(* Visualization *)
p1 = ListLinePlot[
  {Transpose[{tGrid, meanX}], Transpose[{tGrid, ConstantArray[0, Length[tGrid]]}]},
  PlotLegends -> {"Estimated mean", "Theory mean = 0"},
  PlotLabel -> "Random-phase harmonic process: ensemble mean",
  AxesLabel -> {"t", "Mean"},
  GridLines -> Automatic,
  ImageSize -> Large
];

p2 = ListLinePlot[
  Transpose[{tGrid, varY}],
  PlotLabel -> "Nonstationary process variance varies with time",
  AxesLabel -> {"t", "Var[Y(t)]"},
  GridLines -> Automatic,
  ImageSize -> Large
];

Column[{p1, p2}]

      

Library notes: Python uses NumPy and Matplotlib; C++ uses the standard library (<random>, <vector>, <cmath>); Java uses java.util.Random; MATLAB uses base functions and includes optional Simulink block creation; Mathematica uses built-in symbolic/numeric/statistical functions.

Problem 1 (SSS to WSS): Let \( X(t) \) be strictly stationary with finite second moments. Prove that \( m_X(t) \) is constant and \( C_X(t_1,t_2) \) depends only on lag.

Solution: The one-dimensional distributions of \( X(t) \) and \( X(t+\tau) \) are identical for every \( \tau \), so expectations are equal: \( \mathbb{E}[X(t+\tau)] = \mathbb{E}[X(t)] \). Thus the mean is constant. Likewise, the two-dimensional vectors \( (X(t_1),X(t_2)) \) and \( (X(t_1+\tau),X(t_2+\tau)) \) have the same distribution, so \( \mathbb{E}[X(t_1)X(t_2)] = \mathbb{E}[X(t_1+\tau)X(t_2+\tau)] \). Subtract the constant mean product to obtain \( C_X(t_1,t_2)=C_X(t_1+\tau,t_2+\tau) \). Therefore the process is WSS.

Problem 2 (Random-phase cosine autocorrelation): For \( X(t)=A\cos(\omega_0 t+\Theta) \) with \( \Theta \sim \operatorname{Unif}(0,2\pi) \), compute \( m_X(t) \) and \( R_X(t_1,t_2) \).

Solution: The mean is zero because the cosine integrates to zero over one full period:

\[ m_X(t)=\frac{A}{2\pi}\int_{0}^{2\pi}\cos(\omega_0 t+\theta)\,d\theta=0 \]

Using the trigonometric identity \( \cos u \cos v = \tfrac{1}{2}(\cos(u-v)+\cos(u+v)) \),

\[ R_X(t_1,t_2)=\mathbb{E}[X(t_1)X(t_2)] =\frac{A^2}{2}\cos\!\left(\omega_0(t_2-t_1)\right) \]

because the term involving \( \cos(u+v) \) averages to zero over the uniform phase. Hence the process is WSS.

Problem 3 (Stationarity test for affine process): Let \( Y(t)=\alpha t+\beta \) where \( \alpha,\beta \) are random constants. Determine conditions under which \( Y(t) \) is WSS.

Solution: WSS requires a constant mean:

\[ m_Y(t)=\mathbb{E}[\alpha]t+\mathbb{E}[\beta] \]

so we must have \( \mathbb{E}[\alpha]=0 \). Next, the variance must be constant:

\[ \operatorname{Var}(Y(t))=t^2\operatorname{Var}(\alpha)+2t\operatorname{Cov}(\alpha,\beta)+\operatorname{Var}(\beta) \]

For this to be independent of \( t \), the coefficients of \( t^2 \) and \( t \) must vanish. Therefore \( \operatorname{Var}(\alpha)=0 \) and \( \operatorname{Cov}(\alpha,\beta)=0 \). Since \( \operatorname{Var}(\alpha)=0 \), \( \alpha \) is almost surely constant, and combined with \( \mathbb{E}[\alpha]=0 \) this implies \( \alpha=0 \) almost surely. Thus the only WSS case is the degenerate constant-in-time process \( Y(t)=\beta \).

Problem 4 (AR(1) variance recursion): Consider \( x[k]=a x[k-1]+w[k] \) with zero-mean innovation \( w[k] \), variance \( \sigma_w^2 \), independent of \( x[k-1] \). Derive the variance recursion and the stationary variance when \( |a| < 1 \).

Solution: Since the process is zero-mean (if initialized with zero mean),

\[ \operatorname{Var}(x[k]) = \mathbb{E}[x[k]^2] = \mathbb{E}[(a x[k-1]+w[k])^2] \]

\[ \operatorname{Var}(x[k]) = a^2 \operatorname{Var}(x[k-1]) + \sigma_w^2 \]

because the cross term vanishes by independence and zero mean. In the stationary regime, \( \operatorname{Var}(x[k])=\operatorname{Var}(x[k-1])=\sigma_x^2 \), so

\[ \sigma_x^2 = a^2 \sigma_x^2 + \sigma_w^2 \quad \Rightarrow \quad \sigma_x^2 = \frac{\sigma_w^2}{1-a^2} \]

The requirement \( |a| < 1 \) ensures finite variance and convergence to the stationary distribution.

Problem 5 (Practical stationarity-check workflow): You are given ensemble simulation data for a process entering a mechanical vibration model. Outline a practical workflow to classify the process as nonstationary, WSS candidate, or stationary Gaussian.

Solution (workflow):

flowchart TD
  S["Start with ensemble or repeated runs"] --> M["Estimate mean versus time"]
  M --> Q1["Mean approximately constant?"]
  Q1 -->|"No"| NS["Classify as nonstationary"]
  Q1 -->|"Yes"| V["Estimate covariance at several time pairs"]
  V --> Q2["Covariance depends only on lag?"]
  Q2 -->|"No"| NS
  Q2 -->|"Yes"| WSS["WSS candidate"]
  WSS --> G["Check marginal/joint Gaussian assumption"]
  G --> Q3["Gaussian model acceptable?"]
  Q3 -->|"Yes"| SSS["Gaussian + WSS => strict stationary"]
  Q3 -->|"No"| KEEP["Keep as WSS unless higher-order \ntests are required"]
        

This workflow separates second-order verification (sufficient for many linear random vibration calculations) from stronger distribution-level verification needed for strict stationarity claims.

We introduced the probabilistic framework for stochastic system dynamics: random variables, random vectors, and random processes. We defined ensemble moments and covariance, formalized strict-sense and wide-sense stationarity, and proved key implications such as SSS \( \Rightarrow \) WSS (with finite second moments) and Gaussian WSS \( \Rightarrow \) SSS. The examples and code showed how stationarity appears in practice and how to verify it numerically from ensembles. These results prepare the ground for the next lesson on autocorrelation, power spectral density, and white-noise models.