Example: Global Sensitivity Analysis of the Duffing Oscillator

Example: Global Sensitivity Analysis of the Duffing Oscillator#

The preceding theory decomposed the variance of a scalar model output into main effects and interactions. This section applies that decomposition to the Duffing oscillator. A scrambled Sobol sequence supplies the base points, while the column-swapped pick–freeze construction supplies the model evaluations needed to estimate the indices (Saltelli, 2002; Saltelli et al., 2010). We use SALib (Herman and Usher, 2017; Iwanaga et al., 2022), which reports the estimated first-order indices \(\widehat{S}_i\) as S1 and the estimated total-effect indices \(\widehat{S}_{T_i}\) as ST.

Hide code cell source

 
from importlib.metadata import PackageNotFoundError, version
import shutil
import subprocess
import sys

required_versions = {
    "jax": "0.11.1",
    "jaxlib": "0.11.1",
    "diffrax": "0.7.2",
    "SALib": "1.5.2",
}
try:
    needs_install = any(
        version(package) != required_version
        for package, required_version in required_versions.items()
    )
except PackageNotFoundError:
    needs_install = True

if needs_install:
    requirements = [
        f"{package}=={required_version}"
        for package, required_version in required_versions.items()
    ]
    if shutil.which("uv"):
        command = ["uv", "pip", "install", "--python", sys.executable, *requirements]
    else:
        command = [sys.executable, "-m", "pip", "install", *requirements]
    subprocess.check_call(command)

import numpy as np
from jax import vmap
import jax.numpy as jnp
from diffrax import diffeqsolve, ODETerm, SaveAt, Tsit5
from SALib.sample import sobol as sample_sobol
import SALib.analyze.sobol as analyze_sobol

Duffing oscillator#

Let \(t\in[0,50]\) denote time, \(x(t)\) displacement, and \(v(t)=\dot{x}(t)\) velocity. The parameter vector is \(\boldsymbol{\theta}=(\alpha,\beta,\gamma,\delta,\omega)\), where \(\alpha\) and \(\beta\) are the linear and cubic stiffness coefficients, \(\gamma\) is the forcing amplitude, \(\delta\) is the damping coefficient, and \(\omega\) is the forcing angular frequency. The Duffing oscillator satisfies

\[\begin{split} \begin{aligned} \dot{x}(t) &= v(t),\\ \dot{v}(t) &= -\delta v(t)-\alpha x(t)-\beta x(t)^3+\gamma\cos(\omega t),\\ x(0)&=v(0)=0. \end{aligned} \end{split}\]

The state vector is \(\mathbf{y}(t)=(x(t),v(t))^{\mathsf T}\). The nominal parameter vector, denoted \(\boldsymbol{\mu}\) in the local-sensitivity example, is

\[ \boldsymbol{\theta}_0=(1,5,0.37,0.1,1). \]
def vector_field(t, y, theta):
    alpha, beta, gamma, delta, omega = theta
    x = y[0]
    v = y[1]
    return jnp.array(
        [
            v,
            - alpha * x - beta * x ** 3 - delta * v + gamma * jnp.cos(omega * t)
        ]
    )

theta_nominal = jnp.array([
    1.0,  # alpha
    5.0,  # beta
    0.37, # gamma
    0.1,  # delta
    1.0,  # omega
])
solver = Tsit5()
save_times = jnp.linspace(0.0, 50.0, 500)
saveat = SaveAt(ts=save_times)
term = ODETerm(vector_field)

nominal_solution = diffeqsolve(
    term,
    solver,
    t0=0,                       # Initial time
    t1=50,                      # Terminal time
    dt0=0.1,                    # Fixed internal step size
    y0=jnp.array([0.0, 0.0]),
    args=theta_nominal,
    saveat=saveat
)

fig, ax = plt.subplots(figsize=FIGURE_SIZES["full_standard"])
markevery = max(1, len(nominal_solution.ts) // 12)
ax.plot(nominal_solution.ts, nominal_solution.ys[:, 0], color="0.10", linestyle="-", marker="o", markevery=markevery, markersize=3, label="x")
ax.plot(nominal_solution.ts, nominal_solution.ys[:, 1], color="0.45", linestyle="--", marker="s", markevery=markevery, markersize=3, label="v")
ax.set(xlabel="t", ylabel="x(t), v(t)")
ax.title.set_text("Nominal trajectory")
ax.legend(frameon=False)
finalize_axes(keep_box=False);

fig, ax = plt.subplots(figsize=FIGURE_SIZES["full_standard"])
ax.plot(nominal_solution.ys[:, 0], nominal_solution.ys[:, 1], color="0.15", lw=1)
ax.set(xlabel="x", ylabel="v")
ax.title.set_text("Nominal phase portrait")
finalize_axes(keep_box=False);
Nominal Duffing displacement and velocity versus time together with the corresponding displacement-velocity phase portrait. Nominal Duffing displacement and velocity versus time together with the corresponding displacement-velocity phase portrait.

Input distribution#

Sobol indices depend on the assigned input distribution. Let \(\boldsymbol{\Theta}=(\Theta_1,\ldots,\Theta_5)\) denote the random parameter vector in the same parameter order as \(\boldsymbol{\theta}_0\). We model its five components as mutually independent, with

\[ \Theta_i\sim\mathcal{U}\left[0.99(\boldsymbol{\theta}_0)_i,\,1.01(\boldsymbol{\theta}_0)_i\right],\qquad i=1,\ldots,5. \]

Here \(\mathcal{U}[a,b]\) denotes the uniform distribution on \([a,b]\). Each interval has a relative half-width of \(1\%\). At a fixed time \(t>0\), the scalar outputs are \(Y_x(t)=x(t;\boldsymbol{\Theta})\) and \(Y_v(t)=v(t;\boldsymbol{\Theta})\). At \(t=0\), both outputs are fixed at zero, so their variances vanish and their Sobol indices are undefined. The resulting indices are global relative to this product distribution; changing the ranges or distributions can change the sensitivity ranking.

relative_half_width = 0.01
parameter_names = ["alpha", "beta", "gamma", "delta", "omega"]
parameter_labels = [r"$\alpha$", r"$\beta$", r"$\gamma$", r"$\delta$", r"$\omega$"]
theta_nominal_np = np.asarray(theta_nominal)

bounds = np.column_stack(
    (
        (1.0 - relative_half_width) * theta_nominal_np,
        (1.0 + relative_half_width) * theta_nominal_np,
    )
).tolist()

problem = {
    "num_vars": len(parameter_names),
    "names": parameter_names,
    "bounds": bounds,
}

Pick–freeze design#

Let \(N\) be the base sample size and let \(d=5\) be the number of inputs. With pairwise second-order indices disabled, SALib expands the Saltelli cross-sampled design into \(N(d+2)\) parameter vectors: two base vectors and \(d\) column-swapped hybrid vectors for each base point. This one design estimates both S1 and ST; estimating the pairwise interaction indices \(S_{ij}\), reported as S2, would require the larger \(N(2d+2)\) design (Saltelli, 2002; Saltelli et al., 2010). We use a power-of-two base size and a fixed seed for the scrambled Sobol points.

N = 512
d = problem["num_vars"]
sampling_seed = 1729
param_values = sample_sobol.sample(
    problem,
    N,
    calc_second_order=False,
    scramble=True,
    seed=sampling_seed,
)
assert param_values.shape == (N * (d + 2), d)
print(f"{param_values.shape[0]} parameter vectors with {param_values.shape[1]} parameters each")

3584 parameter vectors with 5 parameters each

The design contains \(512(5+2)=3584\) rows and five columns. Each row is one parameter vector at which the model is evaluated, and the full design targets the product distribution of \(\boldsymbol{\Theta}\). The columns follow the parameter order \((\alpha,\beta,\gamma,\delta,\omega)\). One model evaluation maps each row to complete displacement and velocity trajectories.

def run_model(parameter_vector):
    solution = diffeqsolve(
        term,
        solver,
        t0=0,
        t1=50,
        dt0=0.1,  # Fixed internal step size
        y0=jnp.array([0.0, 0.0]),
        args=parameter_vector,
        saveat=saveat,
    )
    return solution.ys


Y = np.asarray(vmap(run_model)(jnp.asarray(param_values)))

Y_x, Y_v = Y[:, :, 0], Y[:, :, 1]

Time-dependent indices#

For each saved time \(t>0\), we analyze the displacement and velocity outputs separately. A single call to SALib returns both S1 and ST. We use 100 bootstrap resamples to calculate 95% pointwise confidence half-widths; the analysis seed is separate from the seed that generated the scrambled design. We plot only the point estimates here. Very small output variance near \(t=0\) can make any normalized sensitivity estimate unstable, even after the exactly zero-variance initial time is removed.

times = np.asarray(nominal_solution.ts)
num_timesteps = len(times)
sensitivity_times = times[1:]
analysis_seed = 2718


def analyze_indices(time_index):
    Si_x = analyze_sobol.analyze(
        problem,
        Y_x[:, time_index],
        calc_second_order=False,
        num_resamples=100,
        conf_level=0.95,
        print_to_console=False,
        seed=analysis_seed,
    )
    Si_v = analyze_sobol.analyze(
        problem,
        Y_v[:, time_index],
        calc_second_order=False,
        num_resamples=100,
        conf_level=0.95,
        print_to_console=False,
        seed=analysis_seed,
    )
    return Si_x["S1"], Si_x["ST"], Si_v["S1"], Si_v["ST"]


results = [analyze_indices(k) for k in range(1, num_timesteps)]
S1_x, ST_x, S1_v, ST_v = (np.asarray(values) for values in zip(*results))

First-order effects#

The estimated first-order curves show the fraction of output variance assigned to each parameter by itself.

line_styles = ["-", "--", "-.", ":", (0, (3, 1, 1, 1))]
fig, axes = new_figure(size="full_tall", nrows=2, ncols=1, sharex=True)

for i, (label, line_style) in enumerate(zip(parameter_labels, line_styles)):
    axes[0].plot(sensitivity_times, S1_x[:, i], label=label, linestyle=line_style)
    axes[1].plot(sensitivity_times, S1_v[:, i], label=label, linestyle=line_style)

axes[0].set(ylabel=r"Estimated first-order index $\widehat{S}_i$", title=r"Displacement $x(t)$")
axes[1].set(xlabel=r"Time $t$", ylabel=r"Estimated first-order index $\widehat{S}_i$", title=r"Velocity $v(t)$")
axes[0].legend(ncol=5, frameon=False, loc="center right")
label_panels(axes, fontweight="bold")
finalize_axes(axes, keep_box=False);
Estimated first-order Sobol indices for each Duffing parameter versus time, shown separately for displacement and velocity.

The estimates indicate that, immediately after \(t=0\), the forcing amplitude \(\gamma\) has the largest first-order contribution. At later times, the forcing frequency \(\omega\) dominates much of the response. This ranking is conditional on the narrow uniform input ranges chosen above.

For the exact indices, the shortfall of \(\sum_i S_i\) from one is the fraction of variance assigned to interactions. The next plot shows the corresponding finite-design estimates.

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.axhline(1.0, color="0.6", linewidth=1.0, label="Exact no-interaction value")
ax.plot(sensitivity_times, np.sum(S1_x, axis=1), label=r"$x(t)$")
ax.plot(sensitivity_times, np.sum(S1_v, axis=1), label=r"$v(t)$", linestyle="--")
ax.set(xlabel=r"Time $t$", ylabel=r"$\sum_i \widehat{S}_i$")
ax.legend(frameon=False)
ax.title.set_text("Sum of estimated first-order indices")
finalize_axes(keep_box=False);
The sum of estimated first-order Sobol indices versus time for Duffing displacement and velocity, compared with the no-interaction reference value one.

Finite-design estimates need not satisfy the population bounds exactly. A small deficit or overshoot can be estimation error, and convergence as \(N\) increases need not be monotone. A persistent shortfall across larger designs and independent scrambles is evidence of interaction variance. Total-effect indices provide a second view by assigning each interaction to every participating input.

fig, axes = new_figure(size="full_tall", nrows=2, ncols=1, sharex=True)

for i, (label, line_style) in enumerate(zip(parameter_labels, line_styles)):
    axes[0].plot(sensitivity_times, ST_x[:, i], label=label, linestyle=line_style)
    axes[1].plot(sensitivity_times, ST_v[:, i], label=label, linestyle=line_style)

axes[0].set(ylabel=r"Estimated total-effect index $\widehat{S}_{T_i}$", title=r"Displacement $x(t)$")
axes[1].set(xlabel=r"Time $t$", ylabel=r"Estimated total-effect index $\widehat{S}_{T_i}$", title=r"Velocity $v(t)$")
axes[0].legend(ncol=5, frameon=False, loc="upper center")
label_panels(axes, fontweight="bold")
finalize_axes(axes, keep_box=False);
Estimated total-effect Sobol indices for each Duffing parameter versus time, shown separately for displacement and velocity.

For each input, \(S_{T_i}\) includes its main effect and every interaction involving that input. Comparing the first-order and total-effect sums gives a weighted summary of interaction variance because an interaction involving \(r\) inputs contributes to \(r\) total-effect indices.

fig, axes = new_figure(size="full_tall", nrows=2, ncols=1, sharex=True)

for ax, S1, ST, title in zip(
    axes,
    (S1_x, S1_v),
    (ST_x, ST_v),
    (r"Displacement $x(t)$", r"Velocity $v(t)$"),
):
    ax.axhline(1.0, color="0.6", linewidth=1.0)
    ax.plot(sensitivity_times, np.sum(S1, axis=1), label=r"$\sum_i \widehat{S}_i$")
    ax.plot(sensitivity_times, np.sum(ST, axis=1), label=r"$\sum_i \widehat{S}_{T_i}$", linestyle="--")
    ax.set(ylabel="Sum of indices", title=title)
    ax.legend(frameon=False, ncol=2)

axes[1].set_xlabel(r"Time $t$")
label_panels(axes, fontweight="bold")
finalize_axes(axes, keep_box=False);
Sums of first-order and total-effect Sobol indices versus time for Duffing displacement and velocity, revealing the contribution of interactions.

For the exact indices, \(\sum_i S_{T_i}\ge 1\), with equality if and only if all interaction variances vanish. An interaction involving \(r\) inputs is counted \(r\) times in this sum, while \(S_{T_i}-S_i\) collects the interactions involving input \(i\). Across most of the time interval, the first-order and total-effect sums are close and remain near one, so main effects dominate. The wider gaps at several isolated times indicate modest but nonzero interaction contributions for the chosen narrow input ranges. Small negative estimates of \(S_{T_i}-S_i\) or other bound violations are finite-design error, not negative interaction variance.

Exercise#

Repeat the calculation with \(N=1024\) and with an independent scramble seed. Compare the first-order and total-effect curves away from \(t=0\) and the very small-variance early times. Which qualitative features persist across both changes?

Local and global sensitivity#

Local sensitivities describe infinitesimal parameter changes near \(\boldsymbol{\theta}_0\). Sobol indices instead allocate output variance over a specified input distribution. The two views are complementary: their rankings can differ because they answer different questions. The next section develops polynomial-chaos expansions, which represent a model response with orthogonal polynomials and can make repeated uncertainty calculations more efficient.