In deterministic simulation, a dynamic model produces a single trajectory for a fixed set of parameters and initial conditions. In stochastic dynamics, at least one of the following is random: \( \boldsymbol{\theta} \) (parameters), \( \mathbf{x}_0 \) (initial state), or \( u(t) \) (input process). Therefore, the response \( \mathbf{x}(t,\omega) \) is itself a random process.

Monte Carlo simulation approximates statistical quantities of this response by repeatedly sampling the uncertainty source and integrating the dynamic model. A generic state model is

\[ \dot{\mathbf{x}}(t,\omega)=\mathbf{f}\!\big(t,\mathbf{x}(t,\omega),u(t,\omega),\boldsymbol{\theta}(\omega)\big), \qquad \mathbf{x}(0,\omega)=\mathbf{x}_0(\omega) \]

A quantity of interest (QoI) may be a scalar at a fixed time, e.g. \( Q(\omega)=x_1(T,\omega) \), a trajectory statistic \( Q(\omega)=\max_{0 \le t \le T}|x_1(t,\omega)| \), or a functional such as energy.

flowchart TD
  A["Define stochastic model: xdot = f(t,x,u,theta)"] --> B["Choose quantity of interest Q"]
  B --> C["Sample random variables / random input path"]
  C --> D["Numerically integrate one trajectory"]
  D --> E["Store Q and time history"]
  E --> F["Repeat for N trials"]
  F --> G["Estimate mean, variance, probability, CI"]
  G --> H["Check convergence and discretization error"]
        

Let \( x^{(r)}(t) \) denote the response from trial \( r=1,\dots,N \). For each fixed time \( t \), the sample mean estimator is \( \hat{\mu}_N(t) \):

\[ \hat{\mu}_N(t)=\frac{1}{N}\sum_{r=1}^{N}x^{(r)}(t) \]

The sample variance estimator is \( \hat{\sigma}_N^2(t) \):

\[ \hat{\sigma}_N^2(t)=\frac{1}{N-1}\sum_{r=1}^{N}\left(x^{(r)}(t)-\hat{\mu}_N(t)\right)^2 \]

For two response components or two time instants, the sample covariance estimator is \( \widehat{\operatorname{Cov}}_N \):

\[ \widehat{\operatorname{Cov}}_N\!\big(x_i(t_a),x_j(t_b)\big)= \frac{1}{N-1}\sum_{r=1}^{N} \Big(x_i^{(r)}(t_a)-\hat{\mu}_{i,N}(t_a)\Big)\Big(x_j^{(r)}(t_b)-\hat{\mu}_{j,N}(t_b)\Big) \]

For reliability-style questions, define an event \( A=\{\omega:\, Q(\omega)>q_{\mathrm{thr}}\} \). Then \( p=\mathbb{P}(A) \) is estimated by the Bernoulli sample mean:

\[ \hat{p}_N=\frac{1}{N}\sum_{r=1}^{N}\mathbf{1}_{A_r}, \qquad \mathbf{1}_{A_r}= \begin{cases} 1, & Q^{(r)}>q_{\mathrm{thr}} \\ 0, & Q^{(r)}\le q_{\mathrm{thr}} \end{cases} \]

These estimators are direct extensions of scalar Monte Carlo estimators, but now the random variable arises from a full numerical trajectory.

Assume the trials are independent and identically distributed (i.i.d.) and \( \mathbb{E}[|x(t)|]<\infty \). Then the sample mean is unbiased:

\[ \mathbb{E}[\hat{\mu}_N(t)] =\mathbb{E}\!\left[\frac{1}{N}\sum_{r=1}^{N}x^{(r)}(t)\right] =\frac{1}{N}\sum_{r=1}^{N}\mathbb{E}[x^{(r)}(t)] =\mu(t) \]

Under finite variance, \( \operatorname{Var}(x(t))=\sigma^2(t)<\infty \), the estimator variance is

\[ \operatorname{Var}\!\big(\hat{\mu}_N(t)\big)=\frac{\sigma^2(t)}{N} \]

Proof sketch: independence eliminates cross-covariance terms:

\[ \operatorname{Var}\!\left(\frac{1}{N}\sum_{r=1}^{N}x^{(r)}(t)\right) =\frac{1}{N^2}\sum_{r=1}^{N}\operatorname{Var}\!\big(x^{(r)}(t)\big) =\frac{N\sigma^2(t)}{N^2} =\frac{\sigma^2(t)}{N} \]

Hence the root-mean-square Monte Carlo sampling error decays as \( N^{-1/2} \), which is the characteristic Monte Carlo convergence rate.

By the Central Limit Theorem (CLT), for large \( N \):

\[ \frac{\hat{\mu}_N(t)-\mu(t)}{\hat{\sigma}_N(t)/\sqrt{N}} \approx \mathcal{N}(0,1) \]

Therefore, an approximate 95% confidence interval for \( \mu(t) \) is

\[ \hat{\mu}_N(t)\pm 1.96\,\frac{\hat{\sigma}_N(t)}{\sqrt{N}} \]

For exceedance probability, because \( \mathbf{1}_{A_r} \) is Bernoulli, \( \operatorname{Var}(\hat{p}_N)=p(1-p)/N \), and we use \( \hat{p}_N \) in place of \( p \) when constructing intervals.

In system dynamics, the response for each trial is usually computed numerically (Euler, RK2, RK4, etc.). Let \( Q_h(\omega) \) denote the QoI produced by a time-step \( h \). The practical estimator is then

\[ \hat{\mu}_{N,h}=\frac{1}{N}\sum_{r=1}^{N}Q_h^{(r)} \]

This estimator has two error sources:

1. Sampling error from finite \( N \).
2. Numerical discretization error from finite \( h \).

A useful decomposition of mean-square error (MSE) is

\[ \mathbb{E}\!\left[(\hat{\mu}_{N,h}-\mu)^2\right] = \Big(\mathbb{E}[Q_h]-\mathbb{E}[Q]\Big)^2 +\operatorname{Var}(\hat{\mu}_{N,h}) \]

and, under i.i.d. sampling,

\[ \operatorname{Var}(\hat{\mu}_{N,h})=\frac{\operatorname{Var}(Q_h)}{N} \]

If the numerical method is order \( p \) for the chosen QoI (e.g., global error \( O(h^p) \)), then the bias term is typically \( O(h^p) \), while the Monte Carlo standard error is \( O(N^{-1/2}) \). A sound simulation plan balances both terms rather than reducing only one.

flowchart TD
  A["Choose step size h"] --> B["Trajectory bias ~ O(h^p)"]
  C["Choose sample count N"] --> D["Sampling std ~ O(N^(-1/2))"]
  B --> E["Total estimator error"]
  D --> E
  E --> F["Refine h or increase N based on dominant term"]
        

Because Monte Carlo convergence is only \( O(N^{-1/2}) \), variance reduction is essential for expensive dynamic models.

5.1 Common Random Numbers (CRN). Suppose we compare two designs (A and B) using the same random samples in each trial. Let \( Q_A^{(r)} \) and \( Q_B^{(r)} \) be the two QoIs, and define \( D^{(r)}=Q_A^{(r)}-Q_B^{(r)} \). Then

\[ \operatorname{Var}(D)=\operatorname{Var}(Q_A)+\operatorname{Var}(Q_B)-2\operatorname{Cov}(Q_A,Q_B) \]

If the two outputs are positively correlated under shared randomness, the covariance term reduces the variance of the comparison estimator.

5.2 Antithetic Sampling. For symmetric inputs, pair each random vector \( \boldsymbol{\xi}^{(r)} \) with \( -\boldsymbol{\xi}^{(r)} \). The paired estimator

\[ \hat{\mu}_{\mathrm{anti}}=\frac{1}{N}\sum_{r=1}^{N/2}\frac{Q(\boldsymbol{\xi}^{(r)})+Q(-\boldsymbol{\xi}^{(r)})}{2} \]

can reduce variance when the QoI is monotone (or approximately monotone) in the driving randomness.

5.3 Control Variates. Let \( C \) be a correlated random variable with known mean \( \mathbb{E}[C] \). Define

\[ Q_{\mathrm{cv}}=Q-\beta\big(C-\mathbb{E}[C]\big) \]

Then \( \mathbb{E}[Q_{\mathrm{cv}}]=\mathbb{E}[Q] \), and the variance is minimized by

\[ \beta^\star=\frac{\operatorname{Cov}(Q,C)}{\operatorname{Var}(C)} \]

In linear stochastic dynamics, an analytically tractable approximation (from Lesson 3) is often an excellent control variate for a nonlinear or numerically integrated QoI.

We consider the second-order model

\[ m\ddot{x}(t)+c(\omega)\dot{x}(t)+k(\omega)x(t)=F(t,\omega) \]

with uncertain damping and stiffness, and sampled random forcing. In state form with \( x_1=x \), \( x_2=\dot{x} \):

\[ \dot{\mathbf{x}}= \begin{bmatrix} \dot{x}_1 \\ \dot{x}_2 \end{bmatrix} = \begin{bmatrix} x_2 \\ -\dfrac{k(\omega)}{m}x_1-\dfrac{c(\omega)}{m}x_2+\dfrac{1}{m}F(t,\omega) \end{bmatrix} \]

In the code examples below, we treat \( F(t,\omega) \) as piecewise-constant Gaussian samples over each time step. This is a practical teaching approximation for random forcing and is easy to implement consistently across programming languages.

The QoIs are:

• \( \mu_x(t)=\mathbb{E}[x(t)] \) and \( \sigma_x^2(t)=\operatorname{Var}(x(t)) \)
• \( p=\mathbb{P}(\max_{0 \le t \le T}|x(t)|>x_{\mathrm{thr}}) \)
• A 95% confidence interval for \( \mathbb{E}[x(T)] \)

Each trajectory is integrated by RK4, which you studied in Chapter 15. This is deliberate: the lesson focus is stochastic estimation, not introducing a new deterministic solver.

The script below performs \( N \) Monte Carlo runs, estimates time-varying mean/variance, computes an exceedance probability, and saves results to CSV. It also plots the mean with a \( \pm 2\sigma \) envelope.

Chapter17_Lesson5.py

# Chapter17_Lesson5.py
import numpy as np
import matplotlib.pyplot as plt

np.random.seed(17)

# ------------------------------------------------------------
# Monte Carlo simulation of a stochastic mass-spring-damper system
# m x'' + c x' + k x = F(t)
# Randomness:
#   - k, c (uncertain parameters)
#   - x(0), v(0) (random initial conditions)
#   - F(t_n) piecewise-constant Gaussian forcing
# ------------------------------------------------------------

T = 10.0
dt = 0.005
N = 2000
m = 1.0

# Nominal parameters and uncertainties
k_mean, k_std = 25.0, 2.0
c_mean, c_std = 1.5, 0.2

# Random initial conditions
x0_mean, x0_std = 0.0, 0.05
v0_mean, v0_std = 0.0, 0.05

# Random forcing (piecewise-constant samples)
sigma_F = 2.0

# Quantity of interest: exceedance probability of peak displacement
x_threshold = 0.75

t = np.arange(0.0, T + dt, dt)
nt = len(t)

sum_x = np.zeros(nt)
sum_x2 = np.zeros(nt)
peak_exceed_count = 0

def rk4_step(x, v, force, k, c, h):
    def f(state):
        xs, vs = state
        dx = vs
        dv = -(k / m) * xs - (c / m) * vs + force / m
        return np.array([dx, dv], dtype=float)

    y = np.array([x, v], dtype=float)
    k1 = f(y)
    k2 = f(y + 0.5 * h * k1)
    k3 = f(y + 0.5 * h * k2)
    k4 = f(y + h * k3)
    yn = y + (h / 6.0) * (k1 + 2*k2 + 2*k3 + k4)
    return yn[0], yn[1]

for trial in range(N):
    # Truncate to keep physical parameters positive
    k = max(1e-6, np.random.normal(k_mean, k_std))
    c = max(1e-6, np.random.normal(c_mean, c_std))

    x = np.random.normal(x0_mean, x0_std)
    v = np.random.normal(v0_mean, v0_std)

    x_hist = np.zeros(nt)
    x_hist[0] = x

    for n in range(nt - 1):
        F_n = np.random.normal(0.0, sigma_F)
        x, v = rk4_step(x, v, F_n, k, c, dt)
        x_hist[n + 1] = x

    sum_x += x_hist
    sum_x2 += x_hist * x_hist

    if np.max(np.abs(x_hist)) > x_threshold:
        peak_exceed_count += 1

mean_x = sum_x / N
var_x = (sum_x2 - N * mean_x**2) / (N - 1)
var_x = np.maximum(var_x, 0.0)
std_x = np.sqrt(var_x)

# 95% CI for the final-time mean displacement using CLT
z975 = 1.959963984540054
se_final = std_x[-1] / np.sqrt(N)
ci_final = (mean_x[-1] - z975 * se_final, mean_x[-1] + z975 * se_final)

p_hat = peak_exceed_count / N
se_p = np.sqrt(max(p_hat * (1.0 - p_hat), 1e-12) / N)
ci_p = (max(0.0, p_hat - z975 * se_p), min(1.0, p_hat + z975 * se_p))

print("Monte Carlo trajectories:", N)
print("Final mean x(T):", mean_x[-1])
print("Final variance x(T):", var_x[-1])
print("95% CI for E[x(T)]:", ci_final)
print("P(max |x| > %.3f): %.6f" % (x_threshold, p_hat))
print("95% CI for probability:", ci_p)

# Save time history statistics for post-processing
data = np.column_stack([t, mean_x, var_x, std_x])
np.savetxt(
    "Chapter17_Lesson5_python_results.csv",
    data,
    delimiter=",",
    header="t,mean_x,var_x,std_x",
    comments=""
)

# Plot mean and ±2 std envelope
plt.figure(figsize=(8, 4.5))
plt.plot(t, mean_x, label="Mean displacement")
plt.plot(t, mean_x + 2*std_x, "--", label="Mean + 2 std")
plt.plot(t, mean_x - 2*std_x, "--", label="Mean - 2 std")
plt.xlabel("Time (s)")
plt.ylabel("Displacement x(t)")
plt.title("Monte Carlo response statistics")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("Chapter17_Lesson5_python_plot.png", dpi=160)
plt.show()

This C++ implementation uses the standard library (<random>, <vector>) and computes the same estimators as the Python version. The output CSV can be plotted using your preferred tool.

Chapter17_Lesson5.cpp

// Chapter17_Lesson5.cpp
#include <iostream>
#include <vector>
#include <random>
#include <fstream>
#include <cmath>
#include <algorithm>
#include <iomanip>

struct State {
    double x;
    double v;
};

State dynamics(const State& s, double force, double k, double c, double m) {
    State ds;
    ds.x = s.v;
    ds.v = -(k / m) * s.x - (c / m) * s.v + force / m;
    return ds;
}

State rk4Step(const State& s, double force, double k, double c, double m, double h) {
    State k1 = dynamics(s, force, k, c, m);

    State s2{ s.x + 0.5 * h * k1.x, s.v + 0.5 * h * k1.v };
    State k2 = dynamics(s2, force, k, c, m);

    State s3{ s.x + 0.5 * h * k2.x, s.v + 0.5 * h * k2.v };
    State k3 = dynamics(s3, force, k, c, m);

    State s4{ s.x + h * k3.x, s.v + h * k3.v };
    State k4 = dynamics(s4, force, k, c, m);

    State out;
    out.x = s.x + (h / 6.0) * (k1.x + 2.0 * k2.x + 2.0 * k3.x + k4.x);
    out.v = s.v + (h / 6.0) * (k1.v + 2.0 * k2.v + 2.0 * k3.v + k4.v);
    return out;
}

int main() {
    const double T = 10.0;
    const double dt = 0.005;
    const int N = 2000;
    const double m = 1.0;

    const double kMean = 25.0, kStd = 2.0;
    const double cMean = 1.5, cStd = 0.2;
    const double x0Mean = 0.0, x0Std = 0.05;
    const double v0Mean = 0.0, v0Std = 0.05;
    const double sigmaF = 2.0;
    const double xThreshold = 0.75;

    const int nt = static_cast<int>(std::round(T / dt)) + 1;

    std::vector<double> t(nt), sumX(nt, 0.0), sumX2(nt, 0.0);
    for (int i = 0; i < nt; ++i) t[i] = i * dt;

    std::mt19937 rng(17);
    std::normal_distribution<double> distK(kMean, kStd);
    std::normal_distribution<double> distC(cMean, cStd);
    std::normal_distribution<double> distX0(x0Mean, x0Std);
    std::normal_distribution<double> distV0(v0Mean, v0Std);
    std::normal_distribution<double> distF(0.0, sigmaF);

    int exceedCount = 0;

    for (int trial = 0; trial < N; ++trial) {
        double k = std::max(1e-6, distK(rng));
        double c = std::max(1e-6, distC(rng));

        State s{ distX0(rng), distV0(rng) };

        std::vector<double> xHist(nt, 0.0);
        xHist[0] = s.x;

        double peakAbs = std::abs(s.x);

        for (int n = 0; n < nt - 1; ++n) {
            double F = distF(rng);
            s = rk4Step(s, F, k, c, m, dt);
            xHist[n + 1] = s.x;
            peakAbs = std::max(peakAbs, std::abs(s.x));
        }

        for (int n = 0; n < nt; ++n) {
            sumX[n] += xHist[n];
            sumX2[n] += xHist[n] * xHist[n];
        }

        if (peakAbs > xThreshold) {
            exceedCount++;
        }
    }

    std::vector<double> meanX(nt), varX(nt), stdX(nt);
    for (int n = 0; n < nt; ++n) {
        meanX[n] = sumX[n] / N;
        varX[n] = (sumX2[n] - N * meanX[n] * meanX[n]) / (N - 1.0);
        if (varX[n] < 0.0) varX[n] = 0.0;
        stdX[n] = std::sqrt(varX[n]);
    }

    const double z975 = 1.959963984540054;
    double seFinal = stdX.back() / std::sqrt(static_cast<double>(N));
    double ciMeanLo = meanX.back() - z975 * seFinal;
    double ciMeanHi = meanX.back() + z975 * seFinal;

    double pHat = static_cast<double>(exceedCount) / N;
    double seP = std::sqrt(std::max(pHat * (1.0 - pHat), 1e-12) / N);
    double ciPLo = std::max(0.0, pHat - z975 * seP);
    double ciPHi = std::min(1.0, pHat + z975 * seP);

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "Monte Carlo trajectories: " << N << "\n";
    std::cout << "Final mean x(T): " << meanX.back() << "\n";
    std::cout << "Final variance x(T): " << varX.back() << "\n";
    std::cout << "95% CI for E[x(T)]: [" << ciMeanLo << ", " << ciMeanHi << "]\n";
    std::cout << "P(max |x| > " << xThreshold << "): " << pHat << "\n";
    std::cout << "95% CI for probability: [" << ciPLo << ", " << ciPHi << "]\n";

    std::ofstream csv("Chapter17_Lesson5_cpp_results.csv");
    csv << "t,mean_x,var_x,std_x\n";
    for (int n = 0; n < nt; ++n) {
        csv << t[n] << "," << meanX[n] << "," << varX[n] << "," << stdX[n] << "\n";
    }
    csv.close();

    return 0;
}

The Java code mirrors the same RK4 + Monte Carlo structure and writes the statistics to a CSV file. This is useful for students who build system dynamics tools in Java-based engineering applications.

Chapter17_Lesson5.java

// Chapter17_Lesson5.java
import java.io.PrintWriter;
import java.io.IOException;
import java.util.Random;

public class Chapter17_Lesson5 {
    static class State {
        double x;
        double v;
        State(double x, double v) { this.x = x; this.v = v; }
    }

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

    static State dynamics(State s, double force, double k, double c, double m) {
        double dx = s.v;
        double dv = -(k / m) * s.x - (c / m) * s.v + force / m;
        return new State(dx, dv);
    }

    static State rk4Step(State s, double force, double k, double c, double m, double h) {
        State k1 = dynamics(s, force, k, c, m);

        State s2 = new State(s.x + 0.5 * h * k1.x, s.v + 0.5 * h * k1.v);
        State k2 = dynamics(s2, force, k, c, m);

        State s3 = new State(s.x + 0.5 * h * k2.x, s.v + 0.5 * h * k2.v);
        State k3 = dynamics(s3, force, k, c, m);

        State s4 = new State(s.x + h * k3.x, s.v + h * k3.v);
        State k4 = dynamics(s4, force, k, c, m);

        double xn = s.x + (h / 6.0) * (k1.x + 2.0 * k2.x + 2.0 * k3.x + k4.x);
        double vn = s.v + (h / 6.0) * (k1.v + 2.0 * k2.v + 2.0 * k3.v + k4.v);
        return new State(xn, vn);
    }

    public static void main(String[] args) throws IOException {
        final double T = 10.0;
        final double dt = 0.005;
        final int N = 2000;
        final double m = 1.0;

        final double kMean = 25.0, kStd = 2.0;
        final double cMean = 1.5, cStd = 0.2;
        final double x0Mean = 0.0, x0Std = 0.05;
        final double v0Mean = 0.0, v0Std = 0.05;
        final double sigmaF = 2.0;
        final double xThreshold = 0.75;

        final int nt = (int)Math.round(T / dt) + 1;

        double[] t = new double[nt];
        double[] sumX = new double[nt];
        double[] sumX2 = new double[nt];
        for (int i = 0; i < nt; i++) {
            t[i] = i * dt;
        }

        Random rng = new Random(17L);
        int exceedCount = 0;

        for (int trial = 0; trial < N; trial++) {
            double k = Math.max(1e-6, gaussian(rng, kMean, kStd));
            double c = Math.max(1e-6, gaussian(rng, cMean, cStd));

            State s = new State(gaussian(rng, x0Mean, x0Std), gaussian(rng, v0Mean, v0Std));

            double[] xHist = new double[nt];
            xHist[0] = s.x;

            double peakAbs = Math.abs(s.x);

            for (int n = 0; n < nt - 1; n++) {
                double F = gaussian(rng, 0.0, sigmaF);
                s = rk4Step(s, F, k, c, m, dt);
                xHist[n + 1] = s.x;
                peakAbs = Math.max(peakAbs, Math.abs(s.x));
            }

            for (int n = 0; n < nt; n++) {
                sumX[n] += xHist[n];
                sumX2[n] += xHist[n] * xHist[n];
            }

            if (peakAbs > xThreshold) exceedCount++;
        }

        double[] meanX = new double[nt];
        double[] varX = new double[nt];
        double[] stdX = new double[nt];

        for (int n = 0; n < nt; n++) {
            meanX[n] = sumX[n] / N;
            varX[n] = (sumX2[n] - N * meanX[n] * meanX[n]) / (N - 1.0);
            if (varX[n] < 0.0) varX[n] = 0.0;
            stdX[n] = Math.sqrt(varX[n]);
        }

        final double z975 = 1.959963984540054;
        double seFinal = stdX[nt - 1] / Math.sqrt((double)N);
        double ciMeanLo = meanX[nt - 1] - z975 * seFinal;
        double ciMeanHi = meanX[nt - 1] + z975 * seFinal;

        double pHat = ((double)exceedCount) / N;
        double seP = Math.sqrt(Math.max(pHat * (1.0 - pHat), 1e-12) / N);
        double ciPLo = Math.max(0.0, pHat - z975 * seP);
        double ciPHi = Math.min(1.0, pHat + z975 * seP);

        System.out.printf("Monte Carlo trajectories: %d%n", N);
        System.out.printf("Final mean x(T): %.6f%n", meanX[nt - 1]);
        System.out.printf("Final variance x(T): %.6f%n", varX[nt - 1]);
        System.out.printf("95%% CI for E[x(T)]: [%.6f, %.6f]%n", ciMeanLo, ciMeanHi);
        System.out.printf("P(max |x| > %.3f): %.6f%n", xThreshold, pHat);
        System.out.printf("95%% CI for probability: [%.6f, %.6f]%n", ciPLo, ciPHi);

        try (PrintWriter pw = new PrintWriter("Chapter17_Lesson5_java_results.csv")) {
            pw.println("t,mean_x,var_x,std_x");
            for (int n = 0; n < nt; n++) {
                pw.printf(java.util.Locale.US, "%.6f,%.10f,%.10f,%.10f%n",
                        t[n], meanX[n], varX[n], stdX[n]);
            }
        }
    }
}

The MATLAB script performs Monte Carlo simulation directly and includes a note on the corresponding block-level organization in Simulink (state integrators, random source, and repeated simulation wrapper).

Chapter17_Lesson5.m

% Chapter17_Lesson5.m
% Monte Carlo simulation of a stochastic mass-spring-damper system
% m*xdd + c*xd + k*x = F(t)
% Randomness in k, c, initial conditions, and sampled forcing.
clear; clc; close all;

rng(17);

T = 10.0;
dt = 0.005;
N = 2000;
m = 1.0;

kMean = 25.0; kStd = 2.0;
cMean = 1.5;  cStd = 0.2;
x0Mean = 0.0; x0Std = 0.05;
v0Mean = 0.0; v0Std = 0.05;
sigmaF = 2.0;
xThreshold = 0.75;

t = 0:dt:T;
nt = numel(t);

sumX = zeros(1, nt);
sumX2 = zeros(1, nt);
exceedCount = 0;

for trial = 1:N
    k = max(1e-6, kMean + kStd * randn());
    c = max(1e-6, cMean + cStd * randn());

    x = x0Mean + x0Std * randn();
    v = v0Mean + v0Std * randn();

    xHist = zeros(1, nt);
    xHist(1) = x;

    peakAbs = abs(x);

    for n = 1:(nt - 1)
        F = sigmaF * randn();
        [x, v] = rk4Step(x, v, F, k, c, m, dt);
        xHist(n + 1) = x;
        peakAbs = max(peakAbs, abs(x));
    end

    sumX = sumX + xHist;
    sumX2 = sumX2 + xHist.^2;

    if peakAbs > xThreshold
        exceedCount = exceedCount + 1;
    end
end

meanX = sumX / N;
varX = (sumX2 - N * meanX.^2) / (N - 1);
varX = max(varX, 0);
stdX = sqrt(varX);

z975 = 1.959963984540054;
seFinal = stdX(end) / sqrt(N);
ciMean = [meanX(end) - z975 * seFinal, meanX(end) + z975 * seFinal];

pHat = exceedCount / N;
seP = sqrt(max(pHat * (1 - pHat), 1e-12) / N);
ciP = [max(0, pHat - z975 * seP), min(1, pHat + z975 * seP)];

fprintf('Monte Carlo trajectories: %d\n', N);
fprintf('Final mean x(T): %.6f\n', meanX(end));
fprintf('Final variance x(T): %.6f\n', varX(end));
fprintf('95%% CI for E[x(T)]: [%.6f, %.6f]\n', ciMean(1), ciMean(2));
fprintf('P(max |x| > %.3f): %.6f\n', xThreshold, pHat);
fprintf('95%% CI for probability: [%.6f, %.6f]\n', ciP(1), ciP(2));

resultTable = table(t(:), meanX(:), varX(:), stdX(:), ...
    'VariableNames', {'t','mean_x','var_x','std_x'});
writetable(resultTable, 'Chapter17_Lesson5_matlab_results.csv');

figure;
plot(t, meanX, 'LineWidth', 1.5); hold on;
plot(t, meanX + 2 * stdX, '--', 'LineWidth', 1.0);
plot(t, meanX - 2 * stdX, '--', 'LineWidth', 1.0);
grid on;
xlabel('Time (s)');
ylabel('Displacement x(t)');
title('Monte Carlo response statistics');
legend('Mean', 'Mean + 2 std', 'Mean - 2 std', 'Location', 'best');

% ---------------- Simulink note ----------------
% Equivalent Simulink model:
% 1) State-Space block or two Integrator blocks for x and v
% 2) Random Number block (or Band-Limited White Noise) for F(t)
% 3) MATLAB Function block to sample random k,c per run (or use parsim)
% 4) To Workspace blocks to collect x(t)
% 5) Monte Carlo wrapper script using repeated sim(...) calls

function [xNext, vNext] = rk4Step(x, v, force, k, c, m, h)
    f = @(state) [state(2);
                  -(k / m) * state(1) - (c / m) * state(2) + force / m];

    y = [x; v];
    k1 = f(y);
    k2 = f(y + 0.5 * h * k1);
    k3 = f(y + 0.5 * h * k2);
    k4 = f(y + h * k3);

    yNext = y + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4);
    xNext = yNext(1);
    vNext = yNext(2);
end

The following Mathematica code performs the same Monte Carlo experiment, exports a CSV file, and plots the mean and \( \pm 2\sigma \) envelopes.

Chapter17_Lesson5.nb

(* Chapter17_Lesson5.nb *)
(* Monte Carlo simulation of a stochastic mass-spring-damper system *)
SeedRandom[17];

T = 10.0;
dt = 0.005;
nSamples = 2000;
m = 1.0;

kMean = 25.0; kStd = 2.0;
cMean = 1.5; cStd = 0.2;
x0Mean = 0.0; x0Std = 0.05;
v0Mean = 0.0; v0Std = 0.05;
sigmaF = 2.0;
xThreshold = 0.75;

tGrid = N[Range[0, T, dt]];
nt = Length[tGrid];

rk4Step[{x_, v_}, force_, k_, c_, m_, h_] := Module[
  {f, y, k1, k2, k3, k4, yNext},
  f[{xs_, vs_}] := {vs, -(k/m) xs - (c/m) vs + force/m};
  y = {x, v};
  k1 = f[y];
  k2 = f[y + 0.5 h k1];
  k3 = f[y + 0.5 h k2];
  k4 = f[y + h k3];
  yNext = y + (h/6.0) (k1 + 2 k2 + 2 k3 + k4);
  yNext
];

simulateOne[] := Module[
  {k, c, x, v, xHist, force, n},
  k = Max[10^-6, RandomVariate[NormalDistribution[kMean, kStd]]];
  c = Max[10^-6, RandomVariate[NormalDistribution[cMean, cStd]]];
  x = RandomVariate[NormalDistribution[x0Mean, x0Std]];
  v = RandomVariate[NormalDistribution[v0Mean, v0Std]];
  xHist = ConstantArray[0.0, nt];
  xHist[[1]] = x;
  Do[
    force = RandomVariate[NormalDistribution[0, sigmaF]];
    {x, v} = rk4Step[{x, v}, force, k, c, m, dt];
    xHist[[n + 1]] = x;
    ,
    {n, 1, nt - 1}
  ];
  xHist
];

trajectories = Table[simulateOne[], {nSamples}];

meanX = Mean[trajectories];
varX = Variance /@ Transpose[trajectories];
stdX = Sqrt[varX];

z975 = 1.959963984540054;
seFinal = stdX[[-1]]/Sqrt[nSamples];
ciMean = {meanX[[-1]] - z975 seFinal, meanX[[-1]] + z975 seFinal};

pHat = N[Count[trajectories, tr_ /; Max[Abs[tr]] > xThreshold]/nSamples];
seP = Sqrt[Max[pHat (1 - pHat), 10^-12]/nSamples];
ciP = {Max[0, pHat - z975 seP], Min[1, pHat + z975 seP]};

Print["Monte Carlo trajectories: ", nSamples];
Print["Final mean x(T): ", NumberForm[meanX[[-1]], {10, 6}]];
Print["Final variance x(T): ", NumberForm[varX[[-1]], {10, 6}]];
Print["95% CI for E[x(T)]: ", ciMean];
Print["P(max |x| > ", xThreshold, "): ", NumberForm[pHat, {10, 6}]];
Print["95% CI for probability: ", ciP];

result = Transpose[{tGrid, meanX, varX, stdX}];
Export["Chapter17_Lesson5_mathematica_results.csv",
  Prepend[result, {"t", "mean_x", "var_x", "std_x"}]
];

ListLinePlot[
  {
    Transpose[{tGrid, meanX}],
    Transpose[{tGrid, meanX + 2 stdX}],
    Transpose[{tGrid, meanX - 2 stdX}]
  },
  PlotLegends -> {"Mean", "Mean + 2 std", "Mean - 2 std"},
  GridLines -> Automatic,
  AxesLabel -> {"Time (s)", "Displacement x(t)"},
  PlotLabel -> "Monte Carlo response statistics"
]

A reliable Monte Carlo study of a dynamic system should report both stochastic and numerical accuracy. Recommended checks:

• Repeat with \( N \) and \( 2N \) to verify confidence interval shrinkage by roughly \( 1/\sqrt{2} \).
• Repeat with \( h \) and \( h/2 \) to ensure the mean trajectory and QoI estimates are numerically converged.
• Compare against analytic mean/variance results (Lesson 3) for a simplified linear case whenever possible.
• Use fixed random seeds for reproducibility, but also verify that conclusions are stable under a different seed.

For expensive models, parallelization across trials is natural because trajectories are independent. However, reproducible parallel Monte Carlo requires careful random number stream management (independent substreams).

Problem 1 (Unbiasedness of the Monte Carlo Mean): Let \( Q^{(1)},\dots,Q^{(N)} \) be i.i.d. QoI samples with \( \mathbb{E}[Q]=\mu \). Prove that the estimator \( \hat{\mu}_N=\frac{1}{N}\sum_{r=1}^{N}Q^{(r)} \) is unbiased, and derive its variance.

Solution:

By linearity of expectation,

\[ \mathbb{E}[\hat{\mu}_N] =\mathbb{E}\!\left[\frac{1}{N}\sum_{r=1}^{N}Q^{(r)}\right] =\frac{1}{N}\sum_{r=1}^{N}\mathbb{E}[Q^{(r)}] =\frac{1}{N}\sum_{r=1}^{N}\mu =\mu \]

so the estimator is unbiased. For the variance, use independence:

\[ \operatorname{Var}(\hat{\mu}_N) =\operatorname{Var}\!\left(\frac{1}{N}\sum_{r=1}^{N}Q^{(r)}\right) =\frac{1}{N^2}\sum_{r=1}^{N}\operatorname{Var}(Q^{(r)}) =\frac{\sigma_Q^2}{N} \]

where \( \sigma_Q^2=\operatorname{Var}(Q) \). Hence the standard error is \( \sigma_Q/\sqrt{N} \).

Problem 2 (Confidence Interval Width and Sample Size): Suppose the sample standard deviation of \( x(T) \) is \( \hat{\sigma}=0.30 \). How many trajectories are needed so that the 95% confidence interval half-width for \( \mathbb{E}[x(T)] \) is at most \( 0.01 \)?

Solution:

Using the normal approximation, the half-width is \( h_{\mathrm{CI}}=1.96\,\hat{\sigma}/\sqrt{N} \). Set \( h_{\mathrm{CI}}\le 0.01 \):

\[ 1.96\frac{0.30}{\sqrt{N}}\le 0.01 \quad \Longrightarrow \quad \sqrt{N}\ge 58.8 \quad \Longrightarrow \quad N\ge 3457.44 \]

Therefore choose \( N=3458 \) trajectories (or a slightly larger rounded value such as \( N=3500 \)).

Problem 3 (Bias–Variance Decomposition with Numerical Time Step): Let \( Q_h \) be the numerical QoI computed with step size \( h \) and \( \mu=\mathbb{E}[Q] \) be the exact QoI mean. Show that

\[ \mathbb{E}\!\left[(\hat{\mu}_{N,h}-\mu)^2\right] = \Big(\mathbb{E}[Q_h]-\mu\Big)^2+\frac{\operatorname{Var}(Q_h)}{N} \]

Solution:

Add and subtract \( \mathbb{E}[Q_h] \):

\[ \hat{\mu}_{N,h}-\mu = \big(\hat{\mu}_{N,h}-\mathbb{E}[Q_h]\big) +\big(\mathbb{E}[Q_h]-\mu\big) \]

Square both sides and take expectations. The cross term vanishes because \( \mathbb{E}[\hat{\mu}_{N,h}-\mathbb{E}[Q_h]]=0 \). Thus:

\[ \mathbb{E}\!\left[(\hat{\mu}_{N,h}-\mu)^2\right] = \operatorname{Var}(\hat{\mu}_{N,h}) +\big(\mathbb{E}[Q_h]-\mu\big)^2 \]

and by i.i.d. sampling of \( Q_h^{(r)} \),

\[ \operatorname{Var}(\hat{\mu}_{N,h})=\frac{\operatorname{Var}(Q_h)}{N} \]

which proves the expression.

Problem 4 (Exceedance Probability Estimator): Let \( A_r=\{\max_{0 \le t \le T}|x^{(r)}(t)|>x_{\mathrm{thr}}\} \). Show that the estimator \( \hat{p}_N=\frac{1}{N}\sum_{r=1}^{N}\mathbf{1}_{A_r} \) is unbiased for \( p=\mathbb{P}(A) \), and derive an estimator for its standard error.

Solution:

Each indicator \( \mathbf{1}_{A_r} \) is Bernoulli with success probability \( p \), so

\[ \mathbb{E}[\mathbf{1}_{A_r}] = p, \qquad \operatorname{Var}(\mathbf{1}_{A_r}) = p(1-p) \]

Therefore

\[ \mathbb{E}[\hat{p}_N]=p, \qquad \operatorname{Var}(\hat{p}_N)=\frac{p(1-p)}{N} \]

Since \( p \) is unknown, we estimate the standard error using

\[ \widehat{\operatorname{SE}}(\hat{p}_N)= \sqrt{\frac{\hat{p}_N(1-\hat{p}_N)}{N}} \]

and use the normal approximation interval \( \hat{p}_N \pm 1.96\,\widehat{\operatorname{SE}}(\hat{p}_N) \) when \( N \) is sufficiently large and \( p \) is not extremely close to 0 or 1.

Problem 5 (Optimal Control Variate Coefficient): For \( Q_{\mathrm{cv}}=Q-\beta(C-\mathbb{E}[C]) \), derive the value of \( \beta \) that minimizes \( \operatorname{Var}(Q_{\mathrm{cv}}) \).

Solution:

Expand the variance:

\[ \operatorname{Var}(Q_{\mathrm{cv}}) = \operatorname{Var}(Q)+\beta^2\operatorname{Var}(C)-2\beta\operatorname{Cov}(Q,C) \]

Differentiate with respect to \( \beta \) and set to zero:

\[ \frac{d}{d\beta}\operatorname{Var}(Q_{\mathrm{cv}}) = 2\beta\operatorname{Var}(C)-2\operatorname{Cov}(Q,C)=0 \]

Hence

\[ \beta^\star=\frac{\operatorname{Cov}(Q,C)}{\operatorname{Var}(C)} \]

Substituting \( \beta^\star \) gives the minimum variance

\[ \operatorname{Var}(Q_{\mathrm{cv}}^\star) = \operatorname{Var}(Q)\big(1-\rho_{QC}^2\big) \]

where \( \rho_{QC} \) is the correlation coefficient. The stronger the correlation, the larger the variance reduction.

Monte Carlo simulation converts uncertainty in parameters, inputs, and initial conditions into statistical estimates of dynamic response quantities. The key mathematical facts are: unbiased sample estimators, \( O(N^{-1/2}) \) sampling convergence, CLT-based confidence intervals, and a clear separation between sampling error and numerical integration bias. In practice, effective studies combine a stable time integrator (e.g., RK4), convergence checks in both \( N \) and \( h \), and variance reduction when trajectory simulation is expensive.