We characterize the complex, heavy-tailed probability density functions (pdfs) describing the response and its local extrema for structural systems subject to random forcing that includes extreme events. Our approach is based on recent probabilistic decomposition-synthesis (PDS) technique (Mohamad, M. A., Cousins, W., and Sapsis, T. P., 2016, “A Probabilistic Decomposition-Synthesis Method for the Quantification of Rare Events Due to Internal Instabilities,” J. Comput. Phys., 322, pp. 288–308), where we decouple rare event regimes from background fluctuations. The result of the analysis has the form of a semi-analytical approximation formula for the pdf of the response (displacement, velocity, and acceleration) and the pdf of the local extrema. For special limiting cases (lightly damped or heavily damped systems), our analysis provides fully analytical approximations. We also demonstrate how the method can be applied to high dimensional structural systems through a two-degrees-of-freedom (TDOF) example system undergoing extreme events due to intermittent forcing. The derived formulas can be evaluated with very small computational cost and are shown to accurately capture the complicated heavy-tailed and asymmetrical features in the probability distribution many standard deviations away from the mean, through comparisons with expensive Monte Carlo simulations.

## Introduction

A large class of physical systems in engineering and science can be modeled by stochastic differential equations. For many of these systems, the dominant source of uncertainty is due to the forcing which can be modeled by a stochastic process. Applications include ocean engineering systems excited by water waves (such as ship motions in large waves [14] or high speed crafts subjected to rough seas [5,6]) and rare events in structural systems (such as beam buckling [7,8], vibrations due to earthquakes [9,10], and wind loads [11,12]). For all of these cases, it is common that hidden in the otherwise predictable magnitude of the fluctuations are extreme events, i.e., abnormally large magnitude forces which lead to rare responses in the dynamics of the system (Fig. 1). Clearly, these events must be adequately taken into account for the effective quantification of the reliability properties of the system. In this work, we develop an efficient method to fully describe the probabilistic response of linear structural systems under general time-correlated random excitations containing rare and extreme events.

Fig. 1
Fig. 1

Systems with forcing having these characteristics pose significant challenges for traditional uncertainty quantification schemes. While there is a large class of methods that can accurately resolve the statistics associated with random excitations (e.g., the Fokker–Planck (FP) equation [13,14] for systems excited by white-noise and the joint response-excitation method [1518] for arbitrary stochastic excitation), these have important limitations for high dimensional systems. In addition, even for low-dimensional systems determining the part of the probability density function (pdf) associated with extreme events poses important numerical challenges. On the other hand, Gaussian closure schemes and moment equation or cumulant closure methods [19,20] either cannot “see” rare events completely or they are very expensive and require the solution of an inverse moment problem in order to determine the pdf of interest [21]. Similarly, approaches relying on polynomial-chaos expansions [22,23] have been shown to have important limitations for systems with intermittent responses [24].

For system reliability, Monte Carlo and parameter estimation based methods are another class of strategies. It is now well established that anomalous diffusion leads to Lévy motion when particle jumps have power law tail probabilities (compared to Brownian motion, where the particle jumps have finite first two moments) for which fractional calculus is the appropriate tool for such problems [25]. Consequently, Lévy stable, Pareto, and Mittag–Leffler distributions have all been utilized in methods to understand and quantify non-Gaussian behavior in fractional-order mechanics [2629], including fatigue modeling in solid mechanics [30], where it has been noted that the relatively simple description of Brownian diffusion is not adequate to capture observed non-Gaussian behavior.

Another popular approach for the study of rare event statistics in systems under intermittent forcing is to represent extreme events in the forcing as identically distributed independent impulses arriving at random times. The generalized FP equation or Kolmogorov–Feller equation is the governing equation that solves for the evolution of the response pdf under Poisson noise [14]. However, exact analytical solutions are available only for a limited number of special cases [31]. Although alternative methods such as the path integral method [3234] and the stochastic averaging method [35,36] may be applied, solving the FP or Kolmogorov–Feller equations is often very expensive even for very low dimensional systems [37,38].

Here, we consider the quantification of the probability distribution of the response of linear systems subjected to stochastic forcing containing extreme events based on recently formulated probabilistic-decomposition synthesis (PDS) method [39,40]. The approach relies on the decomposition of the statistics into a “nonextreme core,” typically Gaussian, and a heavy-tailed component. This decomposition is in full correspondence with a partition of the phase space into a “stable” region where we do not have extreme events and a region where nonlinear instabilities or external forcing lead to rare transitions with high probability. We quantify the statistics in the stable region using a Gaussian approximation approach, while the non-Gaussian distribution associated with the intermittently unstable regions of phase space is performed taking into account the nontrivial character of the dynamics (either due to an instability or external forcing). The probabilistic information in the two domains is analytically synthesized through a total probability argument.

### Central Contributions.

We begin with the simplest case of a linear, single-degree-of-freedom (SDOF) system and then formulate the method for multidegree-of-freedom systems. The main result of our work is the derivation of analytic/semi-analytic approximation formulas for the response pdf and the pdf of the local extrema of intermittently forced systems that can accurately characterize the statistics many standard deviations away from the mean. In brief, the principal contributions of this paper are as follows:

• Analytical (under certain conditions) and semi-analytical (under no restrictions) pdf expressions for the response displacement, velocity, and acceleration for single-degree-of-freedom systems under intermittent forcing.

• Semi-analytical pdf expressions for the value and the local extrema of the displacement, velocity, and acceleration for multidegree-of-freedom systems under intermittent forcing.

The systems considered in this work are linear; however, the general method is directly applicable to nonlinear structural systems (see Ref. [41] for nonlinear systems).

The proposed approach circumvents the challenges that rare and extreme events pose for traditional uncertainty quantification schemes, in particular the computational burden associated when dealing with rare events in systems. We emphasize the statistical accuracy and the computational efficiency of the presented approach, which we rigorously demonstrate through extensive comparisons with direct Monte Carlo simulations.

The paper is structured as follows: In Sec. 2, we provide a general formulation of the probabilistic decomposition-synthesis method for structural systems under intermittent forcing. Next, in Sec. 3, we apply the developed method analytically, which is possible for two limiting cases: underdamped systems with $ζ≪1$ or overdamped with $ζ≫1$, where ζ is the damping ratio. The system we consider is excited by a forcing term consisting of a background time-correlated stochastic process superimposed with a random impulse train (describing the rare and extreme component). We give a detailed derivation of the response pdf of the system (displacement, velocity, and acceleration) and compare the results with expensive Monte Carlo simulations. In Sec. 4, we slightly modify the developed formulation to derive a semi-analytical scheme for the same linear system without restriction on the damping ratio ζ, demonstrating global applicability of the proposed approach. In Sec. 5, we consider multiple-degree-of-freedom systems, and in Sec. 6, include the response pdf results for the local extrema. Concluding remarks are made in Sec. 7.

## The Probabilistic Decomposition-Synthesis Method for Intermittently Forced Structural Systems

We provide a description of the recently formulated PDS method adapted to the case of intermittently forced linear structural systems [40]. Consider the following vibrational system:
$Mẍ(t)+Dẋ(t)+Kx(t)=F(t), x(t)∈ℝn$
(1)
where M is a mass matrix, D is the damping matrix, and K is the stiffness matrix. We assume that $F(t)$ is a stochastic forcing with intermittent characteristics that can be expressed as
$F(t)=Fb(t)+Fr(t)$
(2)
The forcing consists of a background component $Fb$ of characteristic magnitude σb and a rare and extreme component $Fr$ with magnitude $σr≫σb.$ The components $Fb$ and $Fr$ may both be (weakly) stationary stochastic processes, while the sum of the two processes is in general nonstationary. This can be seen if we directly consider the sum of two (weakly) stationary processes x1 and x2, with time correlation functions $Corrx1(τ)$ and $Corrx2(τ)$, respectively. Then for the sum $z=x1+x2$, we have
$Corrz(t,τ)=Corrx1(τ)+E[x1(t)x2(t+τ)]+E[x1(t+τ)x2(t)]+Corrx2(τ)$
(3)

Therefore, the process z is stationary if and only if the cross-covariance terms $E[x1(t)x2(t+τ)]$ and $E[x1(t+τ)x2(t)]$ are functions of τ only or they are zero (i.e., x1 and x2 are not correlated).

For the case where the excitation is given in terms of realizations, i.e., time-series, one can first separate the extreme events from the stationary background by applying time-frequency analysis methods (e.g., wavelets [42]). Then, the stationary background can be approximated with a Gaussian stationary stochastic process (with properly tuned covariance function), while the rare event component can be represented with a Poisson process with properly calibrated parameters that represent the characteristics of the extreme forcing events (frequency and magnitude).

To apply the PDS method, we also decompose the response into two terms
$x(t)=xb(t)+xr(t)$
(4)

where $xb$ accounts for the background state (nonextreme) and $xr$ captures extreme responses (due to the intermittent forcing)—see Fig. 2. More precisely, $xr$ is the system response under two conditions: (1) the forcing is given by $F=Fr$ (i.e., we have an impulse) and (2) the norm of the response is greater than the background response fluctuations according to a given criterion, e.g., $‖x‖>γ$. However, as we will see in Sec. 3, other criteria may be used. These rare transitions occur when we have a forcing impulse and include a period where the system relaxes to the background state $xb$. The background component $xb$ corresponds to the response without rare events $xb=x−xr$, and in this regime, the dynamics are primarily governed by the background forcing term $Fb$.

Fig. 2
Fig. 2

We require that rare events are statistically independent from each other. In the general formulation of the PDS method, we also need to assume that rare events have negligible effects on the background state $xb$, but in the current setup this assumption is not necessary due to the linear character of the examples we consider. However, in order to apply the method for general nonlinear structural systems, we need to satisfy this condition. Also, we require ergodicity since the statistics we approximate correspond to long time averages of trajectories [40].

Next, we focus on the statistical characteristics of an individual mode $u(t)∈ℝ$ of the original system in Eq. (1). The first step of the PDS method is to quantify the conditional statistics of the rare event regime. When the system enters a rare event regime, say, at t = t0, we have an arbitrary background state $u(t0)=ub$ as the initial condition. The problem in this regime can be formulated as
$ür(t)+λu̇r(t)+kur(t)=Fr(t), with ur(t0)=ub and F=Fr for t>t0$
(5)

Under the assumption of independent rare events, we can use Eq. (5) as the starting point for analytical and numerical approximations for the statistical response during rare events.

The background component, on the other hand, can be studied through the equation
$Mẍb(t)+Dxḃ(t)+Kxb(t)=Fb(t)$
(6)

Because of the nonintermittent character of the response in this regime, it is sufficient to obtain the low-order statistics of this system. In the context of vibrations, it is reasonable to assume that $Fb(t)$ is a Gaussian process or can be adequately approximated by a Gaussian process, in which case the response statistics is readily obtained from the associated Fokker–Planck equation. Consequently, this step provides us with the statistical steady-state probability distribution for the mode of interest under the condition that the dynamics “live” in the stochastic background.

After analysis of the two regimes is complete, we can synthesize the results through a total probability argument
$f(q)=f(q|‖u‖>γ,F=Fr)︸rare eventsPr+f(q|F=Fb)︸background(1−Pr)$
(7)
where q may be any function of interest involving the response. In the last equation, $Pr$ denotes the overall rare event probability. This is defined as the probability of the response exceeding a threshold γ due to extreme event excitations
$Pr≡P(‖u‖>γ,F=Fr)=1T∫t∈T1(‖u‖>γ,F=Fr) dt$
(8)

where $1( · )$ is the indicator function. The rare event probability measures the total rare event duration.

The utility of the proposed decomposition is its flexibility and effectiveness in cheaply capturing rare responses, since we account for rare event dynamics directly by connecting their statistical properties with the original system response.

### Problem Formulation for Linear Single-Degree-of-Freedom Systems.

To demonstrate the method we begin with a simple example and consider a single-degree-of-freedom linear system
$ẍ+λẋ+kx=F(t)$
(9)

where k is the stiffness, λ is the damping, and $ζ=λ/2k$ is the damping ratio. For what follows, we adopt the standard definitions: $ωn=k$, $ωo=ωnζ2−1$, and $ωd=ωn1−ζ2$.

The forcing F(t) is a stochastic process with intermittent characteristics, which can be written as
$F(t)=Fb(t)+Fr(t)$
(10)

Here, Fb is the background forcing component that has a characteristic magnitude σb, and Fr is a rare and large amplitude forcing component that has a characteristic magnitude $σr$, which is much larger than the magnitude of the background forcing $σr≫σb$. Despite the simplicity of the system, its response may feature a significantly complicated statistical structure with heavy-tailed characteristics (large kurtosis) and asymmetry (skewness).

For concreteness, we consider a prototype system motivated from ocean engineering applications, which models a base excitation problem of a structural mode due to ocean wave impacts
$ẍ+λẋ+kx=ḧ(t)+∑i=1N(t)αi δ(t−τi), 0
(11)
Here, h(t) denotes the zero-mean background base motion term (having opposite sign from x) with a Pierson–Moskowitz spectrum (ocean wave spectrum, although any suitable spectrum may be utilized)
$Shh(ω)=q1ω5exp(−1ω4)$
(12)

where q controls the magnitude of the forcing.

The second forcing term in Eq. (11) describes rare and extreme events. In particular, we assume that this component is a random impulse train ($δ( · )$ is a unit impulse), where N(t) is a Poisson counting process that represents the number of impulses that arrive in the time interval $0, α is the impulse mean magnitude (characterizing the rare event magnitude σr), which we assume is normally distributed with mean $μα$, variance $σα2$ and independent from the state of the system. In addition, the arrival rate is constant and given by $να$ (or by the mean arrival time $Tα=1/να$ so that impulse arrival times are exponentially distributed $τ∼eTα$).

We take the impulse mean magnitude as being m-times larger than the standard deviation of the excitation velocity $ḣ(t)$
$μα=mσḣ, with m>1$
(13)

where $σḣ$ is the standard deviation of $ḣ(t)$. The factor m sets the severity of the extreme impulse events, with respect to the background amplitude level, so that $σr≫σb$.

The setup above is applicable to numerous applications in other domains (with appropriate modifications on the forcing, e.g., the background excitation spectrum shape), including structures under wind excitations, systems under seismic excitations, and vibrations of high-speed crafts and road vehicles [5,13,14].

## Analytical Probability Density Functions of Single-Degree-of-Freedom Systems for Limiting Cases

In this section, we apply the probabilistic decomposition-synthesis method to derive analytical approximations for the pdf of the displacement, velocity, and acceleration for two special limiting cases: underdamped motion with $ζ≪1$ and overdamped motion with $ζ≫1$. We first perform the analysis for the response displacement and then derive the response pdf for the velocity and acceleration.

### Background Response Probability Density Functions.

Consider the statistical response of the system to the background forcing component
$xb¨+λxḃ+kxb=ḧ(t)$
(14)
This is a linear system that is excited by a Gaussian process. The response statistics is thus also Gaussian and is fully characterized by the response spectrum. The spectral density of the displacement, velocity, and acceleration is given by
$Sxbxb(ω)=ω4Shh(ω){k−ω2+λ(jω)}2, Sẋbẋb(ω)=ω2Sxbxb(ω), Sẍbẍb(ω)=ω4Sxbxb(ω)$
(15)
Using the spectral response, we compute the corresponding variance values, which fully characterize the zero mean Gaussian response in the background regime
$σxb2=∫0∞Sxbxb(ω) dω, σẋb2=∫0∞Sẋbẋb(ω) dω, σẍb2=∫0∞Sẍbẍb(ω) dω$
(16)
Furthermore, the corresponding response envelopes are Rayleigh distributed [43]
$ub∼R(σxb), u̇b∼R(σẋb), üb∼R(σẍb)$
(17)

where $R(σ)$ denotes a Rayleigh distribution with parameter σ.

### Analytical Probability Density Functions for the Underdamped Case ζ≪1.

In the strongly underdamped case with $ζ≪1$, due to the highly oscillatory nature of the response, we focus on deriving the statistics of the local extrema. Consequently, we present results in terms of the envelope of the response statistics.

#### Rare Events Response.

To estimate the rare event response, we must take into account the nonzero background velocity of the system $ẋb$ at the moment of the impulse impact, as well as the magnitude of the impact, α. The actual value of the response xb is considered negligible. For the present case, $ζ≪1$, we have that the envelopes of the response (displacement, velocity, and acceleration) during the rare event are given by (see Appendix  A for details)
$ur(t)≃|ẋb+α|ωde−ζωnt, u̇r(t)≃|ẋb+α|e−ζωnt, ür(t)≃ωd|ẋb+α|e−ζωnt$
(18)
Above, in Eq. (18), the two contributions $ẋb$ and α in the term $ẋb+α$ are both Gaussian distributed and independent, and therefore, their sum is also Gaussian distributed
$η≡ẋb+α∼N(μα, σẋb2+σα2)$
(19)
Therefore, the distribution of the quantity $|η|$ is given by a folded normal distribution
$f|η|(n)=1σ|η|2π{ exp(−(n−μα)22σ|η|2)+exp(−(n+μα)22σ|η|2)}, 0
(20)

where $σ|η|=σẋb2+σα2$.

#### Rare Event Probability.

Next, we compute the rare event probability: the total rare event period over a fixed time interval, as defined in Eq. (8). This is done by employing an appropriate definition of an extreme event in terms of a threshold value. One possible option is to set an absolute threshold γ. However, in the current context, it is far more convenient to set this threshold relative to the value of the local maximum of the extreme event response. Specifically, the time duration τe a rare response takes to return back to the background state will be given by the duration starting from the initial impulse event time ($t0$) to the point where the response has decayed back to $ρc$ (or $100ρc%$) of its absolute maximum; here and throughout this paper, we set $ρc=0.1$. This is a value that we considered without tuning. We emphasize that the derived approximation is not sensitive to the exact value of ρc as long as this has been chosen within reasonable values.

In the current context, this means that the rare event duration τe is defined by
$ur(τe+t0)=ρcur(t0)$
(21)
We solve the above using the analytically derived response envelope equations from Sec. 3.2.1
$τe=−1ζωnlog ρc$
(22)
Due to the linear character of the system, the typical duration τe is independent of the background state or the impulse impact intensity. Using this rare event duration τe, we can compute the total rare event probability from the arrival frequency $να$ (recall $να=1/Tα$) of the Poisson process that describes the frequency of impulse events
$Pr=νατe=τe/Tα$
(23)

Note, since we naturally assume that extreme events are rare enough to be statistical independent, the above probability is smaller than 0.5.

#### Conditional Probability Density Functions for Rare Events.

We now proceed with the derivation of the pdf in the rare event regime. Consider again the response displacement during a rare event
$ur(t)∼|η|ωde−ζωnt′$
(24)
Here $t′$ is a random variable uniformly distributed between the initial impulse event time and the end time τe (from Eq. (22)) when the response has relaxed back to the background dynamics
$t′∼uniform(0, τe)$
(25)
We condition the rare event distribution as follows:
$fur(r)=∫fur||η|(r|n)f|η|(n) dn$
(26)

where we have already derived the pdf for $f|η|$ in Eq. (20). What remains is the derivation of the conditional pdf for $fur∣|η|$.

By conditioning on $|η|=n$, we find the derived distribution for the conditional pdf given by
$fur||η|(r|n)=1rζωnτe{s(r−nωde−ζωnτe)−s(r−nωd)}$
(27)

where $s( · )$ denotes the step function, which is equal to 1 when the argument is greater than or equal to 0 and 0 otherwise. Refer to Appendix  B for a detailed derivation of the above.

Using Eqs. (27), (19), and (26), we obtain the final result for the rare event distribution for the response displacement
$fur(r)=∫fur||η|(r|n)f|η|(n) dn$
(28)

$=1rζωnσ|η|τe2π∫0∞{exp(−(n−μα)22σ|η|2)+exp(−(n+μα)22σ|η|2)}×{s(r−nωde−ζωnτe)−s(r−nωd)}dn$
(29)

#### Summary of Results for the Underdamped Case

##### Displacement envelope.
Finally, combining the results of Secs. 3.1, 3.2.3, and 3.2.2 using the total probability law
$fu(r)=fub(r)(1−Pr)+fur(r)Pr$
(30)
We obtain the desired envelope distribution for the displacement of the response
$fu(r)=rσxb2exp(−r22σxb2)(1−νατe)+νατerζωnσ|η|τe2π∫0∞{ exp(−(n−μα)22σ|η|2)+exp(−(n+μα)22σ|η|2)}×{s(r−nωde−ζωnτe)−s(r−nωd)}dn$
(31)

where $τe=−log(ρc)/(ζωn)$ and $s( · )$ denotes the step function.

##### Velocity envelope.
Similar to the derivation for the displacement, we can also obtain the envelope distribution of the velocity of the system. The statistical response of the background dynamics for the velocity response was obtained in Eq. (17). Noting that from Eq. (18), $u̇r=ωdur$, the rare event pdf is modified by a constant factor and thus
$fu̇(r)=fu̇b(r)(1−Pr)+ωd−1fur(r/ωd)Pr$
(32)
The final formula for the velocity envelope pdf is then given by
$fu̇(r)=rσẋb2exp(−r22σẋb2)(1−νατe)+ναrζωnσ|η|2π∫0∞{ exp(−(n−μα)22σ|η|2)+exp(−(n+μα)22σ|η|2)}×{s(r−ne−ζωnτe)−s(r−n)}dn$
(33)
##### Acceleration envelope.
Finally, we also obtain the envelope distribution of the acceleration. Noting that from Eq. (18), $ür=ωd2ur$, the rare event pdf for acceleration is also modified by a constant factor
$fü(r)=füb(r)(1−Pr)+ωd−2fur(r/ωd2)Pr$
(34)
The final formula for the acceleration envelope pdf is then
$fü(r)=rσẍb2exp(−r22σẍb2)(1−νατe)+ναrζωnσ|η|2π∫0∞{ exp(−(n−μα)22σ|η|2)+exp(−(n+μα)22σ|η|2)}×{s(r−nωde−ζωnτe)−s(r−nωd)}dn$
(35)

#### Comparison With Monte Carlo Simulations.

For Monte Carlo simulations, the excitation time series is generated by superimposing the background and rare event components together. The background excitation, described by a stationary stochastic process with a Pierson–Moskowitz spectrum (Eq. (12)), is simulated through a superposition of cosines over a range of frequencies with corresponding amplitudes and uniformly distributed random phases. The intermittent component is the random impulse train, and each impact is introduced as a velocity jump with a given magnitude at the point of the impulse impact. For each of the comparisons, we generate ten realizations of the excitation time series, each with a train of 100 impulses. Once each time series for the excitation is computed, the governing ordinary differential equations are solved using a 4th/5th-order Runge–Kutta method (we carefully account for the modifications in the momentum that an impulse imparts by integrating up to each impulse time and modifying the initial conditions that the impulse imparts before integrating the system to the next impulse time). We ensured that the time series length and ensemble size lead to converged response statistics (displacement, velocity, and acceleration).

Furthermore, we utilize a shifted Pierson–Moskowitz spectrum $Shh(ω−1)$ in order to avoid resonance conditions. The other parameters and resulting system statistical properties are given in Table 1. As seen in Fig. 3, the analytical approximations compare well with the Monte Carlo simulations many standard deviations away from the zero mean. The results are robust to different parameters as long as the assumption of independent (nonoverlapping) extreme events is satisfied.

Fig. 3
Fig. 3
Table 1

Parameters and relevant statistical quantities for SDOF system 1

 λ 0.01 k 1 $Tα$ 5000 ζ 0.005 ωn 1 ωd 1 $μα=7×σh˙$ 0.1 q $1.582×10−4$ $σα=σh˙$ 0.0143 σh 0.0063 $σx˙b$ 0.0179 $σxb$ 0.0082 $σ|η|$ 0.0229 $Pr$ 0.0647
 λ 0.01 k 1 $Tα$ 5000 ζ 0.005 ωn 1 ωd 1 $μα=7×σh˙$ 0.1 q $1.582×10−4$ $σα=σh˙$ 0.0143 σh 0.0063 $σx˙b$ 0.0179 $σxb$ 0.0082 $σ|η|$ 0.0229 $Pr$ 0.0647
In Table 2, we include two quantitative error measures between our analytical approximation and Monte Carlo simulations. They include the Kullback–Lieber divergence (i.e., the relative entropy) and the maximum absolute error, computed using
$EKL(x)=∫fx,MC(q)log(fx,MC(q)fx(q)) dq, Eabs(x)=maxq|fx(q)−fx,MC(q)|fx,MC(q)$
(36)
respectively.
Table 2

Error measures for SDOF system 1

 $EKL(u)$ $1.1×10−4$ $Eabs(u)$ 0.013 $EKL(u)$ $5.3×10−3$ $Eabs(u)$ 0.024 $EKL(u)$ $1.0×10−3$ $Eabs(u)$ 0.065
 $EKL(u)$ $1.1×10−4$ $Eabs(u)$ 0.013 $EKL(u)$ $5.3×10−3$ $Eabs(u)$ 0.024 $EKL(u)$ $1.0×10−3$ $Eabs(u)$ 0.065

Some of the discrepancies observed between the Monte Carlo simulations and analytical results can be attributed to the envelope approximation used to derive the conditionally rare event pdf. Indeed, these discrepancies are significantly reduced (results not shown) if we utilize the semi-analytical method presented in Sec. 4, where no simplifications are made for the form of the rare event response.

### Analytical Probability Density Functions for the Overdamped Case ζ≫1.

In Sec. 3.2, we derived the analytical response pdf under the underdamped assumption $ζ≪1$. Here, we briefly summarize the results for the response pdf in the overdamped case with $ζ≫1$. The derivation follows the exact same steps as in Secs. 3.1 and 3.2, but instead uses the corresponding formulas for the rare event response in overdamped conditions $ζ≫1$ (see Appendix  A). An important difference in the overdamped scenario is that the system does not exhibit oscillatory response, hence we directly operate on the response instead of the envelope of the response.

#### Displacement.

The total probability law becomes
$fx(r)=fxb(r)(1−Pr,dis)+fxr(r)Pr,dis$
(37)
and we obtain the following pdf for the displacement of the system:
$fx(r)=1σxb2πexp(−r22σxb2)(1−νατe,dis)+νατe,disr(ζωn−ωo)ση2π(τe,dis−τs)∫−∞∞ exp (−(n−μα)22ση2)×{s(r−n2ωoe−(ζωn−ωo)τe,dis)−s(r−n2ωo)}dn$
(38)

where $τe,dis=(π/2ωo)−log ρc/(ζωn+ωo)$.

#### Velocity.

Similarly, we derive the total probability law for the response velocity
$fẋ(r)=fẋb(r)(1−Pr,vel)+fẋr(r)Pr,vel$
(39)
The final result for the velocity pdf is given by
$fẋ(r)=1σẋb2πexp(−r22σẋb2)(1−νατe,vel)+νατe,velr(ζωn+ωo)ση2πτe,vel∫−∞∞ exp (−(n−μα)22ση2)×{s(r−ne−(ζωn+ωo)τe,vel)−s(r−n)}dn$
(40)

where $τe,vel=−(1/ζωn+ωo)log ρc$.

#### Acceleration.

The total probability law for the response acceleration is
$fẍ(r)=fẍb(r)(1−Pr,acc)+fẍr(r)Pr,acc$
(41)
and this gives the following result for the acceleration pdf:
$fẍ(r)=1σẍb2πexp(−r22σẍb2)(1−νατe,acc)+νατe,accr(ζωn+ωo)ση2πτe,acc∫−∞∞ exp (−(n−μα)22ση2)×{s(r−n(ζωn+ωo)e−(ζωn+ωo)τe,acc)−s(r−n(ζωn+ωo))}dn$
(42)

where $τe,acc=−log ρc/(ζω+ωo)$.

Note that here, we do not have the simple scaling as in the underdamped case for the conditionally rare component of the pdf.

#### Comparison With Monte Carlo Simulations.

We confirm the accuracy of the analytical results given in Eqs. (38), (40), and (42) for the strongly overdamped case through comparison with direct Monte Carlo simulations. The parameters and resulting statistical quantities of the system are given in Table 3. The analytical expressions for the response show excellent agreement compared to the Monte Carlo simulations, see Fig. 4 and Table 4.

Fig. 4
Fig. 4
Table 3

Parameters and relevant statistical quantities for SDOF system 2

 λ 6 k 1 $Tα$ 1000 ζ 3 ωn 1 ωd 2.828 $μα=7×σh˙$ 0.1 q $1.582×10−4$ $σα=σh˙$ 0.0143 σh 0.0063 $σx˙b$ 0.0056 $σxb$ 0.0022 $ση$ 0.0154 $Pr,dis$ 0.0140 $Pr,vel$ 0.0004 $Pr,acc$ 0.0004
 λ 6 k 1 $Tα$ 1000 ζ 3 ωn 1 ωd 2.828 $μα=7×σh˙$ 0.1 q $1.582×10−4$ $σα=σh˙$ 0.0143 σh 0.0063 $σx˙b$ 0.0056 $σxb$ 0.0022 $ση$ 0.0154 $Pr,dis$ 0.0140 $Pr,vel$ 0.0004 $Pr,acc$ 0.0004
Table 4

Error measures for SDOF system 2

 $EKL(x)$ $2.6×10−4$ $Eabs(x)$ 0.016 $EKL(x˙)$ $6.5×10−5$ $Eabs(x˙)$ 0.024 $EKL(x¨)$ $1.9×10−3$ $Eabs(x¨)$ 0.028
 $EKL(x)$ $2.6×10−4$ $Eabs(x)$ 0.016 $EKL(x˙)$ $6.5×10−5$ $Eabs(x˙)$ 0.024 $EKL(x¨)$ $1.9×10−3$ $Eabs(x¨)$ 0.028

## Semi-Analytical Probability Density Functions of the Response of Single-Degree-of-Freedom Systems

We now formulate a semi-analytical approach to quantify the response pdf for arbitrary system parameters (including the severely underdamped or overdamped cases considered previously). The approach here adapts the numerical scheme described in Ref. [1] for systems undergoing internal instabilities.

Whereas for the limiting cases that were previously studied, where (time series) knowledge of the system trajectory ($xr(t)$ or $ur(t)$) could be translated to information about the corresponding pdf ($fxr$ or $fur$), this is not always possible. In addition, nonlinear structural systems, in general, have analytical expressions for the rare event transitions. For these cases, we can compute the rare event statistics by numerically approximating the corresponding histogram, using either analytical or numerically generated trajectories in the rare event regime.

### Numerical Computation of Rare Events Statistics.

Consider the same SDOF system introduced in Sec. 3. Recall that we quantify the response pdf by the PDS method via a total probability law argument
$fx(r)=fxb(r)(1−Pr)+fxr(r)Pr$
(43)

In Sec. 3, the derivation consisted of estimating all three unknown quantities: the background density $fxb$, the rare event density $fxr$, and the rare event probability $Pr$ analytically. However, in the semi-analytical scheme, we will obtain the rare event density $fxr$ and the rare event probability $Pr$ by taking a histogram of the numerically simulated analytical form of the rare response. The background density $fxb$ will still be obtained analytically as in Sec. 3.1.

Recall that the rare event density is given by
$fxr(r)=∫fxr|η(r|n)fη(n) dn$
(44)
where $fη(n)$ is known analytically (Eq. (19)). It is the conditional pdf $fxr∣η(r∣n)$ that we estimate by a histogram
$fxr|η(r|n)=Hist{xr|η(t|n)}, t=[0,τe,dis]$
(45)
where we use the analytical solution of the oscillator with nonzero initial velocity n
$xr|η(t|n)=n2ωo(e−(ζωn−ωo)t−e−(ζωn+ωo)t)$
(46)
The histogram is taken from t = 0 (the beginning of the rare event) until the end of the rare event at $t=τe$. The conditional density of rare event response for the velocity and acceleration are, similarly, given by
$fẋr|η(r|n)=Hist{ẋr|η(t|n)}, t=[0,τe,vel]$
(47)

$fẍr|η(r|n)=Hist{ẍr|η(t|n)}, t=[0,τe,acc]$
(48)

### Numerical Estimation of the Rare Events Probability.

In order to compute the histogram of a rare impulse event, the duration of a rare response needs to be obtained numerically. Recall that we have defined the duration of a rare response by
$xr(τe)=ρc max{|xr|}$
(49)
where $ρc=0.1$. In the numerical computation of τe, the absolute maximum of the response is obtained from the numerical simulations of the trajectories. Once the rare event duration has been determined, we compute the total rare event probability, just as before, from
$Pr=νατe=τe/Tα$
(50)
which is independent of the conditional background magnitude value. The above procedure is applied for the rare event response displacement $τe,dis$, velocity $τe,vel$, and acceleration $τe,acc$.

### Semi-Analytical Probability Density Functions.

With the modifications just described, we can readily compute the response pdf using this semi-analytical strategy. For the response displacement, this gives the formula
$fx(r)=1−νατe,disσxb2πexp(−r22σxb2)+νατe,dis∫0∞Hist{xr|η(t|n)}fη(n) dn$
(51)

The corresponding pdf for the velocity $fẋ$ and acceleration $fẍ$ can be computed with the same formula but with the appropriate variance for the Gaussian core ($σẋb$ or $σẍb$), rare event duration ($τe,vel$ or $τe,acc$), and histograms ($ẋr∣η(t∣n)$ or $ẍr∣η(t∣n)$).

#### Comparison With Monte Carlo Simulations.

For illustration, a SDOF configuration is considered with critical damping ratio, ζ = 1. The detailed parameters and relevant statistical quantities of the system are given in Table 5. This is a regime where the analytical results derived in Sec. 3 are not applicable. Even for this ζ value, the semi-analytical pdf for the response shows excellent agreement with direct simulations (Fig. 5 and Table 6). We emphasize that the computational cost of the semi-analytical scheme is comparable with that of the analytical approximations (order of seconds) and both are significantly lower than the cost of Monte Carlo simulation (order of hours).

Fig. 5
Fig. 5
Table 5

Parameters and relevant statistical quantities for SDOF system 3

 λ 2 k 1 $Tα$ 400 ζ 1 ωn 1 ωd 0 $μα=7×σh˙$ 0.1 q $1.582×10−4$ $σα=σh˙$ 0.0143 σh 0.0063 $σx˙b$ 0.012 $σxb$ 0.0052 $σα=ση$ 0.0187 $Pr,dis$ 0.0122 $Pr,vel$ 0.0075 $Pr,acc$ 0.0032
 λ 2 k 1 $Tα$ 400 ζ 1 ωn 1 ωd 0 $μα=7×σh˙$ 0.1 q $1.582×10−4$ $σα=σh˙$ 0.0143 σh 0.0063 $σx˙b$ 0.012 $σxb$ 0.0052 $σα=ση$ 0.0187 $Pr,dis$ 0.0122 $Pr,vel$ 0.0075 $Pr,acc$ 0.0032
Table 6

Error measures for SDOF system 3

 $EKL(x)$ $2.8×10−4$ $Eabs(x)$ 0.018 $EKL(x˙)$ $2.0×10−4$ $Eabs(x˙)$ 0.023 $EKL(x¨)$ $7.2×10−3$ $Eabs(x¨)$ 0.048
 $EKL(x)$ $2.8×10−4$ $Eabs(x)$ 0.018 $EKL(x˙)$ $2.0×10−4$ $Eabs(x˙)$ 0.023 $EKL(x¨)$ $7.2×10−3$ $Eabs(x¨)$ 0.048

## Semi-Analytical Probability Density Functions for the Response of Multidegree-of-Freedom Systems

An important advantage of the semi-analytical scheme is the straightforward applicability of the algorithm to multidegree-of-freedom systems. In this section, we demonstrate how the extension can be made for a two-degrees-of-freedom (TDOF) linear system.

Consider the system (see Fig. 6)
$mẍ+λẋ+kx+λa(ẋ−ẏ)+ka(x−y)=F(t)$
(52)

$maÿ+λa(ẏ−ẋ)+ka(y−x)=0$
(53)

where the stochastic forcing $F(t)=Fb(t)+Fr(t)$ is applied to the first mass (mass m) and x, y are displacements of the two masses. As before, $Fb(t)$ is the background component and $Fr(t)=∑i=1N(t)αi δ(t−τi)$ is the rare event component.

Fig. 6
Fig. 6

The background statistics are obtained by analyzing the response spectrum of the TDOF system subject to the background excitation component. Details of this derivation are provided in Appendix  C. The histograms for the rare event transitions can be computed through standard analytical expressions that can be derived for linear systems, under the set of initial conditions: $x0=y0=ẏ0=0$ and $ẋ0=n.$

Once the impulse response has been obtained, we numerically quantify the rare event distribution, as well as the rare event duration, using the semi-analytical decomposition in Sec. 4. The pdf is then given by
$fz(r)=1−νατezσzb2πexp(−r22σzb2)+νατez∫0∞Hist{zr|η(t|n)}fη(n) dn$
(54)

where z can be either of the degrees-of-freedom (x or y) or the corresponding velocities or accelerations, while $τez$ is the typical duration of the rare events, which is estimated numerically. For each desired quantity z, the corresponding variance under background forcing, temporal duration of the rare events, and the rare event histogram need to be utilized.

Results are presented for the pdf of the displacement, velocity, and acceleration of each degree-of-freedom, see Fig. 7. These compare favorably with the direct Monte Carlo simulations. The parameters and resulting statistical quantities of the system are provided in Table 7 with error measures in Table 8. Further numerical simulations (not presented) demonstrated strong robustness of the approach to different system parameters.

Fig. 7
Fig. 7
Table 7

Parameters and relevant statistical quantities for the TDOF system

 m 1 ma 1 λ 0.01 k 1 λa 1 ka 0.1 $Tα$ 1000 $ση$ 0.0199 $μα$ 0.1 q $1.582×10−4$ $σα$ 0.0143 $σFb$ 0.0351 $Pr,disx$ 0.0177 $Pr,disy$ 0.0190 $Pr,velx$ 0.0098 $Pr,vely$ 0.0209 $Pr,accx$ 0.0066 $Pr,accy$ 0.0082
 m 1 ma 1 λ 0.01 k 1 λa 1 ka 0.1 $Tα$ 1000 $ση$ 0.0199 $μα$ 0.1 q $1.582×10−4$ $σα$ 0.0143 $σFb$ 0.0351 $Pr,disx$ 0.0177 $Pr,disy$ 0.0190 $Pr,velx$ 0.0098 $Pr,vely$ 0.0209 $Pr,accx$ 0.0066 $Pr,accy$ 0.0082
Table 8

Error measures for TDOF system

 $EKL(x)$ $2.5×10−4$ $Eabs(x)$ 0.033 $EKL(x˙)$ $3.6×10−4$ $Eabs(x˙)$ 0.019 $EKL(x¨)$ $2.9×10−3$ $Eabs(x¨)$ 0.014 $EKL(y)$ $7.4×10−4$ $Eabs(y)$ 0.033 $EKL(y˙)$ $2.6×10−4$ $Eabs(y˙)$ 0.031 $EKL(y¨)$ $3.2×10−3$ $Eabs(y¨)$ 0.019
 $EKL(x)$ $2.5×10−4$ $Eabs(x)$ 0.033 $EKL(x˙)$ $3.6×10−4$ $Eabs(x˙)$ 0.019 $EKL(x¨)$ $2.9×10−3$ $Eabs(x¨)$ 0.014 $EKL(y)$ $7.4×10−4$ $Eabs(y)$ 0.033 $EKL(y˙)$ $2.6×10−4$ $Eabs(y˙)$ 0.031 $EKL(y¨)$ $3.2×10−3$ $Eabs(y¨)$ 0.019

## Semi-Analytical Probability Density Functions of Local Maxima

It is straightforward to extend the semi-analytical framework for other quantities of interest, such as the local extrema (maxima and minima) of the response, which is of interest in reliability applications. In this case, the numerically simulated analytical trajectories of rare responses can be used to compute the numerical histogram for the local extrema instead of the full response. For the conditionally nonextreme component, we can use known results from the theory of stationary Gaussian stochastic processes to describe the corresponding background pdf analytically.

### Distribution of Local Maxima Under Background Excitation.

For a stationary Gaussian process with arbitrary spectral bandwidth ϵ, the probability density function of positive extrema (maxima) is given by [44,45]
$fm+(ξ)=ϵ2πe−ξ2/2ϵ2+1−ϵ2ξe−ξ2/2Φ(1−ϵ2ϵξ), −∞≤ξ≤∞$
(55)
where $ξ=x/μ0$, x is the magnitude of the maxima, spectral bandwidth $ϵ=(1−μ22/(μ0μ4))$, and $Φ(x)=(1/2π)∫−∞xe−u2/2du$ is the standard normal cumulative distribution function. The spectral moments for the background response displacement xb are also defined as
$μn=∫0∞ωnSxbxb(ω) dω$
(56)

We note that that in the limit of an infinitesimally narrowbanded signal $(ϵ=0)$, the pdf converges to a Rayleigh distribution. On the other hand, for an infinitely broadbanded signal $(ϵ=1)$, the distribution converges to the Gaussian pdf. For a signal with in between spectral bandwidth $(0≤ϵ≤1)$, the pdf has a blended structure with the form in Eq. (55).

Considering the asymmetric nature of the intermittent excitation, we need to consider both positive and negative extrema of rare events. While for the background response, the pdf is analytically given by
$fx̂b(x)=12μ0{fm+(xμ0)+fm+(−xμ0)}, −∞≤x≤∞$
(57)

where the $x̂$ notation denotes the local extrema of x. Similar expressions can be obtained for the velocity and acceleration extrema.

### Statistics of Local Extrema During Rare Transitions.

The conditional pdf $fx̂r∣η$ for the local extrema can be numerically estimated through the histogram
$fx̂r|η(r|n)=Hist{M(xr|η(t|n))}, t=[0,τe,dis]$
(58)

where $M( · )$ is the operator that provides all the positive/negative extrema. The positive/negative extrema are defined as the points where the derivative of the signal is zero.

### Semi-Analytical Probability Density Function for Local Extrema.

The last step consists of applying the decomposition of the pdf. This takes the form
$fx̂(r)=(1−νατe,dis)fx̂b(r)+νατe,dis∫0∞Hist{M(xr|η(t|n))}fη(n) dn$
(59)

The same decomposition can be used for the velocity and acceleration local maxima. We compare the semi-analytical decomposition with Monte Carlo simulations. In Fig. 8, we present results for the two-degrees-of-freedom system. The pdf is shown for local extrema of the displacement, velocity, and acceleration for each degree-of-freedom. We emphasize the nontrivial structure of the pdf and especially the tails. We observe from these comparisons that the semi-analytical scheme accurately estimates both the tails of the response and the non-Gaussian structure of core of the distribution.

Fig. 8
Fig. 8

## Summary and Conclusions

We have formulated a robust approximation method to quantify the probabilistic response of structural systems subjected to stochastic excitation containing intermittent components. The foundation of our approach is the recently developed probabilistic decomposition-synthesis method for the quantification of rare events due to internal instabilities to the problem where extreme responses are triggered by external forcing. The intermittent forcing is represented as a background component, modeled by a colored stochastic processes with energy distributed across a range of frequencies, and additionally a rare and extreme component that can be represented by impulses that are Poisson distributed with large interarrival time. Owing to the nature of the forcing, even the probabilistic response of a linear system can be highly complex with asymmetry and complex non-Gaussian tail behavior, which would be expected if the forcing did not consist of intermittently extreme events.

The main result of this work is the derivation of analytical and semi-analytical expressions for the pdf of the response and its local extrema for structural systems (including the displacement, velocity, and acceleration pdf). These expressions decompose the pdf into a statistical core, capturing the response under background forcing, as well as a heavy-tailed component associated with extreme transitions due to rare impulse impacts. We performed a thorough analysis for linear SDOF systems under various system parameters and also derived analytical formulas for two special parameter ranges (lightly damped or heavily damped systems). The general semi-analytical decomposition is applicable for arbitrary system parameters and we have demonstrated its validity through comprehensive comparisons with Monte Carlo simulations. The general framework is also directly applicable to multidegree-of-freedom systems, as well as systems with nonlinearities. We demonstrated applicability to a 2DOF linear system consisting of two coupled masses, where the first mass is excited. Modifications of the method to compute statistics of local extrema have also been presented.

The developed approach allows for computation of the response pdf of structural systems many orders of magnitude faster than a direct Monte Carlo simulation, which is currently the only reliable tool for such computations. The rapid evaluation of response pdfs for systems excited by extreme forcing by the proposed method paves the way for enabling robust design of structural systems subjected to extreme events of a stochastic nature. In such cases, it is usually not feasible to run Monte Carlo simulations for various parameter sets during the design process owing to the computational costs associated with low probability rare events, let alone perform parametric optimization. In different work [41], we demonstrate application of the method to design and optimization in ocean engineering applications and show that the approach is well suited to predict optimal extreme-event mitigating designs and for system reliability.

## Acknowledgment

H. K. J. and M. A. M. have been supported through the first and third grants as graduate students. We are also grateful to the Samsung Scholarship Program for support of H. K. J. as well as the MIT Energy Initiative for support under the grant Nonlinear Energy Harvesting From Broad-Band Vibrational Sources By Mimicking Turbulent Energy Transfer Mechanisms.

## Funding Data

• Air Force Office of Scientific Research (Grant No. FA9550-16-1-0231).

• MITEI (Nonlinear Energy Harvesting From Broad-Band Vibrational Sources By Mimicking Turbulent Energy Transfer Mechanisms).

• Office of Naval Research (Grant Nos. N00014-14-1-0520 and N00014-15-1-2381).

### Appendix A: Impulse Response of Single-Degree-of-Freedom Systems

The response of the system
$ẍr(t)+λẋr(t)+kxr(t)=0$
(A1)
under an impulse α at an arbitrary time t0, say $t0=0$, and with zero initial value but nonzero initial velocity ($(xr,ẋr)=(0,ẋr0)$ at $t=0−$) are given by the following equations under the two limiting cases of damping:
Severely underdamped case $ζ≪1$: With the approximation of $ζ≪1$ (or $ωd≈ωn$), we can simplify responses as
$xr(t)=α+ẋr0ωde−ζωnt sin ωdt$
(A2)

$ẋr(t)=(α+ẋr0)e−ζωnt cos ωdt$
(A3)

$ẍr(t)=−(α+ẋr0)ωde−ζωnt sin ωdt$
(A4)
Severely overdamped case $ζ≫1$: Similarly, with the approximation of $ζ≫1$ (or $ωo≈ζωn$), we can simplify responses as
$xr(t)=(α+ẋr0)2ωoe−(ζωn−ωo)t$
(A5)

$ẋr(t)=(α+ẋr0)e−(ζωn+ωo)t$
(A6)

$ẍr(t)=−(ζωn+ωo)(α+ẋr0)e−(ζωn+ωo)t$
(A7)

### Appendix B: Probability Density for an Arbitrarily Exponentially Decaying Function

Consider an arbitrary time series in the following form:
$x(t)=Ae−αt, where t∼uniform(τ1,τ2)$
(B1)
where $A,α>0$ and $τ1<τ2$. The cumulative distribution function of x(t) is
$Fx(x)=P(Ae−αt
(B2)

$=P(t>1αlog(A/x))$
(B3)

$=1−P(t<1αlog(A/x))$
(B4)

$=1−∫−∞1αlog(A/x)fT(t) dt$
(B5)
where $fT(t)$ is expressed using the step function $s( · )$ as
$fT(t)=1τ2−τ1{s(t−τ1)−s(t−τ2)}, τ1<τ2$
(B6)
The pdf of the response x(t) can then be derived by differentiation
$fx(x)=ddxFx(x)$
(B7)

$=1αxfT(1αlog(Ax))$
(B8)

$=1αx(τ2−τ1){s(x−Ae−ατ2)−s(x−Ae−ατ1)}$
(B9)
We utilize the above formula for deriving analytical pdfs. Note that the step function with respect to x in the above has been derived using
$τ1
(B10)

$−ατ2<−αt<−ατ1$
(B11)

$Ae−ατ2
(B12)

### Appendix C: Background Response for Two-Degrees-of-Freedom System

Consider the statistical response of the system to the background forcing component
$mẍb+λẋb+kxb+λa(ẋb−ẏb)+ka(xb−yb)=Fb(t)$
(C1)

$maÿb+λa(ẏb−ẋb)+ka(yb−xb)=0$
The spectral densities are given by
$Sxbxb(ω)=ω4SFb(ω){A(ω)−B(ω)2C(ω)}{A(−ω)−B(−ω)2C(−ω)}$
(C2)

$Sẋbẋb(ω)=ω2Sxbxb(ω)$
(C3)

$Sẍbẍb(ω)=ω4Sxbxb(ω)$
(C4)

$Sybyb(ω)=ω4SFb(ω){A(ω)C(ω)B(ω)−B(ω)}{A(−ω)C(−ω)B(−ω)−B(−ω)}$
(C5)

$Sẏbẏb(ω)=ω2Sybyb(ω)$
(C6)

$Sÿbÿb(ω)=ω4Sybyb(ω)$
(C7)
where $SFb(ω)$ is the spectral density of $Fb(t)$, and
$A(ω)=(λa+λ)(jω)+(ka+k)−mω2$
(C8)

$B(ω)=λa(jω)+ka$
(C9)

$C(ω)=λa(jω)+ka−maω2$
(C10)
Thus, we can obtain the following conditionally background variances:
$σxb2=∫0∞Sxbxb(ω)dω, σẋb2=∫0∞Sẋbẋb(ω)dω, σẍb2=∫0∞Sẍbẍb(ω)dω$
(C11)

$σyb2=∫0∞Svbvb(ω)dω, σẏb2=∫0∞Sẏbẏb(ω)dω, σÿb2=∫0∞Sÿbÿb(ω)dω$
(C12)

## References

References
1.
,
M. A.
, and
Sapsis
,
T. P.
,
2016
, “
Probabilistic Response and Rare Events in Mathieu's Equation Under Correlated Parametric Excitation
,”
Ocean Eng.
,
120
, pp.
289
297
.
2.
Belenky
,
V. L.
, and
Sevastianov
,
N. B.
,
2007
,
Stability and Safety of Ships: Risk of Capsizing
,
The Society of Naval Architects and Marine Engineers
, Alexandria, VA.
3.
Cousins
,
W.
, and
Sapsis
,
T. P.
,
2016
, “
Reduced Order Precursors of Rare Events in Unidirectional Nonlinear Water Waves
,”
J. Fluid Mech.
,
790
, pp.
368
388
.
4.
Cousins
,
W.
, and
Sapsis
,
T. P.
,
2014
, “
Quantification and Prediction of Extreme Events in a One-Dimensional Nonlinear Dispersive Wave Model
,”
Physica D
,
280–281
, pp.
48
58
.
5.
Riley
,
M. R.
,
Coats
,
T.
,
Haupt
,
K.
, and
Jacobson
,
D.
,
2011
, “
Ride Severity Index a New Approach to Quantifying the Comparison of Acceleration Responses of High-Speed Craft
,” 11th International Conference on Fast Sea Transportation
, Honolulu, HI, Sept. 26–29, pp.
693
699
.
6.
Riley
,
M. R.
, and
Coats
,
T. W.
,
2012
, “
A Simplified Approach for Analyzing Accelerations Induced by Wave-Impacts in High-Speed Planing Craft
,”
Third Chesapeake Power Boat Symposium
, Annapolis, MD, June 15–16, pp.
14
15
.
7.
Abou-Rayan
,
A.
, and
Nayfeh
,
A.
,
1993
, “
Stochastic Response of a Buckled Beam to External and Parametric Random Excitations
,”
AIAA
Paper No. 93-1425-CP.
8.
Lin
,
Y. K.
, and
Cai
,
C. Q.
,
1995
,
Probabilistic Structural Dynamics
,
McGraw-Hill
, New York.
9.
Lin
,
Y. K.
,
1963
, “
Application of Nonstationary Shot Noise in the Study of System Response to a Class of Nonstationary Excitations
,”
ASME J. Appl. Mech.
,
30
(4), pp.
555
558
.
10.
Branstetter
,
L. J.
,
Jeong
,
G. D.
,
Yao
,
J. T.
,
Wen
,
Y.
, and
Lin
,
Y.
,
1988
, “
Mathematical Modelling of Structural Behaviour During Earthquakes
,”
Probab. Eng. Mech.
,
3
(
3
), pp.
130
145
.
11.
Lin
,
Y.
,
1996
, “
Stochastic Stability of Wind-Excited Long-Span Bridges
,”
Probab. Eng. Mech.
,
11
(
4
), pp.
257
261
.
12.
Spence
,
S.
, and
Gioffrè
,
M.
,
2012
, “
Large Scale Reliability-Based Design Optimization of Wind Excited Tall Buildings
,”
Probab. Eng. Mech.
,
28
, pp.
206
215
.
13.
Soong
,
T. T.
, and
Grigoriu
,
M.
,
1993
,
Random Vibration of Mechanical and Structural Systems
,
Prentice Hall
,
Englewood Cliffs, NJ
, p.
14690
.
14.
Sobczyk
,
K.
,
2001
,
Stochastic Differential Equations: With Applications to Physics and Engineering
, Vol.
40
,
Springer
, Berlin.
15.
Sapsis
,
T. P.
, and
Athanassoulis
,
G. A.
,
2008
, “
New Partial Differential Equations Governing the Joint, Response–Excitation, Probability Distributions of Nonlinear Systems, Under General Stochastic Excitation
,”
Probab. Eng. Mech.
,
23
(
2–3
), pp.
289
306
.
16.
Venturi
,
D.
,
Sapsis
,
T. P.
,
Cho
,
H.
, and
,
G. E.
,
2012
, “
A Computable Evolution Equation for the Joint Response-Excitation Probability Density Function of Stochastic Dynamical Systems
,”
Proc. R. Soc. A
,
468
(
2139
), p.
759
.
17.
Joo
,
H. K.
, and
Sapsis
,
T. P.
,
2016
, “
A Moment-Equation-Copula-Closure Method for Nonlinear Vibrational Systems Subjected to Correlated Noise
,”
Probab. Eng. Mech.
,
46
, pp. 120–132.
18.
Athanassoulis
,
G.
,
Tsantili
,
I.
, and
Kapelonis
,
Z.
,
2016
, “
Beyond the Markovian Assumption: Response-Excitation Probabilistic Solution to Random Nonlinear Differential Equations in the Long Time
,”
Proc. R. Soc. A
,
471
, p.
201505
.
19.
Beran
,
M.
,
1968
,
Statistical Continuum Theories
,
Interscience Publishers
, Geneva, Switzerland.
20.
Wu
,
W. F.
, and
Lin
,
Y. K.
,
1984
, “
Cumulant-Neglect Closure for Non-Linear Oscillators Under Random Parametric and External Excitations
,”
Int. J. Non Linear Mech.
,
19
(
4
), pp.
349
362
.
21.
Athanassoulis
,
G. A.
, and
,
P. N.
,
2002
, “
The Truncated Hausdorff Moment Problem Solved by Using Kernel Density Functions
,”
Probab. Eng. Mech.
,
17
(
3
), pp.
273
291
.
22.
Xiu
,
D.
, and
,
G.
,
2002
, “
The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations
,”
SIAM J. Sci. Comput.
,
24
(
2
), pp.
619
644
.
23.
Xiu
,
D.
, and
,
G.
,
2003
, “
Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos
,”
J. Comput. Phys
,
187
(
1
), pp.
137
167
.
24.
Majda
,
A. J.
, and
Branicki
,
M.
,
2012
, “
Lessons in Uncertainty Quantification for Turbulent Dynamical Systems
,”
Discrete Contin. Dyn. Syst.
,
32
(
9
), pp.
3133
3221
.
25.
Meerschaert
,
M. M.
, 2016,
Fractional Calculus, Anomalous Diffusion, and Probability
,
World Scientific
, Singapore, pp.
265
284
.
26.
Klafter
,
J.
,
Lim
,
S. C.
, and
Metzler
,
R.
, eds., 2016,
Fractional Dynamics
,
World Scientific
, Singapore.
27.
Aban
,
I. B.
,
Meerschaert
,
M. M.
, and
Panorska
,
A. K.
,
2006
, “
Parameter Estimation for the Truncated Pareto Distribution
,”
J. Am. Stat. Assoc.
,
101
(
473
), pp.
270
277
.
28.
Liang
,
Y.
, and
Chen
,
W.
,
2015
, “
A Relative Entropy Method to Measure Non-Exponential Random Data
,”
Phys. Lett. A
,
379
(
3
), pp.
95
99
.
29.
Liang
,
Y.
, and
Chen
,
W.
,
2015
, “
Reliability Analysis for Sluice Gate Anti-Sliding Stability Using Lévy Stable Distributions
,”
Signal Process.
,
107
, pp.
425
432
.
30.
Liang
,
Y.
, and
Chen
,
W.
,
2016
, “
A Regularized Miner's Rule for Fatigue Reliability Analysis With Mittag-Leffler Statistics
,”
Int. J. Damage Mech.
,
25
(
5
), pp.
691
704
.
31.
Vasta
,
M.
,
1995
, “
Exact Stationary Solution for a Class of Non-Linear Systems Driven by a Non-Normal Delta-Correlated Process
,”
Int. J. Non-Linear Mech.
,
30
(
4
), pp.
407
418
.
32.
Köylüoglu
,
H. U.
,
Nielsen
,
S. R. K.
, and
Iwankiewicz
,
R.
,
1995
, “
Response and Reliability of Poisson-Driven Systems by Path Integration
,”
J. Eng. Mech.
,
121
(
1
) pp.
117
130
.
33.
Iwankiewicz
,
R.
, and
Nielsen
,
S. R. K.
,
2000
, “
Solution Techniques for Pulse Problems in Non-Linear Stochastic Dynamics
,”
Probab. Eng. Mech.
,
15
(
1
), pp.
25
36
.
34.
Barone
,
G.
,
Navarra
,
G.
, and
Pirrotta
,
A.
,
2008
, “
Probabilistic Response of Linear Structures Equipped With Nonlinear Damper Devices (PIS Method)
,”
Probab. Eng. Phys.
,
23
(
2–3
), pp.
125
133
.
35.
Zhu
,
W. Q.
,
1988
, “
Stochastic Averaging Methods in Random Vibration
,”
Appl. Mech. Rev.
,
41
(
5
), pp.
189
199
.
36.
Zeng
,
Y.
, and
Zhu
,
W. Q.
,
2011
, “
Stochastic Averaging of Strongly Nonlinear Oscillators Under Poisson White Noise Excitation
,”
IUTAM
Symposium on Nonlinear Stochastic Dynamics and Control, Hangzhou, China, May 10–14, pp.
147
155
.
37.
Masud
,
A.
, and
Bergman
,
L. A.
,
2005
, “
Solution of the Four Dimensional Fokker-Planck Equation: Still a Challenge
,”
ICOSSAR
2005, Rome, Italy, June 19–23, pp.
1911
1916
.http://web.engr.illinois.edu/~amasud/Papers/Masud-Bergman-ICOSSAR-OS0506.pdf
38.
Di Matteo
,
A.
,
Di Paola
,
M.
, and
Pirrotta
,
A.
,
2014
, “
Probabilistic Characterization of Nonlinear Systems Under Poisson White Noise Via Complex Fractional Moments
,”
Nonlinear Dyn.
,
77
(
3
), pp.
729
738
.
39.
,
M. A.
, and
Sapsis
,
T. P.
,
2015
, “
Probabilistic Description of Extreme Events in Intermittently Unstable Systems Excited by Correlated Stochastic Processes
,”
SIAM/ASA J. Uncertain. Quantif.
,
3
(
1
), pp.
709
736
.
40.
,
M. A.
,
Cousins
,
W.
, and
Sapsis
,
T. P.
,
2016
, “
A Probabilistic Decomposition-Synthesis Method for the Quantification of Rare Events Due to Internal Instabilities
,”
J. Comput. Phys.
,
322
, pp.
288
308
.
41.
Joo
,
H. K.
,
,
M. A.
, and
Sapsis
,
T. P.
, “
Extreme Events and Their Optimal Mitigation in Nonlinear Structural Systems Excited by Stochastic Loads: Application to Ocean Engineering Systems
,”
Ocean Eng.
,
142
, pp.
145
160
.
42.
Srirangarajan
,
S.
,
Allen
,
M.
,
Preis
,
A.
,
Iqbal
,
M.
,
Lim
,
H. B.
, and
Whittle
,
A. J.
,
2013
, “
Wavelet-Based Burst Event Detection and Localization in Water Distribution Systems
,”
J. Signal Process. Syst.
,
72
(
1
), pp.
1
16
.
43.
Langley
,
R. S.
,
1986
, “
On Various Definitions of the Envelope of a Random Process
,”
J. Sound Vib.
,
105
(
3
), pp.
503
512
.
44.
Ochi
,
M.,K.
,
1990
,
Applied Probability and Stochastic Processes: In Engineering and Physical Sciences
, Vol.
226
,
Wiley-Interscience
, Hoboken, NJ.
45.