Abstract

Monte Carlo sampling is used to characterize parameter uncertainty in linear single- and two-degree-of-freedom free-vibration benchmarks with fixed initial conditions. Each model uses 100,000 trials: uniform mass, stiffness, and damping inputs for SDOF and independent, unfiltered normal inputs for 2DOF. The SDOF undamped natural-frequency mean and sample standard deviation are 1.38036 and 0.02817 rad/s. An analytical audit reveals that the original forward Euler calculation gives 67.384% peak-displacement exceedance at 13.5 m, compared with 39.842% for the analytical response with identical samples. The original figure is retained as a documented numerical approximation; the analytical estimate supersedes its threshold result. The lower and upper 2DOF natural frequencies have means 1.80020 and 5.72171 rad/s and sample standard deviations 0.06309 and 0.18961 rad/s. Nominal damped eigencomponents have amplitude half-lives of 7.88 and 0.385 s. The 2DOF population retains 4,805 negative first-damper draws, and its central response bands do not describe an entirely passive population. Undamped modal illustrations are distinguished from damped transients. The study demonstrates uncertainty propagation and numerical verification, rather than validated manufacturing yield or service reliability. Differences between the two models do not isolate distribution-family effects.

Keywords

Monte Carlo Simulation, Uncertainty Quantification, Vibration Analysis, MDOF Systems, Natural Frequency, Transient Response

Share and Cite:

1. Introduction

Monte Carlo simulation is a numerical method that uses repeated random sampling to investigate systems whose inputs or governing conditions are uncertain. Rather than evaluating a model for one fixed parameter set, the analyst repeatedly solves it using samples drawn from specified probability distributions or admissible ranges. The resulting ensemble can be used to estimate moments, quantiles, exceedance probabilities, and other characteristics of the output. Classical illustrations include estimating an area by random point sampling and Buffon’s needle, in which repeated random needle drops provide an estimate of π [1]. These examples capture the central principle: repeated evaluations can approximate quantities that may be difficult to obtain analytically.

Monte Carlo methods are useful in engineering because real components and operating environments are not perfectly deterministic. Mass, stiffness, damping, geometry, material properties, and loads may vary because of manufacturing tolerances, assembly differences, environmental conditions, aging, and measurement error [2]-[4]. A deterministic model evaluated at nominal values provides a reference response, but it does not describe the distribution of responses across a population of nominally similar systems. By assigning defensible probability models to uncertain inputs and repeatedly solving the governing equations, Monte Carlo simulation propagates input variability into quantities directly relevant to design and operation [5]-[7].

This paper examines uncertainty propagation and numerical verification in linear SDOF and 2DOF free-vibration benchmarks. Reliability analysis asks whether a system meets a defined performance requirement over a specified duration [8], whereas quality control concerns conformance to manufacturing requirements [9]. The present results describe modeled frequency and transient-response distributions and illustrative threshold exceedance. They do not establish manufacturing yield or service reliability for a physical product.

Uncertainty in structural properties is a major source of variation in natural frequencies and vibration response. A stochastic finite-element formulation based on Neumann expansion treats random structural parameters together with random excitation, using direct Monte Carlo calculations to assess dynamic estimates [2]. For composite plates, material-property variability affects the statistics of nonlinear free vibration, and the sensitivity of frequency variability changes with plate thickness and oscillation amplitude [3]. For damaged laminated beams, spatially variable micromechanical properties and cracking introduce further modal-frequency variability; surrogate predictions have been checked against finite-element and Monte Carlo calculations, including sensitivity and noise effects [4]. These examples motivate describing a manufactured population by response distributions rather than by a single nominal frequency.

Several complementary approaches address the information required to characterize uncertain systems. Stochastic model updating uses response observations within an inverse Monte Carlo procedure to infer uncertain model parameters, as illustrated for a three-degree-of-freedom system [10]. Interval-process sampling instead provides response bounds when a probability distribution cannot be justified [11]. These methods answer different questions: probabilistic updating estimates uncertainty from observations and a model, whereas interval analysis describes admissible ranges. In particular, a response bound is not itself a failure probability; probability estimates require a defensible probabilistic description of the inputs.

Nonlinear vibration introduces mechanisms that are not represented by constant-stiffness linear models. Pinching hysteresis can incorporate stiffness, strength, and pinching degradation, with Monte Carlo simulation serving as a reference for approximate response statistics [12]. Contact loss in switching apparatus creates a piecewise-linear two-degree-of-freedom problem in which vibration can interrupt intended mechanical contact [13]. Parametric random vibration of an axially moving laminated shape-memory-alloy beam extends the problem to moving structures with material-dependent dynamics [14]. These studies illustrate how uncertainty analysis can be adapted to physical mechanisms, such as degradation, intermittent contact, and axial motion.

The evolution of response uncertainty is also important. Radial-basis-function neural networks have been used to approximate transient response probability distributions in nonlinear stochastic systems [15]. Here, random mechanical parameters generate an ensemble of deterministic free responses from fixed initial states. The examples involve parameter uncertainty rather than random external excitation. Mass and stiffness determine the undamped natural frequencies; damping and the initial state determine the transient response.

Repeated dynamic calculations can become expensive, particularly when the quantity of interest is a rare failure event [16]. Sampling procedures such as Russian roulette and splitting direct computational effort toward low-probability response regions [16]. A different strategy approximates the mapping from uncertain inputs to outputs: polynomial-chaos models of natural frequencies, damping ratios, and other modal quantities can support predictions of uncertain frequency-response-function bounds [5], while multi-fidelity co-kriging combines Euler-Bernoulli and Timoshenko beam calculations to estimate uncertain natural frequencies with fewer expensive evaluations [17]. These approaches reduce computational burden through targeted sampling, modal approximation, or the combination of models with different fidelity.

Accuracy analysis is especially important when the output is a small failure probability, because modest errors near a limit-state boundary may materially affect the estimate. Optimized equivalent linearization combines a physical linear surrogate with Monte Carlo correction methods, including control variates and importance sampling, to estimate a specified nonlinear response quantity [18]. This issue is directly relevant to reliability analysis: computational efficiency is useful only when the approximation remains accurate in the response region governing failure. Direct Monte Carlo simulation is therefore a practical benchmark for the relatively small one- and two-degree-of-freedom models considered here [16]-[18].

Reliability over a finite duration requires determining whether a response crosses a prescribed limit at any time during that interval. Random excitation and random system parameters can be treated jointly by combining a Girsanov transformation with subset simulation [19]. Simulation-based reliability procedures also extend to laboratory testing, including automotive vibration under road-load conditions on a four-post rig [20]. Subset simulation with Markov-chain splitting provides another experimental route to rare-event estimation, illustrated for a building frame under biaxial random support motion [21]. These examples connect computational estimates with dynamic quantities that can, in principle, be measured and checked experimentally.

Response variability, exceedance probability at a particular instant, and failure probability over an entire event are distinct quantities. A narrow response distribution at one time does not establish satisfactory performance throughout an operating interval. A time-dependent reliability calculation therefore needs an explicit limit state, a defined observation duration, and a trial-by-trial check of the maximum or other relevant response functional [19]-[21]. Experiments and application-specific evidence are also needed to establish whether the assumed input distributions, model form, and failure criterion represent the system of interest.

Vibration-control studies have compared viscous dampers, tuned mass dampers, and tuned inerter dampers using reliability metrics under Gaussian and non-Gaussian loading [22]. Pounding tuned mass dampers introduce nonlinear impact behavior [23], and surrogate reliability methods extend to systems with multiple degrees of freedom [24]. In rotating machinery, uncertain misalignment, damping, and support stiffness affect response [25], while transient balancing has been investigated using Monte Carlo and interval methods [26]. Probabilistic critical-speed analysis connects uncertain dynamic parameters with operating limits [6]. Related studies address probabilistic design constraints [7], combined wear and vibration failure [27], and fatigue-damage estimation under stochastic aerospace loading [28]. These applications require excitation, degradation, or failure models beyond the unforced benchmarks considered here.

Passive damping and degradation are distinct: damping dissipates energy during a vibration event, whereas degradation changes properties such as stiffness or strength over service life. This study keeps each sampled parameter constant during its trial and includes no degradation model. Its contribution is a reproducible numerical benchmark that distinguishes frequency variability, transient response, and threshold exceedance, and demonstrates how integration error can materially change an exceedance estimate even when modal-frequency statistics are accurate.

2. Application in Reliability Engineering and Quality Control

Monte Carlo simulation provides a practical method for examining how uncertainty influences the performance of engineering systems. Mechanical models often begin with fixed values for mass, stiffness, damping, material strength, and applied loading. These values describe a nominal system, but actual components and operating conditions vary. Manufacturing processes introduce differences between products, while changing loads, environmental conditions, and degradation introduce variation during service. Monte Carlo simulation incorporates these uncertainties by repeatedly sampling input values and calculating the resulting system response. The collection of results describes a range of possible outcomes and allows engineers to estimate the probability of exceeding a performance limit.

2.1. Application in Reliability Engineering

Reliability engineering evaluates whether a component or system performs its intended function under specified conditions for a specified duration. A reliability calculation therefore requires a limit-state function, a sign convention, a probabilistic model for uncertain inputs, and a clearly defined observation event. For example, the illustrative strength-load equation could be written as:

$g = S - L$ (1)

where S is strength and L is the applied load in the same units. The convention used here is $g>0$ for satisfactory performance and $g\le 0$ for failure. For an instantaneous event the limit state is evaluated at the specified time; for a duration-based vibration event the limit state can instead be written $g=x_{\lim }-{\max }_{t\in \left[0,T\right]}\left|x\left(t\right)\right|$ , so failure occurs if the peak response over the full interval reaches or exceeds the allowable limit $x_{\text{lim}}$ .

In Monte Carlo analysis, strength and applied load are assigned probability distributions based on measurements, relevant evidence, or stated assumptions. Each trial samples these variables and evaluates the limit state. Their dependence must be retained when shared physical factors influence both quantities.

After $N$ independent trials, the estimated probability of failure and reliability are

$\widehat{p}=\frac{N_f}{N},\widehat{R}=1-\widehat{p}$(2)

Here $N_f$ is the number of trials classified as failures. For example, if 500 failures occur in 100,000 simulated trials,

$\widehat{p}=\frac{500}{100000}=0.005,\widehat{R}=0.995$(3)

For this illustrative count, the estimated failure probability is 0.5% and the estimated reliability is 99.5%. An approximate normal 95% binomial confidence interval is 0.456% - 0.544%. This example is not a result of the vibration simulations. Increasing the trial count reduces sampling uncertainty but cannot correct an unsuitable input model, an unjustified failure criterion, or numerical integration bias.

Rare failure events can require many trials before enough failures occur for a stable estimate. Advanced Monte Carlo procedures and importance sampling address this difficulty [16] [18]; subset simulation has also been used for time-dependent reliability [19] [21]. Such methods are especially useful when each dynamic calculation is expensive.

2.2. Application to Mechanical Vibration

Vibration analysis extends the reliability problem because system response changes throughout an operating event. A suspension assembly, for example, may contain uncertain masses, spring stiffnesses, and damping coefficients while experiencing variable road excitation. Each Monte Carlo trial represents one possible assembly and one possible loading history. Solving its equations of motion produces displacement, velocity, acceleration, or other response histories that can be compared with specified limits.

A vibration failure criterion may concern displacement clearance, acceleration, transmitted force, or stress, depending on the engineering function. Peak-displacement exceedance is distinct from fatigue damage, which requires stress-cycle and material-fatigue information. Monte Carlo vibration-damage estimation provides an example of the latter [28]. The thresholds in this paper are numerical illustrations rather than application-specific failure limits.

Uncertainty in mass and stiffness changes natural frequencies and mode shapes. Under external excitation, these changes affect proximity to resonance, while damping influences amplification and free-response decay. Studies of uncertain frequency-response functions, composite structures, and rotating machinery illustrate these effects [3] [5] [6]. The present models are unforced; resonance under operating loads is not simulated.

Monte Carlo results can provide more than a single pass/fail estimate. Response distributions characterize the range of predicted vibration levels, and sensitivity analyses can identify which uncertain inputs most strongly affect a selected performance measure. In manufacturing, this information can help evaluate tolerance allocation and inspection strategies; in service, it can help identify operating conditions or components warranting further investigation. Such uses remain conditional on model adequacy and on validation of the input distributions and response predictions. When repeated dynamic calculations are costly, surrogate models may improve efficiency, but their accuracy should be assessed over the region relevant to the engineering decision [17] [24].

2.3. Application in Quality Control

Quality control evaluates whether manufactured products satisfy defined requirements. Even for a common design, machining, material processing, assembly, and measurement introduce variation in dimensions and mechanical properties. Monte Carlo simulation propagates this variation through a physical model to estimate the fraction of products meeting dimensional or functional specifications. For vibration-sensitive products, functional specifications may include a natural-frequency band, a maximum transient displacement, or a settling-time requirement, rather than dimensional tolerances alone.

Consider a manufactured dimension D that follows a normal distribution with a mean of 10.00 mm and a standard deviation of 0.05 mm. If acceptable values range from 9.90 to 10.10 mm, the limits are two standard deviations below and above the mean. Under this assumed distribution, approximately 95.45% of the population lies within specification. A sufficiently large Monte Carlo simulation should approach this analytical yield. This simple example also provides a useful check that the sampling and acceptance calculations have been implemented correctly.

As an illustrative manufacturing model, a spring stiffness could be sampled uniformly between 950 and 1050 N/m about a nominal value of 1000 N/m. These bounds are a hypothetical quality-control example and are separate from the SDOF and 2DOF input assumptions used below.

Each trial samples stiffness and evaluates the resulting response against a defined functional requirement. If mass and damping also vary, their joint distribution must be represented. Repeating the calculation estimates a modeled conformance fraction; it does not establish production yield without justified inputs and acceptance limits.

The assumed distribution is a substantive modeling choice. A uniform distribution assigns equal density throughout a specified interval, but specification limits alone do not establish that values are uniformly distributed. Production measurements may instead support a concentration near the target, skewness, batch effects, or dependence among component properties. Shared material lots or manufacturing operations can induce correlations; treating such variables as independent when they are not may misrepresent the output variability. Input models should therefore be justified by data where possible and identified as assumptions otherwise.

Quality-oriented simulation can be used to examine proposed tolerances and acceptance rules before implementation. A narrower stiffness tolerance may reduce natural-frequency variability, but whether that reduction is beneficial depends on the functional requirement and on the cost or feasibility of achieving the tolerance. Probabilistic design studies illustrate how uncertain responses can be incorporated into constraints [7]. Monte Carlo analysis makes it possible to compare alternatives in terms of predicted yield and response distributions, provided that the alternatives are evaluated under consistent assumptions.

2.4. Connecting Manufacturing Quality with Service Reliability

Manufacturing conformance and service reliability describe different stages of a product’s life. A component can pass factory inspection and later experience wear, stiffness loss, or other changes that alter its dynamic response. Conversely, some initial variability may have little consequence under the actual service conditions. Wear-dependent reliability analysis emphasizes the need to model how properties and failure risk evolve during operation [27]. Thus, an acceptance decision should not automatically be interpreted as evidence of reliability over the full service life.

For products accepted at the factory, the relevant reliability measure is $R\left(T|A\right)=P\left(\text{no failure through service age}T|\text{factory}\text{acceptance}A\right)$ . Here A denotes factory acceptance and $T$ denotes the service age being evaluated. A Monte Carlo study can estimate this quantity by first simulating the manufactured population, applying an acceptance rule, and then evaluating the service performance of the accepted units. Each unit retains its sampled manufacturing properties while its condition evolves according to the assumed degradation model.

This distinction motivates the proposed research direction. Factory screening based on natural frequencies can be compared with screening that also uses vibration-decay information. To isolate the effect of the screening metric, the comparison should maintain the same manufacturing acceptance yield; service reliability can then be estimated for the accepted populations under a specified degradation and operating model. The current numerical examples motivate this comparison but do not themselves establish that combined screening reduces later failure. That claim requires a completed conditional-reliability analysis and, ultimately, experimental validation.

Figure 1. Original Monte Carlo workflows for (a) SDOF and (b) MDOF vibration models. Panel (a) describes the explicit Euler calculation underlying the retained SDOF response plots; an analytical audit of that calculation is reported in Section 3. Panel (b) describes matrix assembly, modal-frequency calculation, and an RK45 implementation. In its second-order equation, the displacement vector should be denoted x: $M\ddot{x}+C\dot{x}+Kx=0$ . The interleaved displacement-velocity state z has 2n components and instead satisfies ż = Az. The figure is retained without graphical modification; this caption resolves the notation used in panel (b). N denotes the number of trials and n the number of mechanical degrees of freedom.

As shown in Figure 1, the workflows separate parameter sampling, dynamic solution, and ensemble post-processing. The RK45 implementation and the exact state-transition calculation used in the local 2DOF summaries are distinguished in Section 4. Neither the workflow nor a large trial count alone establishes numerical convergence.

3. Single-Degree-of-Freedom Vibration

A single-degree-of-freedom (SDOF) system is the simplest lumped-parameter vibration model. A degree of freedom is an independent coordinate required to describe the system configuration; one coordinate is therefore sufficient to describe the motion of an SDOF system.

A common SDOF representation consists of a mass connected to a spring and a damper as shown in Figure 2(a). The mass represents the moving component, the spring supplies a restoring force, and the damper dissipates mechanical energy. A simplified suspension or an isolated machine component may be approximated this way when motion in one dominant direction is of interest. The model is an idealization: its usefulness depends on whether neglected modes and nonlinear effects are unimportant for the question being studied.

Figure 2. Vibration-model schematics. (a) SDOF system with mass m (kg), stiffness k (N/m), damping coefficient c (N·s/m), and displacement x(t) (m). (b) Coupled 2DOF system: spring k₁ and damper c₁ connect mass m₁ to the fixed support; spring k₂ and damper c₂ connect m₁ and m₂. Mass m₂ has no separate support connection. (c) Generalized chain of n masses coupled by adjacent spring-damper elements, with the first mass connected to the fixed support. The schematics define connectivity; their drawn geometry does not set the benchmark displacement scale.

The three panels in Figure 2 establish the physical models to which the workflow is applied: Figure 2(a) is the single-coordinate mass-spring-damper benchmark, Figure 2(b) is the coupled two-mass system used for the detailed 2DOF simulations, and Figure 2(c) depicts the generalized MDOF arrangement. In the 2DOF schematic, the connecting spring and damper transmit forces between the masses, so their motions are coupled rather than independent.

After an initial displacement and velocity are imposed, the unforced system oscillates and, for positive damping, dissipates energy. Its response follows from the governing equation of motion.

The spring and damping forces act in the direction opposite to the displacement and velocity, respectively. According to Newton’s second law,

${\displaystyle \sum F}=ma$ (4)

The restoring force produced by the spring is described by Hooke’s law as

$F_k=-kx$ (5)

where k is the spring stiffness. The negative sign indicates that the spring force acts opposite to the direction of displacement. Similarly, the viscous damping force is

$F_c=-c\dot{x}$ (6)

Here c is the viscous damping coefficient and $\dot{x}$ is velocity. For positive c, the damping force opposes motion. With no external force, Newton’s second law gives

$m\ddot{x}=-c\dot{x}-kx$ (7)

Equivalently,

$m\ddot{x}+c\dot{x}+kx=0$ (8)

Equation (8) is a linear second-order homogeneous differential equation with constant coefficients, describing a linear time-invariant system. Substituting the trial solution x(t) = exp(rt) gives

$m{r}^2{\text{e}}^{rt}+cr{\text{e}}^{rt}+k{\text{e}}^{rt}=0$ (9)

Since exp(rt) is nonzero, the characteristic equation is

$m{r}^2+cr+k=0$ (10)

Its roots are

$r=\frac{-c\pm \sqrt{c^2-4mk}}{2m}$ (11)

All SDOF samples are underdamped and have positive mass, stiffness, and damping. Thus c² < 4mk, and the characteristic roots form a complex-conjugate pair:

$r=\frac{-c}{2m}\pm i{\omega }_d$ (12)

The damped natural angular frequency is

${\omega }_d=\sqrt{\frac{k}{m}-{\left(\frac{c}{2m}\right)}^2}$ (13)

Hence, the displacement response can be written as

$x\left(t\right)={\text{e}}^{-\frac{c}{2m}t}\left[C_1\cos \left({\omega }_{d}t\right)+C_2\sin \left({\omega }_{d}t\right)\right]$ (14)

The constants are C₁ = x₀ and C₂ = (v₀ + αx₀)/ω_d, where α = c/(2m). The exponential term exp(−αt) describes amplitude decay. This expression is valid for the underdamped case; the SDOF ranges below satisfy that condition.

A deterministic vibration calculation assumes that the parameters retain their nominal values. Actual systems may depart from those values: spring stiffness can vary with manufacturing tolerances, effective mass can change with payload, and damping can depend on temperature, wear, assembly, or operating amplitude. Representing these quantities as uncertain inputs allows the model to describe a population of possible systems rather than only a nominal specimen.

The SDOF benchmark uses nominal mass 1050 kg, stiffness 2000 N/m, and damping coefficient 50 N·s/m, with fixed initial displacement x(0) = 5 m and velocity $\dot{x}\left(0\right)$ = −18 m/s. These deliberately large values define an illustrative mathematical benchmark, not a calibrated suspension or machine-mount application. The initial kinetic energy allows the subsequent absolute displacement to exceed 5 m while passive damping reduces total mechanical energy.

For N = 100,000 trials, independent uniform draws use m ∈ [997.5, 1102.5] kg, c ∈ [47.5, 52.5] N·s/m, and k ∈ [1900, 2100] N/m. NumPy default_rng with seed 0 supplies draws in the order m, c, k for each trial. The original calculation advances displacement and velocity by simultaneous forward Euler updates with step 0.01 s over 0 - 120 s; displacement summaries are stored every 0.1 s. Figure 3 retains this original Euler output. The analytical audit uses the same sampled parameters and fixed initial state.

Figure 3. Original SDOF outputs retained for comparison with the analytical audit. (a) Euler displacement median and pointwise 5% - 95% population band over 0 - 120 s. (b) Empirical CDF of peak absolute Euler displacement with illustrative threshold 13.5 m; the displayed 67.4% is the original Euler exceedance estimate, not the corrected analytical result. (c) Undamped natural-angular-frequency histogram and fitted normal density. The analytical exceedance fraction is 39.842% for the same samples, as reported in the text. Frequency statistics in (c) depend only on sampled mass and stiffness and are unaffected by the transient integration error.

Figure 3 separates response variability, peak-displacement exceedance, and frequency variability. Its pointwise bands summarize the central 90% of Euler responses at each stored time; they are neither confidence intervals for the median nor simultaneous coverage bands for entire trajectories. The retained transient plots document the original numerical approximation and are not used as a converged reference solution.

For the nominal SDOF parameters, ω_d = 1.3799 rad/s and the analytical amplitude envelope is 13.8891 exp(−0.0238095t) m, giving 6.7993 m at 30 s. Agreement with this envelope must be assessed quantitatively; an envelope calculation alone does not validate a numerical trajectory. The trial-by-trial audit below shows that the original Euler peak estimates are biased upward.

The original Euler calculation gives 67,384 exceedances among 100,000 trials at 13.5 m, or 67.384%. The analytical solution gives 39,842 exceedances, or 39.842%, on the same 0 - 120 s interval and for the same parameter samples. For each analytical trial, the extrema are obtained by setting $\dot{x}\left(t\right)$ = 0 in Equation (14) and comparing all relevant extrema with the initial displacement. Positive damping ensures decreasing amplitudes of successive extrema. Relative to the analytical peaks, the Euler peaks have a mean relative upward bias of 1.337%, which changes the exceedance fraction by 27.542 percentage points. The analytical result therefore supersedes the Euler value as the model-based SDOF threshold estimate. The unchanged Figure 3(b) remains explicitly labeled in its caption as the original Euler result.

The SDOF undamped natural frequencies have mean 1.38036 rad/s and sample standard deviation 0.02817 rad/s, with a sampled range of 1.31298 - 1.45061 rad/s. Standard deviations reported here and below use the sample convention with denominator N − 1; the normal-density fits in the original figures use their fitted scale parameters. The nominal frequency sqrt (2000/1050) is 1.38013 rad/s. Because frequency depends nonlinearly on mass and stiffness, the frequency at the mean inputs need not equal the mean sampled frequency.

Wilson 95% binomial confidence intervals are 67.093% - 67.674% for the original Euler exceedance estimate and 39.539% - 40.146% for the analytical estimate. These intervals quantify Monte Carlo sampling uncertainty conditional on the assumed model; they exclude uncertainty in the input distributions, model form, and threshold. Their separation shows that increasing N cannot resolve the integration bias. Neither result is a validated physical failure probability.

4. Multi-Degree-of-Freedom Vibration

Many mechanical and civil-engineering systems cannot be represented adequately by a single moving mass. Vehicles contain components that move relative to one another, while buildings and bridges can exhibit motion at multiple locations and in several modes. Such systems are represented by multi-degree-of-freedom (MDOF) models, in which the coordinates are dynamically coupled (Figure 2(c)).

A two-degree-of-freedom system has two independent coordinates, x₁(t) and x₂(t), while their motions are dynamically coupled through the connecting spring and damper (Figure 2(b)).

The motion of one mass influences the other through the connecting elements. Consequently, the model contains coupled differential equations rather than one independent equation for each isolated mass. This coupling produces system-level modes whose frequencies and shapes depend on the complete mass, damping, and stiffness arrangement.

For the two-mass connectivity in Figure 2(b), the first equation is

$m_1{\ddot{x}}_1+\left(c_1+c_2\right){\dot{x}}_1-c_2{\dot{x}}_2+\left(k_1+k_2\right)x_1-k_2{x}_2=0$ (15)

For the second mass,

$m_2{\ddot{x}}_2+c_2{\dot{x}}_2-c_2{\dot{x}}_1+k_2{x}_2-k_2{x}_1=0$ (16)

These equations differ from the SDOF equation because the displacement and velocity of one mass contribute to the forces acting on the other. For example, motion of the first mass deforms the connecting spring and damper, transmitting forces to the second mass. The response of either coordinate therefore depends on the coupled system, not solely on its own mass or support.

For numerical simulation, these second-order differential equations can be rewritten as a system of first-order differential equations. New state variables can be defined as

$\{\begin{array}{c}{z}_1=x_1\\ z_2={\dot{x}}_1\\ z_3=x_2\\ z_4={\dot{x}}_2\end{array}$ (17)

Therefore,

${ż}_1=z_2$ ${ż}_3=z_4$ (18)

Using the newly defined variables, the acceleration of the first mass and the second mass becomes

${ż}_2=\frac{-\left(c_1+c_2\right)z_2+c_2{z}_4-\left(k_1+k_2\right)z_1+k_2{z}_3}{m_1}$ (19)

${ż}_4=\frac{-c_2{z}_4+c_2{z}_2-k_2{z}_3+k_2{z}_1}{m_2}$ (20)

The coupled equations are written in state-space form to describe the evolution of the system through its current displacements and velocities. These variables constitute the state; together with the mechanical parameters and applied forces, they determine the subsequent motion under the model. This representation is convenient for numerical integration and for repeatedly solving the system with different Monte Carlo samples.

The support spring k₁ and damper c₁ act on mass m₁, while k₂ and c₂ transmit forces between the two masses. This connectivity gives the sums k₁ + k₂ and c₁ + c₂ in the first acceleration equation and only coupling terms in the second. The state order in Equations (17)-(21) is [x₁, v₁, x₂, v₂]ᵀ.

The first-order equations can therefore be written as

$\frac{\text{d}}{\text{d}t}\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ z_4\end{array}\right]=\left[\begin{array}{cccc}0& 1& 0& 0\\ -\frac{k_1+k_2}{m_1}& -\frac{c_1+c_2}{m_1}& \frac{k_2}{m_1}& \frac{c_2}{m_1}\\ 0& 0& 0& 1\\ \frac{k_2}{m_2}& \frac{c_2}{m_2}& -\frac{k_2}{m_2}& -\frac{c_2}{m_2}\end{array}\right]\left[\begin{array}{c}{z}_1\\ z_2\\ z_3\\ z_4\end{array}\right]$ (21)

This formulation preserves the two mechanical degrees of freedom and uses four state variables. For fixed sampled coefficients, the linear response can be calculated either numerically with an ODE solver or by $z\left(t\right)=\exp \left[A\left(t-t_0\right)\right]z_0$ . The latter is an exact solution formula; evaluating the matrix exponential still involves floating-point computation.

Figure 4 documents the input distribution families. The SDOF output uses the uniform distributions in panel (b), and the 2DOF output uses the independent normal distributions in panel (c) and Table 1. Panels (a) and (d) illustrate alternative distribution shapes with matched means and standard deviations; they are not separate transient-response simulations. N (μ, σ²) denotes a normal distribution with mean μ and variance σ².

The 2DOF benchmark samples six uncertain parameters: m₁, m₂, k₁, k₂, c₁, and c₂. The initial displacement and velocity of each mass remain fixed.

Each uncertain parameter is assigned a probability distribution or physically justified range. A Monte Carlo trial draws one parameter combination and solves the coupled equations; subsequent trials use new combinations. The ensemble provides response histories for both masses and permits investigation of whether particular coordinates experience large or persistent motion. Potential applications include simplified vehicle suspensions, machinery supported by compliant mounts, and reduced-order representations of building or bridge components. These are motivating examples; the present two-mass model is not, by itself, a validated representation of any particular vehicle or civil structure.

The SDOF and MDOF models serve different purposes. The SDOF model isolates the basic relationships among mass, damping, stiffness, frequency, and displacement. The MDOF model captures interactions among coordinates and the resulting multiple modes. Monte Carlo sampling can be applied to either formulation to characterize response variability, but conclusions remain specific to the chosen model, uncertain inputs, and loading or initial conditions.

The physical parameter vector is θ = [m₁, m₂, k₁, k₂, c₁, c₂]ᵀ. For every trial, one independent draw of these six variables defines a time-invariant system. In the seeded code, standard-normal columns are mapped in the implementation order [k₁, k₂, c₁, c₂, m₁, m₂]; this order must be retained to reproduce the exact seeded ensemble.

Table 1 lists prescribed numerical-example means and standard deviations, and the six inputs are assumed independent. These values, the SDOF ranges, and the independence assumptions are used for illustration; they are not fitted to measurements or empirically calibrated manufacturing data. The accompanying

Figure 4. Input probability-density illustrations. (a) Illustrative normal distributions for the SDOF parameters. (b) Uniform inputs actually used for the SDOF model. (c) Unfiltered normal inputs actually used for the 2DOF model. (d) Illustrative uniform alternatives for the 2DOF parameters. Alternative uniform bounds are $\mu \pm \sqrt{3}\sigma$ , preserving the corresponding mean and standard deviation. The c₁ normal model in (c), and its illustrative uniform alternative in (d), allow negative damping values. Panels (a) and (d) are conceptual alternatives only. Density units are the reciprocal of the corresponding parameter units; these panels do not isolate distribution-family effects on vibration response.

Table 1. Input probability distributions for the two-degree-of-freedom model.

code contains no measurement input, distribution fitting, or correlation estimation, and does not document a separate physical basis for selecting the values.

The six parameters are sampled independently. This is a numerical modeling assumption rather than a measured absence of correlation. For each trial, the mass, damping, and stiffness matrices are

$M=\left[\begin{array}{cc}{m}_1& 0\\ 0& m_2\end{array}\right]$ (22)

$C=\left[\begin{array}{cc}{c}_1+c_2& -c_2\\ -c_2& c_2\end{array}\right]$ (23)

$K=\left[\begin{array}{cc}{k}_1+k_2& -k_2\\ -k_2& k_2\end{array}\right]$ (24)

The undamped modal properties of each sampled system are then determined by solving the generalized eigenvalue problem

$K{\phi }_r={\lambda }_{r}M{\phi }_r$ (25)

where λ_r is the eigenvalue associated with mode r, and φ_r is the corresponding eigenvector. In vibration analysis, the eigenvector represents the relative displacement pattern of the system and is therefore interpreted as the mode shape, or eigenmode.

The natural angular frequency associated with each eigenvalue is

${\omega }_r=\sqrt{{\lambda }_r}$ (26)

The two positive eigenvalues give the lower and upper undamped natural angular frequencies of the coupled assembly. Their eigenvectors describe undamped mode shapes. Damping is absent from this eigenproblem and therefore does not affect these frequencies, although it affects damped oscillation frequencies and transient decay.

The initial state of the specific 2DOF model is written as

$z_0={\left[x_1\left(0\right),v_1\left(0\right),x_2\left(0\right),v_2\left(0\right)\right]}^{\text{T}}$ (27)

For every 2DOF Monte Carlo trial, the fixed initial state is z₀ = [15, 5, −5, 5]ᵀ in the order [x₁(0), v₁(0), x₂(0), v₂(0)]ᵀ. Thus, x₁(0) = 15 m, v₁(0) = 5 m/s, x₂(0) = −5 m, and v₂(0) = 5 m/s. Differences between trial responses arise from the sampled mechanical parameters; the initial state is not sampled.

The equations of motion for each sampled system can be written in matrix form as

$M\ddot{x}+C\dot{x}+Kx=0$ (28)

where

$x=\left[\begin{array}{c}{x}_1\\ x_2\end{array}\right]$ (29)

For numerical integration, the second-order system is converted into first-order state-space form. The state vector is defined as

$z=\left[\begin{array}{c}{x}_1\\ v_1\\ x_2\\ v_2\end{array}\right]$ (30)

and the state equation can be written symbolically as

$\dot{z}=Az$ (31)

The original 2DOF script uses solve_ivp with RK45, relative tolerance 10⁻³, and absolute tolerance 10⁻⁶ over 0 - 120 s, recording 1201 times at 0.1 s spacing and raising an error on solver failure. The local plotting summaries associated with the retained figures were rebuilt from the same seed-0 parameter draws using the linear state transition exp(0.1A) instead of RK45. Consequently, the plotted-summary method is distinguished from the original solver implementation rather than treating the two as identical. No full-ensemble tolerance-convergence result is claimed for RK45. The reproduction materials specify both calculation paths.

The post-processing stores undamped frequency samples, displacement summaries, and peak displacement on the recorded time grid. The modal illustrations in Figure 5(c) & Figure 5(d) are obtained separately from the undamped mass-stiffness eigenproblem with the same fixed initial displacements and velocities. They are not a decomposition of the nonproportionally damped transient response.

The normal draws are retained without truncation, rejection, or resampling. For c₁ ~ N(50, 30²) N·s/m, the probability of c₁ < 0 is approximately 4.78%; the seeded 100,000-trial ensemble contains 4805 negative c₁ values and no nonpositive c₂ values. The seeded masses and stiffnesses are positive, although their unbounded normal models do not enforce positivity in general. A negative c₁ makes the damping matrix nonpassive; dynamic stability must be assessed from the state-matrix eigenvalues, not inferred solely from the sign of one coefficient. This study retains the population as a mathematical benchmark rather than representing an entirely passive mechanical population.

All 100,000 trials contribute to the pointwise medians and percentiles. The central bands exclude the most extreme 5% on each side at a given time and therefore can conceal extreme growth in a small fraction of responses.

Figure 5. Original 2DOF response and modal illustrations. (a), (b) Damped response medians and pointwise 5% - 95% population bands from the unfiltered parameter ensemble. (c) Undamped modal displacement summaries at 1 s for the lower mode, upper mode, and trial-level superposition. The labels report population medians rather than nominal mode-shape amplitudes. (d) Undamped modal time histories: lines are the displayed population medians, darker bands are the central 90%, and lighter bands are finite-sample minimum-maximum ranges. Panels (c), (d) deliberately omit damping and must not be interpreted as the damped free-decay solution in (a), (b). (e) Recorded-grid peak-displacement CDFs for both masses and the either-mass event, with illustrative threshold 16 m and exceedance 2.020%. The displayed horizontal range is restricted and does not contain the largest stored peaks from the nonpassive population. None of these thresholds is a validated failure limit.

Figure 5(a) & Figure 5(b) presents median displacement histories and pointwise 5% - 95% population bands for the two masses. The bands describe the central response population at each time, not the full range of sampled trajectories. In particular, a narrowing central band does not show that every trial decays, because the unfiltered damping population includes nonpassive cases.

The lower undamped mode is primarily in phase and the upper mode primarily out of phase. Trial-level modal contributions add linearly to give total motion, but medians of those contributions need not add to the median of their sum. This explains small differences between the rounded modal summaries and the total summary in Figure 5(c). Because the nominal damping is nonproportional, the undamped illustration is separate from the damped response and its decay rates.

At the prescribed mean parameter values, the two undamped natural angular frequencies are approximately 1.7981 rad/s and 5.7058 rad/s. The lower mode involves both masses moving in the same direction, whereas the upper mode involves opposing motion. These frequencies describe modes of the complete assembly rather than separate frequencies assigned to individual masses.

Figure 6(a) displays the marginal distributions of the lower and upper natural angular frequencies, including the fitted normal curves, and Figure 6(b) shows their joint distribution across trials. The marginal panels describe the spread of each mode separately; the joint panel additionally reveals how the two modal frequencies co-vary within the same sampled assembly. These results concern undamped modal frequencies and should be distinguished from the damping-dependent transient decay discussed above.

The lower and upper frequencies have means 1.80020 and 5.72171 rad/s and sample standard deviations 0.06309 and 0.18961 rad/s, respectively. Their coefficients of variation are 3.50% and 3.31%, so proportional uncertainty is similar despite different absolute spreads. Frequency variability describes differences

Figure 6. Undamped natural-angular-frequency results for the 2DOF ensemble. (a) Marginal density histograms and fitted normal curves for the lower and upper modes. (b) Joint frequency histogram rendered with hexagonal bins; shading represents trial count on a logarithmic color scale, rather than an individual trajectory or a normalized joint density. Each trial contributes one frequency pair. The upward trend indicates positive association between the two modal frequencies under the specified independent parameter inputs.

between sampled assemblies; the parameters and undamped frequencies remain fixed within each trial.

At the mean mass and stiffness values, the undamped frequencies are 1.79813 and 5.70580 rad/s. Their difference from the ensemble means reflects nonlinear transformation and finite sampling. For the nominal damped state matrix, the eigenvalues are approximately −0.08801 ± 1.79684i and −1.79883 ± 5.41196i s⁻¹. Each conjugate pair has an amplitude decay half-life ln(2)/|Re(s)|, where s is the damped state-matrix eigenvalue, giving 7.88 s and 0.385 s. These rates describe the nominal damped eigencomponents; they do not describe every sampled system or the undamped curves in Figure 5(d).

For each 2DOF trial, the stored peak is the maximum of |x₁| or |x₂| across the 1201 recorded times. At 16 m, the either-mass event occurs in 2020 of 100,000 trials, giving 2.020% with a Wilson 95% binomial interval of 1.935% - 2.109%. A recorded-grid maximum can underestimate the continuous-time peak, and this interval excludes time-grid error and model uncertainty. Extreme responses outside the restricted horizontal range in Figure 5(e) remain included in the calculation. These results apply only to the unfiltered mathematical population; introducing nonnegative damping would change the input model and require a new transient analysis.

5. Influence of Input Probability Distributions

The probability distribution assigned to an uncertain parameter determines how values are sampled during a Monte Carlo simulation. A uniform distribution assigns equal probability density to all values within specified lower and upper bounds, whereas a normal distribution concentrates values near its mean and assigns progressively lower probability to values farther from the mean. The appropriate distribution depends on the physical parameter being represented, the available measurements, and any known physical limits. Figure 4 illustrates the shapes of the normal and uniform input distributions considered in this study and shows how the uncertain mechanical parameters are represented under each distribution family.

The SDOF simulation uses uniform distributions for mass, stiffness, and damping. The mass is sampled between 997.5 and 1102.5 kg, the stiffness between 1900 and 2100 N/m, and the damping coefficient between 47.5 and 52.5 N·s/m. Because these distributions are uniform, every value within the specified interval has the same probability density. Figure 4(b) shows the input distributions actually used in the SDOF Monte Carlo simulation. Figure 4(a) is included only to illustrate how the same types of SDOF parameters would appear under a normal-distribution representation and is not used to generate the reported SDOF response results.

The 2DOF simulation uses normal distributions for six uncertain parameters. The masses m1 and m2 have means of 190 and 100 kg with standard deviations of 20 and 5 kg, respectively. The stiffnesses k₁ and k₂ have means of 1000 and 2000 N/m with standard deviations of 30 and 80 N/m, respectively. The damping coefficients c₁ and c₂ have means of 50 and 230 N·s/m, each with a standard deviation of 30 N·s/m; all draws, including negative values, are retained in the reported transient results. In a normal distribution, values close to the mean occur more frequently, while values farther from the mean occur less frequently. Figure 4(c) shows the unfiltered normal input distributions used for the 2DOF Monte Carlo simulation. Figure 4(d) is included only as an illustrative uniform-distribution representation and is not used to generate the reported 2DOF response results.

Mass and stiffness determine the undamped natural frequencies, while damping controls dissipation or energy input and affects transient displacement and damped oscillation frequencies. Parameter uncertainty is held fixed within each trial, so the response variability is across assemblies rather than a random fluctuation of parameters during an event.

Figure 4 is intended to document and illustrate the input probability assumptions rather than to provide a controlled comparison of the effects of normal and uniform distributions on vibration response. The SDOF and 2DOF simulations use different mechanical models, different parameter values, and different probability distributions. As a result, differences between their response results cannot be attributed solely to the choice of probability-distribution family. A controlled distributional comparison would require the same mechanical model, parameter values, variability, dependence structure, and physical constraints to be retained while changing only the probability distribution being tested. Such a matched comparison is not performed in the present study, and no differences between the SDOF and 2DOF outputs are therefore interpreted as being caused only by uniform or normal sampling.

6. Frequency Variability and Transient Response in Reliability Screening

Natural-frequency and transient-response statistics answer different questions. Mass and stiffness determine the undamped modal frequencies; damping and the fixed initial state determine decay, phase, and displacement amplitude. A frequency distribution alone cannot determine peak displacement, settling time, or whether a performance limit is exceeded over an interval.

For the SDOF benchmark, the nominal frequency is 1.38013 rad/s and the sampled mean is 1.38036 rad/s. Their agreement concerns modal-frequency variability and does not validate the Euler time integration. The analytical audit shows that accurate frequency sampling can coexist with biased peak-response estimates.

For the 2DOF benchmark, the lower and upper frequency statistics describe variability in the assembly modes. The damped transient statistics additionally depend on the damping assumptions, including the retained negative c₁ draws. The nominal decay half-lives are specific to the mean-parameter damped model, whereas the population response bands summarize many different systems.

Peak displacement, settling time, and residual vibration are response functionals rather than undamped modal properties. A settling-time criterion must specify its displacement tolerance and the subsequent interval during which the response must stay within that tolerance; a single instantaneous percentile is insufficient.

A future screening study could compare acceptance based on frequency alone with acceptance based on frequency and a transient metric. A meaningful comparison would define the same manufactured population, dependence assumptions, acceptance yield, operating conditions, service interval, and failure event. Accepted units would retain their sampled properties while any service-age changes were modeled explicitly.

No factory-acceptance or service-reliability comparison is performed here. The thresholds 13.5 m and 16 m are illustrative benchmark values, and the model contains no calibrated degradation law, service loading, or experimental validation. The results therefore support numerical uncertainty characterization and motivate future screening research; they do not show that combined frequency-and-decay screening improves service reliability.

7. Conclusions

This study examines Monte Carlo uncertainty propagation and numerical verification in linear SDOF and 2DOF free-vibration benchmarks with 100,000 trials per model and fixed initial states. For SDOF, the undamped frequency mean is 1.38036 rad/s and sample standard deviation 0.02817 rad/s. The original Euler output retained in Figure 3 gives 67.384% exceedance at 13.5 m, whereas the analytical solution for the same samples gives 39.842%. The analytical value supersedes the original value as the model-based threshold estimate and demonstrates that numerical bias must be distinguished from Monte Carlo sampling uncertainty.

For 2DOF, the lower and upper undamped frequencies have means 1.80020 and 5.72171 rad/s and sample standard deviations 0.06309 and 0.18961 rad/s. Their coefficients of variation are 3.50% and 3.31%. At the mean parameters, damped state-matrix eigenvalues give decay half-lives of 7.88 s and 0.385 s. These damped quantities are separate from the undamped modal illustrations in Figure 5(c) & Figure 5(d). The recorded-grid either-mass exceedance at 16 m is 2.020%, conditional on the stated unfiltered population and time grid.

The 2DOF ensemble includes 4805 negative c₁ draws and therefore does not represent an entirely passive mechanical population. Central percentile bands can conceal extreme responses, and recorded-grid peaks can miss larger continuous-time maxima. A passive-system assessment would require physically admissible inputs and a new transient analysis. Frequency statistics are unaffected by damping because the undamped eigenproblem uses only mass and stiffness.

The models use different mechanical parameters and different input distribution families, so their outputs do not establish an isolated effect of normal versus uniform sampling. Initial conditions are fixed, not random, and every trial has time-invariant mechanical coefficients.

The benchmarks demonstrate the complementary roles of modal-frequency and time-domain statistics and the importance of checking numerical response calculations before interpreting threshold exceedance. Application-specific reliability conclusions require justified input distributions, performance limits, operating and degradation models, and experimental or field validation. Such conclusions are outside the present study.

Supporting Information

The SDOF and 2DOF simulation code is publicly available at https://github.com/angelasun0210-del/Monte-Carlo-Simulation-Vibration-simulation/releases/tag/v1.0.0, version v1.0.0, commit 1f24ddcc149c0805d0502c7856c5468585c8e260. Both models were executed with 100,000 trials. The repository includes software requirements, random seed settings, numerical solver settings, reproduction instructions, and reference outputs.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References