A probabilistic decomposition-synthesis method for the quantification of rare events due to internal instabilities
AA probabilistic decomposition-synthesis method for thequantification of rare events due to internal instabilities
Mustafa A. Mohamad, Will Cousins, Themistoklis P. Sapsis ∗ Department of Mechanical Engineering,Massachusetts Institute of Technology,77 Massachusetts Ave., Cambridge, MA 02139October 5, 2018
Abstract
We consider the problem of probabilistic quantification of dynamical systems that haveheavy-tailed characteristics. These heavy-tailed features are associated with rare transientresponses due to the occurrence of internal instabilities. Systems with these properties canbe found in a variety of areas including mechanics, fluids, and waves. Here we develop acomputational method, a probabilistic decomposition-synthesis technique, that takes intoaccount the nature of internal instabilities to inexpensively determine the non-Gaussianprobability density function for any arbitrary quantity of interest. Our approach relieson the decomposition of the statistics into a ‘non-extreme core’, typically Gaussian, anda heavy-tailed component. This decomposition is in full correspondence with a partitionof the phase space into a ‘stable’ region where we have no internal instabilities, and aregion where non-linear instabilities lead to rare transitions with high probability. Wequantify the statistics in the stable region using a Gaussian approximation approach, whilethe non-Gaussian distributions associated with the intermittently unstable regions of thephase space are inexpensively computed through order-reduction methods that take intoaccount the strongly nonlinear character of the dynamics. The probabilistic informationin the two domains is analytically synthesized through a total probability argument. Theproposed approach allows for the accurate quantification of non-Gaussian tails at more than10 standard deviations, at a fraction of the cost associated with the direct Monte-Carlosimulations. We demonstrate the probabilistic decomposition-synthesis method for rareevents for two dynamical systems exhibiting extreme events: a two-degree-of-freedom systemof nonlinearly coupled oscillators, and in a nonlinear envelope equation characterizing thepropagation of unidirectional water waves.
Keywords intermittency; heavy-tails; rare events; stochastic dynamical systems; rogue waves;uncertainty quantification. ∗ Corresponding author: [email protected], Tel: (617) 324-7508, Fax: (617) 253-8689 a r X i v : . [ phy s i c s . c o m p - ph ] N ov Introduction
Quantifying extreme or rare events is a central issue for many technological processes and naturalphenomena. As extreme events, we consider rare transient responses that push the system awayfrom its statistical steady state, which often lead to catastrophic consequences. Complex systemsexhibiting rare events include (i) dynamical systems found in nature, such as the occurrence ofrare climate events [22, 19, 15] and turbulence [30, 32, 44, 63, 31], formation of freak water wavesin the ocean [10, 40, 24, 64, 7, 6]; but also (ii) dynamical systems in engineering applicationsinvolving mechanical components subjected to stochastic loads [38, 39, 41, 52], ship rolling andcapsizing [2, 26, 1], critical events in power grids [25, 42, 61, 27], as well as chemical reactionsand conformation changes in molecules [68, 43, 3, 14].For many systems of practical interest like those above, it has now been well established thatrare transitions occur frequently enough that they are of critical importance. These intermittentevents are randomly triggered while the system evolves on the (stochastic) background attractor,and they are subsequently governed, primarily, by the spatiotemporally local and stronglynonlinear dynamics associated with finite-time instabilities. Systems with these properties pose asignificant challenge for uncertainty quantification schemes [29] and in recent years a wide rangeof research efforts has taken place in various fields towards the quantification and short-termprediction of rare events in complex dynamical systems.The quantification of rare events is one of the most fundamental problems in chemistry.Chemical reactions, conformation changes of molecules, and quantum tunneling are examples ofrare events [68, 43, 3, 14]. These events are rare because the system has to overcome certainbarriers of energetic or entropic nature in order to move from one stable state to the other. Theusual setup for modeling such systems is their formulation in terms of a Langevin equation, i.e.a dynamical system with some non-quadratic potential that has multiple equilibria excited bywhite noise (see for example [14]). The goal then is to study barrier-crossing events by computingtransition rates as well as shortest paths between states. For such systems the classical transitionstate-theory (TST) [18, 68] has been successful in providing the foundation for the developmentof computational tools that determine transition trajectories between different states. However,important limitations for transition state theory may occur when the system potential is notsmooth. In this case it is essential to seek for transition tubes (i.e. ensembles of transitiontrajectories). The statistical framework to analyze such transition-path ensemble is known asthe transition-path theory (TPT) [13, 35, 34] and it has been applied successfully for applied tointeresting and challenging problems in a variety of areas, for example, in chemistry, biology,and material science.Although successful on quantifying transitions between different states, path theory canhave limitations when considering dynamical systems that exhibit rare responses due to theoccurrence of intermittent instabilities (as opposed to multiple equilibria) that lead to strongenergy transfers between modes, as is the case in turbulence or nonlinear waves. In such cases,the rare event is not the result of a transition that takes the system from one state to the other.Rather, it is the result of intermittent instabilities that can ‘push’ the system away from itsstatistical steady state to a dynamical regime with a strongly transient character. The periodthat the system spends away from its statistical steady state attractor as well as the distancefrom the attractor usually takes on a continuous range of values, depending on the intensityand duration of the instability. This situation is completely different from the setup involvingtransitions between discrete states, for which TST and TPT theories have been developed.Large deviations theory [65, 66, 9, 60] is a powerful method for the probabilistic quantificationof extreme events in sequences of probability distributions. It has also been applied in thecontext of stochastic differential equations, known as Freidlin-Wentzell theory [20], as well as forstochastic partial differential equations [5, 58, 59]. In this case, the method essentially provides2s with rates of convergence to probabilistic limits. For example, for a dynamical system excitedby very low intensity noise, large deviations theory gives closed form expressions boundingthe probability for a big deviation of the stochastic solution from the completely unperturbedsolution, i.e. a probabilistic characterization of the stochastic solution relative to the asymptotic(deterministic) limit. Despite its importance, it is not straightforward to apply this frameworkin order to quantify extreme events that rise out of the steady state attractor of the system dueto the occurrence of intermittent instabilities, which is the problem that we are interested here.For the probabilistic description of rare events in phenomena characterized by intermittentinstabilities, the analysis is usually limited to the statistical examination of observed statistics.For example, in ocean engineering, where it is important to analyze the probability of upcrossingsand maxima for various quantities of interests (e.g. wave elevation or mechanical stresses), thestandard setup involves the adoption of globally stable dynamics, for which many techniqueshave been developed (see e.g. [37, 57, 38]). There are numerous technical steps involved in thiscase that lead to elegant and useful results, but the starting point is usually the assumption ofstationarity in the system response, which is not a valid hypothesis for intermittently unstablesystems.Extreme value theory [39, 28, 45, 16, 21] is also a widely applicable method which focuseson thoroughly analyzing the extreme properties of stationary stochastic processes followingvarious distributions. However, even in this case the analysis does not take into account anyinformation about the unstable character of the system dynamics and is usually restricted tovery specific forms of correlation functions for the response statistics [39, 28]. To this end, it isnot surprising that for a large range of complex systems exhibiting intermittent characteristics,the Monte-Carlo method is the only reliable computational approach to arrive at accurateestimates for the tails. However, for high-dimensional systems this direct approach is usuallyprohibitively expensive for practical purposes.Recently, there have been efforts to quantify the heavy tail statistical structures of systemsundergoing transient instabilities. In [36] a probabilistic decomposition was utilized to obtainanalytical approximations for the full probability density function (pdf), including the tailstructure, of systems subjected to parametric (or multiplicative) noise with correlated charac-teristics. In [63] the intermittent behavior of turbulent diffusion models with a mean gradientwas rigorously quantified, while in [44] the capacity of imperfect models to capture intermittentbehavior of turbulent systems was studied.The goal of this work is to present a computational framework for efficiently quantifying thestatistical characteristics of extreme events. We focus on systems where uncertainty interactswith system dynamics to produce intermittent extreme events , i.e. sporadically occurring largeamplitude responses, which give rise to non-Gaussian statistics. Uncertainty could be due tothe initial data, parameters, or the dynamical system itself. The core idea of our probabilisticdecomposition-synthesis method is the separation of intermittent events from the backgroundstochastic attractor, in the spirit of the work done in [36], for low-dimensional systems. Thisdecomposition allows us to apply different uncertainty quantification schemes for the tworegimes (the background and the intermittent component). The background component, althoughpotentially very high-dimensional, can be efficiently described by uncertainty quantificationschemes that resolve low-order statistics. On the other hand, the intermittent component, canbe described in terms of a low-dimensional representation through a small number of localizedmodes. The probabilistic information from these two regimes is synthesized according to atotal probability decomposition argument, in order to approximate the heavy-tailed probabilitydistributions for functionals of interest. Thus, the core of the approach relies on the assumptionthat heavy-tails are primarily due to the action of intermittent instabilities, whereas the ‘mainmass’ of the probability density is due to the the background component.To illustrate the method, we apply the developed framework on two systems that exhibit3ntermittently extreme responses and estimate their statistical distributions. The first example isa nonlinear system of coupled random oscillators, which serves to illustrate the various steps ofthe method in a simple prototype, where the mechanism behind intermittent instabilities is easyto understand although the statistical characteristics of the response are highly non-trivial. Thesecond example, is a nonlinear envelope equation characterizing the propagation of unidirectionalwater waves in deep water. The benefits of the method are highlighted in this complex example,where we are able to capture the statistics at a fraction of the cost of direct numerical simulations.In both cases, we compare our estimates for the statistical distribution with direct Monte-Carloresults, and illustrate the performance of the approach.The paper is structured as follows. In section 2 we describe the problem setup. In section 3we detail the various steps for the proposed decomposition-synthesis method. Section 4 weanalyze and detail the various order reduction schemes for the statistical quantification for rareevents and the background attractor. In section 5 we illustrate the method to the first exampleof a coupled, nonlinear system of random oscillators. In section 6, we demonstrate the methodon the second example to nonlinear waves in deep water.
Let (Ω , B , P ) be a probability space, where Ω is the sample space with ω ∈ Ω denoting anelementary event of the sample space, B the associated σ -algebra of the sample space, and P aprobability measure. (In the following, P X will denote the probability measure associated witha variable X and ρ X the associated probability density function, pdf, if appropriate).Let the dynamical system of interest be governed by the following stochastic Partial Differ-ential Equation (SPDE): ∂u ( x, t ) ∂t = N [ u ( x, t ); ω ] , x ∈ D, t ∈ [0 , T ] , ω ∈ Ω , (1)where N is a general (nonlinear) differential operator with appropriate boundary conditions.We assume the initial state at t = t is random and described by u ( x, t ) = u ( x ; ω ). In whatfollows we utilize a spatial inner product, denoted for two arbitrary functions u ( x ) and v ( x ) as u · v . Extreme events are meant in terms of a norm k · k (e.g. the spatial supreme norm). Inaddition, the ensemble average of a random quantity f ( ω ) is denoted by ¯ f .Here we are interested in determining the statistical distribution for a quantity of interest given by a functional of the solution u ( x, t ) or as a solution of another dynamical systemsubjected to u ( x, t ): q = q [ u ( x, t )] , or dqdt = M [ q, u ( x, t )] . (2)Examples of such quantities could be properties of inertial tracers in turbulent flows, stressesfor ocean structures subjected to water waves, or strains of mechanical components subjected toparametric and/or additive excitations. The computational cost associated with the estimation ofthe heavy-tailed statistics for such quantities is vast given the fact that the dynamical system (1)is characterized by the occurrence of rare events, which define the heavy-tail properties for q .Therefore, application of direct methods, such as Monte-Carlo simulations are prohibitivelyexpensive for these problems. We describe a systematic method to quantify the non-Gaussian response statistics (due torare events caused by intermittent instabilities) in a computationally efficient manner. More4 e Region associated with instabilities that lead to extreme eventsTypical trajectory without extreme eventsRare event creates divergence from background dynamics 𝜁 u u u u r = span{ u } u b = span{ u , u } u b = span{ u , u , u } Figure 1: The conditional decomposition (3) partitions the system response only when a rare eventoccurs due to an instability. This happens when the state of the system enters the instabilityregion R e . In this example the subspace associated with rare events due to instabilities is V s = span { u } .specifically, in this section we first give an overview of the steps involved for the decomposition-synthesis method; more details are provided in the following section. Decomposition and assumptions
We assume that all rare event states due to internal instabilities, defined by the condition k u k > ζ, with ζ being the rare event threshold, ‘live’ in a low dimensional subspace V s . We thendecompose the response of the system as: u ( x, t ) = u b ( x, t ) + u r ( x, t ) , with u r = Π V s [ u ] , if k u k > ζ, and u b = u − u r , (3)where Π V s denotes the linear projection to the subspace V s . A similar decomposition ontoa subspace of interest and a background component, but without taking into account theconditioning on rare events, has been utilized successfully for the uncertainty quantification andfiltering of turbulent systems [50, 33]. This conditional decomposition will allow for the study ofthe two components separately (but taking into account mutual interactions), using different uncertainty quantification methods that (i) take into consideration the possibly high-dimensional(broad spectrum) character of the stochastic background, and (ii) the nonlinear and unstablecharacter of rare events. In this work we are interested in rare events that occur due to transient instabilities . To thisend, we denote as R e the set of all background states that trigger an instability that lead to arare event. In figure 1 we demonstrate the adopted decomposition in a three-dimensional systemwhere the rare event subspace is defined by the linear span of u .The application of this decomposition onto a stochastic background and rare events relies onthe following assumptions:1. The existence of intermittent events have negligible effects on the statistical characteristicsof the stochastic attractor and can be ignored when analyzing the background state u b .2. Rare events are statistically independent from each other.5. Rare events are characterized by low-dimensional dynamics.The first assumption allows for the application of closure models or representation methodsthat can deal with the high dimensional character of the stochastic background attractor. Itexpresses the property that the rare events, although of large magnitude, are localized in timeand space and can induce only negligible modifications to the statistics of the background state.The second assumption follows from the rare character of extreme events. As for the thirdassumption, related to the low-dimensionality of the dynamics of rare events, follows naturallyfrom their spatially or temporally localized character. We emphasize that all these assumptionsdo not imply any restrictions on the dimensionality of the stochastic background state. Analysis of the various regimes
The analysis of the two regimes will consist of the following steps:1.
Order-reduction in the subspace V s in order to model the rare event dynamics, ex-pressed through u r . Then using the approximation u ( x, t ) ’ u r ( x, t ) we will compute theconditional pdf ρ ( q | k u k > ζ, u b ∈ R e ), under the condition that an extreme event occursdue to an internal instability in R e .2. Quantification of the instability region R e using the reduced-order model, by ana-lyzing the conditions that lead to a rare event.3. Description of the background dynamics , expressed through the statistics of u b ,which is not influenced by any internal instabilities in R e . Thus, when the response isdominated only by the background dynamics, we have u ( x, t ) = u b ( x, t ) and the pdf forthe quantity of interest is given by ρ ( q | u = u b ) . Probability for rare events due to internal instabilities P ( k u k > ζ, u b ∈ R e ), whichquantifies the total time/space that the response spends in the rare event regime due tothe occurrence of instabilities. Probabilistic Synthesis
The next step of our technique is to probabilistically synthesize the information obtained fromthe previous analysis. Using a total probability argument, in the spirit of [36], we obtain thestatistics for the quantity of interest q by ρ ( q ) = ρ ( q | k u k > ζ, u b ∈ R e ) P r + ρ ( q | u = u b ) (1 − P r ) , (4)where P r = P ( k u k > ζ, u b ∈ R e ) is the probability of a rare event due to an instability. The firstterm expresses the contribution of rare events due to internal instabilities and is the heavy-tailedpart of the distribution for q . The second term expresses the contribution of the backgroundstate and is the main main probability mass in the pdf for q .This total probability decomposition separates the full response into the conditionallyextreme response and the conditionally background response, weighted by their appropriateprobabilities. The decomposition (4) separates statistical quantities according to the totalprobability law through conditioning on dynamical regimes . In this manner, our approachconnects the statistical quantities that we are interested in with important dynamical regimesthat determine the dominant statistical features (e.g. a Gaussian core due to the backgroundstate and exponential like heavy-tails due to intermittent bursts). An outline of all the stepsinvolved is presented in figure 2. 6 onditional statistics aways from R e Probability of rare events ∥ u ∥ > 𝜁 Conditional statistics during extreme eventsQuantify instability region R e DECOMPOSITION ANALYSIS SYNTHESIS
Total probability lawOrder reduction in V s R e 𝜁 u u u u r = span{ u } u b = span{ u , u } u b = span{ u , u , u } u ( x , t ) = u b ( x , t )+ u r ( x , t ) with u r = Π Vs [ u ] if ∥ u ∥ > 𝜁 and u b = u – u r x P D F x P D F Figure 2: Outline of the steps of the decomposition-synthesis method.For low-dimensional systems, this decomposition can lead to analytical results for theprobability distribution of the response, which would otherwise be impossible to obtain fromthe dynamical system and properties of the noise. For more complex systems, the primarybenefits are computational. We can resolve tail statistics with far fewer simulations than directMonte-Carlo simulation would require, since the decomposition allows for the evaluation of thetails by targeted simulations of rare events , as opposed to the random sampling of Monte-Carlosimulations.We emphasize that the expensive computation in the proposed algorithm is the conditionalstatistics for rare events ρ ( q | k u k > ζ, u b ∈ R e ) . Once this has been obtained, it is simple tocompute rare event statistics for different configurations of the background state, since the other(non-rare) quantities are easy to compute multiple times. This contrasts with the Monte-Carlomethod, where sampling must be performed anew when noise and/or system parameters arechanged. Furthermore, with the proposed method the statistics of arbitrary functionals ofthe response can be easily obtained with minimal additional computational expense, whereasMonte-Carlo simulations may be prohibitively expensive in this case.
We now provide the details for the analysis of each component starting from the reduced-orderdescription of rare events and continuing with the statistical quantification of the backgrounddynamics and the probability of rare events.
The first step for the effective description of the dynamics of rare events is the selection ofa reduced-order basis or a reduced-order subspace denoted as V s . These events are spatiallylocalized structures, to this end we employ localized basis around a neighborhood of an arbitrarypoint x c at time t c : ˆ v n ( x − x c , t − t c ) , n = 1 , . . . , s, (5)where s is the dimension of the subspace. Note that the subspace explicitly depends on thespatiotemporal location ( x c , t c ) of the rare event, which is arbitrary. There are numerous methodswith variable complexity that can be utilized to choose or compute the subspace of rare events.7he simplest choice in this case is a steady localized basis, such as Gabor basis or wavelets [6].However, if there is important spatial translation of the rare event during its lifetime then itmay be beneficial to utilize adaptive methods for the evolution of the basis elements ˆ v n , such asthe dynamically orthogonal field equations [47, 46].Apart from the selection of an appropriate set of basis functions, there is a variety of optionsfor performing the order-reduction of the dynamics. Here we discuss two different approachesfor obtaining the local dynamics, which we later illustrate through two specific problems. The simplest strategy to study the reduced-order dynamics of rare events is to perform aGalerkin projection ignoring the background state u b ( x, t ) for the evolution of the rare eventsinside the reduced-order subspace V s . Such an assumption is valid if the background state plays arole only for the initial triggering of a nonlinear instability , while its small magnitude (comparedwith the intense local responses associated with the nonlinear dynamics in V s ) has very smallinfluence on the reduced-order dynamics of the rare event.Therefore, with the proposed order-reduction strategy we study rare events assuming thattheir evolution is isolated from the stochastic background. The coupling with the backgroundstate is introduced only through the background conditions at t = t c , which is the initial conditionfor the reduced-order dynamics within V s . Based on this setup, we express the intermittentevent in terms of the localized basis that span the subspace V s : u r ( x, t ) = s X n =1 a n ( t )ˆ v n ( x − x c , t − t c ) . (6)Projecting the local (in the sense of x c ) dynamics of the full nonlinear equation on the subspace V s , we obtain the following local dynamical system:˙ a i = N (cid:20) s X n =1 a n ˆ v n ; ω (cid:21) · ˆ v i , i = 1 , . . . , s,a i = ( u b · ˆ v i ) | t = t c , i = 1 , . . . , s, (7)where we again emphasize that the reduced-order system contains information about thebackground state only through the initial condition at t = t c . Through this reduced-orderdescription we can obtain the conditions (in terms of a i and ω ) that define R e , the domain ofattraction to rare events. A more accurate approach for quantifying localized instabilities is the formulation of reduced-order models around the background attractor, instead of a zero background state. This can becritical if the background state not only triggers intermittent instabilities but also determinestheir evolution forward in time.Here we utilize the assumption that the presence of the localized instability u r does nothave important influence on the evolution of the background state u b . We then express the fullsolution as in decomposition (3): u ( x, t ) = Π V ⊥ s [ u ( x, t )] + s X n =1 a n ( t )ˆ v n ( x − x c , t − t c ) , (8)8here the projection operator onto the orthogonal complement is defined as:Π V ⊥ s [ u ] ≡ u − s X n =1 ( u · ˆ v n )ˆ v n . (9)Using this setup we project the original equation on the basis elements ˆ v i to obtain:˙ a i = N (cid:20) Π V ⊥ s [ u ] + s X n =1 a n ˆ v n ; ω (cid:21) · ˆ v i − Π V ⊥ s (cid:20) ∂u∂t (cid:21) · ˆ v i , i = 1 , . . . , s, (10)where the last term on the right hand side vanishes identically. Therefore, the localized dynamicalsystem takes the form:˙ a i = N (cid:20) Π V ⊥ s [ u ] + s X n =1 a n ˆ v n ; ω (cid:21) · ˆ v i , i = 1 , . . . , s,a i = ( u b · ˆ v i ) | t = t c , i = 1 , . . . , s. (11)Clearly, if we set the background state to zero, the dynamical system above reduces to theformulation derived previously. The description of the background state is obtained by projectingthe full equation in V ⊥ s and taking into account that the evolution of u r has negligible effecton the projected dynamics of the background component. However, the infinite dimensionalcharacter of the dynamics in the orthogonal complement, V ⊥ s , makes it impractical to utilize thefull equations for Π V ⊥ s [ u ]. To this end, the reduced-order dynamical system (11) should be seenas a starting point, where appropriate finite-dimensional truncations of Π V ⊥ s [ u ] can be utilizedthat capture the essential effects of the stochastic background on the evolution of the localizedrare events. Such an approach will be demonstrated in the problem involving water waves. To simplify the analysis we consider the case where the reduced-order dynamics have noimportant dependence on ω . We define R e as the set of background states for which the reduced-order dynamics in V s posses at least one finite-time Lyapunov exponent that is positive. Morespecifically, let the flow map φ tt c : V s → V s that maps any point a = Π V s [ u ] to its position attime time t under the effect of the dynamical flow (7): φ tt c : V s → V s , a a ( t, t c , a ; u b ) . (12)Note that this dynamical flow may also depend on the background state u b ∈ V ⊥ s . For each initialcondition a we define the maximum finite-time Lyapounov exponent over some finite-interval τ : λ τ = max i =1 ,...,s log l i (cid:16)(cid:2) ∇ a φ t c + τt c (cid:3) ∗ ∇ a φ t c + τt c (cid:17) , (13)where l i denotes the eigenvalue of the Cauchy-Green tensor, which is by definition symmetricand positive-definite. The domain of attraction to rare events R e is then defined as the set ofbackground states for which we have expansion in the reduced-order subspace: R e = (cid:8) a ∈ V s , u b ∈ V ⊥ s (cid:12)(cid:12) λ τ ( a ; u b ) > (cid:9) . (14) Here we discuss methods for the representation of the statistics of the background stochasticattractor. This refers to the part of the decomposition (3) u b for which instabilities have no role.There are numerous ways to approach the problem and here we review some of the techniquesthat can be used. 9 .3.1 Gaussian closure Based on our setup, the stochastic background u b does not contain intermittent events due toinstabilities; as a consequence, it is reasonable to assume that its statistics can be approximatedby a Gaussian distribution. This assumption can be the starting point for the application ofclosure schemes.For systems that are characterized by a stable mean state u , the finite variance of the steadystate attractor is caused by the external stochastic excitation (see [49] for more details). In thiscase partial linearization of the dynamics or a Gaussian closure of the infinite system of momentequations (see e.g. [17]) is an effective option to capture the conditional statistics of the systemin the state R ce , where we have no rare events occurring. For a such situation the governingequations can be linearized to give ∂u b ∂t = N [ u ; ω ] + L u u b , ω ∈ Ω , (15)where the first term on the right-hand side contains deterministic and stochastic noise termsthat are independent of the state u b , while the second term denotes the linearization of thesystem, L u = δ N δu (cid:12)(cid:12)(cid:12)(cid:12) u (16)around the mean u, where all the eigenvalues have a negative real part (stable mean state). Forsuch system it is straightforward to formulate the second order moment equations and obtainan expression for the Gaussian measure that characterizes the statistics of the backgroundattractor.Such an approach will not be effective if the system under consideration has persistentinstabilities that lead to nonlinear energy transfers between modes, i.e. the mean state is notstable. For such cases other methods may be used to obtain representations of the backgroundattractor statistics, such mean stochastic models [51] or quasi-linear Gaussian closures [48, 51]. A more direct approach involves the numerical simulation of the system and the conditionalsampling of the second-order statistics of the background state using the representation (3). Itis also often the case that the spectrum of the full random field has been measured or estimated.Such an approach is typical, for example, in water waves or other geophysical systems.
For systems having background states that can be modeled by dynamical systems of a specialform, analytical descriptions of their stationary probability measure may be available. Inparticular, for systems excited by white noise one can formulate the corresponding Fokker-Planck-Kolmogorov (FPK) equation defining the evolution of the state pdf [54, 57]. For thespecial case of vibrational systems possessing a Hamiltonian structure perturbed by additiveand/or parametric white noise under linear damping, for special conditions, the stationarymeasure has an explicit form in terms of the Hamiltonian of the system [55, 56]. Furthermore,in [4, 70, 67] an analytical approach is utilized for determining the stationary pdf of more genericvibrational systems, where the steady FPK equation is solved by splitting to simpler partialdifferential equations.To demonstrate the analytical approach, we consider the case of a background state describedby a collection of decoupled vibrational modes, each governed by an equation of the following10orm: ¨ u + δ ˙ u + ∇ V ( u ) = σ ˙ W ( t ; ω ) . (17)Then, the pdf in the statistical steady state will be given in terms of the Hamiltonian H ( u, ˙ u ) ofthe system as [54]: ρ u ˙ u ( u, ˙ u ) = C exp (cid:18) − δσ H ( u, ˙ u ) (cid:19) , (18)where H ( u, ˙ u ) = ˙ u + V ( u ), and C is a normalization constant. We utilize this approach in thefirst example involving a nonlinear vibrational system. P ( k u k > ζ, u b ∈ R e ) Having determined the conditional statistics of the background dynamics, the final step is todetermine the probability of rare transitions, i.e. the likelihood of k u k > ζ when an instabilityoccurs, i.e. u b ∈ R e . This will be defined for an arbitrary point in space, x , through the integral: P ( k u k > ζ, u b ∈ R e ) = 1 T Z t ∈ T ( k u ( x , t ) k > ζ, u b ∈ R e ) dt, (19)where is the indicator function. This integral measures the duration of rare events comparedwith the overall time interval.Note that the above probability is not directly equal to the probability of the backgroundstate crossing into the instability region R e , which is a condition just for the occurrence of arare event and does not contain information regarding the duration of the rare event, which canlast even if the system background has moved outside the instability region.The probability that we are interested in will be found using information obtained from thereduced-order model developed previously. In particular, using the reduced-order model we willobtain the temporal extent T e ( u ) for each rare event that corresponds to any unstable point u ∈ R e . Note that for states u not associated with instabilities T e ( u ) = 0.Then the probability of rare transitions due to instabilities will be approximated by P ( k u k > ζ, u b ∈ R e ) = 1 T Z u T e ( u ) ρ u b ( u ) du, (20)where ρ u b = ρ ( u | u = u b ) is the conditional pdf of the background state, which contains noinformation on rare events. We first illustrate the decomposition-synthesis method on a two-degree-of-freedom systemcomposed of a linear oscillator coupled quadratically to a nonlinear oscillator with cubic stiffness,both under the action of white-noise excitation. The system is given by d xdt + c x dxdt + (cid:0) k x + a x y ( t ) (cid:1) x + sx = σ f x ˙ W x ,d ydt + c y dydt + k y y = σ f y ˙ W y , (21)where x, y are the two state variables, c x , c y > k x , k y > s = 0 is the nonlinearity parameter. The constants11 f x , σ f y > W x , W y are two independent scalar Wiener processes.With a x = 0, the system (21) possesses regimes with intermittent instabilities that lead tofat-tailed equilibrium pdfs in the variable x , whereas the variable y converges to a Gaussian pdf(see figure 3). This system is introduced such that y can trigger an intermittent response in x through a large deviation from its mean value.System (21) is prototypical of the action of intermittency in more complex systems wheresimilar interactions are at play between system modes, e.g. internal instabilities associated withtransfers of energy between modes in turbulence, buckling of beams and plates under stochasticexcitation, ship rolling under parametric stochastic resonance, just to mention a few. Althoughsystem (21) is low-dimensional, quantifying the stationary pdf for the variable x is challenging.For one, the finite time-correlated parametric excitation term due to y on the variable x (which isthe mechanism behind transient instabilities), makes application of the Fokker-Planck equation,unpractical and computationally prohibitive (the resulting equations become high-dimensional).Moreover, there is a nonlinear restoring force ( s = 0) that has to be taken into consideration;this nonlinear term has a significant impact on the pdf of x .Besides its prototypical character, part of the motivation behind the current example is toapply our proposed method to a system where the mechanism behind rare events is transparent.As a result, we do not need to apply the order-reduction schemes described in section 4 forthe rare event component (this will be done in the next application). This will allow for thedemonstration of the decomposition-synthesis procedure in a manner that can be adapted todifferent problems that share similar characteristics. x ( t ) -5 -2 time -20-1001020 y ( t ) -5 -2 Figure 3: Realization of the intermittent variable x ( t ) (top, left) with large amplitude burstsand a fat-tailed equilibrium pdf (top, right), alongside the exciting signal y ( t ) (bottom, left)with a Gaussian response pdf (bottom, right). R e In this case the decomposition step is trivial since all the rare events involve the variable x ,thus V s = span { (1 , T } . In particular, intermittency in the variable x is a consequence of κ ( t ) := k x + a x y ( t ) switching signs from positive to negative values, during which x ( t ) transitionsfrom its regular response to a domain where the likelihood of an instability is high or guaranteed.This switching in κ ( t ) is the triggering mechanism behind localized instabilities in the variable12igure 4: Joint density ρ αξ ( α, ξ ) (left), alongside the marginalized densities for amplitude α (top right) and length scale ξ (bottom right) for c y = 0 . , k y = 1, σ f y = 2 . x ( t ). Therefore, we define the instability region R e = { y | κ = k x + a x y < } . (22) The intensity of the rare events depends both of the average magnitude of κ when this becomesnegative, but also on the duration of the downcrossing. Therefore, we will choose the rare eventdescription presented in section 4.1.2, where the background state is taken into account forthe evolution of the rare events. We choose to parameterize the instability region R e with theaverage amplitude α and average duration ξ of each downcrossing event of κ , which expressesthe background state information.We compute a few realization of the background state d ydt + c y dydt + k y y = σ f y ˙ W y , (23)and for each realization we identify all the regions where κ ( t ) = k x + a x y ( t ) <
0. For each ofthese events, that is, a zero downcrossing followed by a subsequent zero upcrossing, we takethe duration of time that κ < ξ and assign a characteristic amplitude α by taking themean of κ over the duration that κ < α, ξ ). Since we are interested in the probability at a given temporallocation we have an event κ < ξ and amplitude α , we scale the amplitudesin the resulting histogram by their corresponding durations ξ in order to correctly weight thesamples. This gives ρ αξ ( α, ξ | R e ): the probability of finding a rare event characterized by anexcitation duration ξ and amplitude α . We display this pdf in figure 4.Next, we proceed with the computation of the conditional statistics for rare events, ρ x ( x | α, ξ ).This will be based on a few simulation of the equation determining x taking into account the13ackground state (expressed through ( α, ξ )). More specifically, for a given α i , ξ i we evolve therare event dynamics according to d xdt + c x dxdt + ˜ κ ( t ; α i , ξ i ) x + sx = σ f x ˙ W x , (24)where ˜ κ ( t ; α i , ξ i ) is a local representation of the function κ ( t ) = k x + a x y ( t ) in the critical regimewith the same mean amplitude α i (over 0 ≤ t ≤ ξ i ) and duration ξ i . For simplicity, we choosethe following function to approximate κ :˜ κ ( t, α i , ξ i ) = πα i πt/ξ i ) , for 0 ≤ t ≤ ξ i ,k x , otherwise . (25)Initial conditions for x are chosen from the background state, which we will describe later.In this regime, we are only interested in extreme responses for which | x ( t ) | > ζ = σ b , where σ b = ( x | κ> ) is the variance of the background stochastic state. To perform the computationwe simulate an instability according to (24) for a length of time that is long enough to ensurethat x ( t ) returns back to the background regime, where | x ( t ) | < ζ , and sample only the pointsbefore the response relaxes back, i.e. the points for which | x ( t ) | > ζ . Following this procedurewe obtain the quantity ρ x ( x | α, ξ ). Next we determine the statistical characteristics of the background attractor ρ x ( x | R ce ). Inthis regime, by definition, we have no rare events, and we choose to derive the statisticaldistribution analytically. This quantity could also be obtained more directly through theconditional Monte-Carlo approach described in section 4.3.2 by constraining the sampling tothe regime κ = k x + a x y ( t ) > κ > R ce , we choose to approximate κ ( t ) by its conditionally average value κ | κ> , d xdt + c x dxdt + κ | κ> x + sx = σ f x ˙ W x . (26)The steady state pdf can be found by solution of the associated FPK equation: ρ x ˙ x ( q, ˙ q | R ce ) = C exp (cid:18) − c x σ f x H ( q, ˙ q ) (cid:19) , where H ( q, ˙ q ) = 12 ˙ q + 12 κ | κ> q + 14 sq , (27)and C is a normalization constant. Marginalizing out ˙ x in the equation above, gives theconditional pdf for x ( t ) in R ce : ρ x ( q | R ce ) = C exp (cid:20) − c x σ f x (cid:18) κ | κ> q + 14 sq (cid:19)(cid:21) , (28)where C is the normalization constant.To determine κ | κ> we utilize the steady state pdf of y ( t ) which is Gaussian distributed: ρ y ( q ) = s c y k y πσ f y exp (cid:18) − c y k y σ f y q (cid:19) . (29)with variance σ y = σ f y / (2 c y k y ). Therefore, we obtain κ | κ> = k x + a x σ φ (cid:0) − k x a x σ y (cid:1) − Φ (cid:0) − k x a x σ y (cid:1) , (30)14here φ ( · ) is the normal probability distribution function and Φ( · ) the normal cumulativedistribution function. We now proceed to compute the probability of rare events due to instabilities in R e . This canbe found by analyzing the duration of the rare transitions using (24). Denoting by T e ( α i , ξ i ) theduration of the rare event corresponding to an instability in R e of average magnitude α i andduration ξ i , according to (20), we have that P ( x > ζ, y ∈ R e ) = 1 T Z T e ( α, ξ ) ρ αξ ( α, ξ | κ < P ( κ < dαdξ, (31)where P ( κ <
0) = P ( y < − k x /a x ) = Φ (cid:18) − k x a x σ (cid:19) . (32)The probability (31) can be approximated using the analytical argument presented [36], wherethe typical temporal duration of the growth and relaxation phase of a rare event are examined.More specifically, consider a single representative extreme response with an average growth Λ + and decay rate Λ − . During the growth phase we have that u p = u e Λ + T κ< , (33)where T κ< is the temporal duration of the growth event and u p is the peak value of the responseand u an arbitrary initial condition. Similarly, over the decay phase we have u = u p e − Λ − T decay , (34)and by connecting these equations (33) and (34), T decay T κ< = Λ + Λ − . (35)The average duration of an extreme transition is then given by T e = (1+Λ + / Λ − ) T κ< . Therefore,by dividing over the total time duration T we have P ( x > ζ, y ∈ R e ) ’ (cid:18) + Λ − (cid:19) P ( κ < . (36)To leading order (neglecting the nonlinear restoring term) Λ + = p − κ | κ< , since dissipation andadditive forcing have a small role during the growth phase of an extreme response. Therefore,the average growth exponent is Λ + = E (cid:0)p − κ | κ< (cid:1) , which is straightforward to compute using ρ αξ . Similarly, we find that during the decay phase is Λ − = c/
2, and thus Λ − = c/ We have now determined all the components required to construct the distribution of theresponse according to the decomposition-synthesis method. We synthesize the results of theprevious sections by the total probability law: ρ x ( q ) = ρ x ( q | R ce ) (1 − P r ) + P r Z Z ρ x ( q | α, ξ ) ρ αξ ( α, ξ ) dαdξ, (37)15here P r = P ( x > ζ, y ∈ R e ).In figure 5 we compare the decomposition-synthesis results for the equilibrium pdf of x ( t )alongside the ‘true’ density from Monte Carlo simulations, for increasing values of the cubicstiffness parameter for a hardening spring s >
0. Monte-Carlo results are computed using 5000realizations of (21), integrated using the Euler–Maruyama method with time step ∆ t = 0 .
002 to T = 500 time units, discarding the first t = 60 data to ensure a statistical steady state. Overall,we have very good quantitative agreement for both the tails and the core of the distributionbetween the decomposition-synthesis approach and the Monte Carlo results. For larger valuesof the nonlinearity parameter s , tail values are suppressed due to larger restoring forces. Wesee that this effect is also accurately captured for increasing s values in the decomposition-synthesis approach. Here we emphasize the non-uniform decay character of the tails whichcan still be captured very accurately. In particular we have favorable comparison for all threequalitatively different regimes: the core of the distribution, the exponential like heavy-tails atextreme values, and the subsequent sub-exponential decay at very extreme values of distribution(where nonlinearity is very important).We point out that computing the response pdf via the decomposition-synthesis for differentbackground parameters c x , s, σ x parameters is extremely cheap since ρ αξ ( α, ξ ) , which involvesthe rare event dynamics, remains fixed. Also, if we are interested in the response pdf for different k x , a x , this also has minimal computational cost because we can store rare event realizations of y and use them to determine ρ αξ ( α, ξ ), as required. Moreover, the computation to determine ρ αξ ( α, ξ ) is easy and fast, since we do not need a large number of realizations to give goodresults when used in (37) and only requires the simulation of the background variable y ( t ). The second application involves the problem of local extremes for nonlinear dispersive waterwaves. In particular, our goal is to quantify the non-Gaussian, heavy-tailed distribution ofthe local wave field maxima in a nonlinear envelope equation characterizing the propagationof unidirectional water waves. This example has many ocean engineering applications, sincequantifying extreme water waves is critical for ocean structures and naval operations due theircatastrophic consequences. Indeed, extreme waves, termed freak or rogue waves, have causedconsiderable damage to ships, oil rigs, and human life [10, 40, 24, 64]. A wave is termed a freakif its crest-to-trough height exceeds twice the significant wave height H s , with H s being equalto four times the standard deviation of the surface elevation (figure 6). Therefore, such wavesare extreme responses that ‘live’ in the tails of the distribution of the wave elevation.We consider waves on the surface of a fluid of infinite depth and work with approximateequations that govern the dynamics of the wave envelope. In particular, the evolution of aunidirectional, narrow-banded wave field is described well by the modified Nonlinear Schrodinger(MNLS) equation of Dysthe [11], a high order approximation of the fully nonlinear model: ∂u∂t + 12 ∂u∂x + i ∂ u∂x − ∂ u∂x + i | u | u + 32 | u | ∂u∂x + 14 u ∂u ∗ ∂x + iu ∂φ∂x (cid:12)(cid:12)(cid:12)(cid:12) z =0 = 0 , (38)where x is space, t is time, and u is the envelope of the modulated carrier wave. The velocitypotential φ at the surface may be expressed explicitly in terms of u , giving ∂φ/∂x | z =0 = −F − (cid:2) | k |F (cid:2) | u | (cid:3)(cid:3) /
2, where F is the Fourier transform. Equation (38) has been nondimension-alized with x = k ˜ x , t = ω ˜ t , u = k ˜ u , where ˜ x , ˜ t , ˜ u are physical space, time, and envelopevariables, with k the dominant spatial frequency of the surface and ω = √ gk . To lead-ing order, the nondimensionalized surface elevation around the undistributed level is given16 = 0.40 -6 -4 -2 0 2 4 6 x -4 -2 P D F Monte-CarloDecomposition s = 0.60 x -4 -2 P D F Monte-CarloDecomposition s = 0.80 -4 -2 0 2 4 x -4 -2 P D F Monte-CarloDecomposition s = 1.00 -4 -2 0 2 4 x -4 -2 P D F Monte-CarloDecomposition
Figure 5: Probability distribution function approximation with the decomposition-synthesismethod for the variable x ( t ) for the nonlinear system of coupled oscillators (21). The blueline denotes the Monte Carlo simulation and the red line denotes the approximation by thedecomposition-synthesis procedure. The vertical dashed lines denotes 4 standard deviations of x ( t ). Parameters are: c x = 1 , c y = 0 . , k x = 4 , k y = 1, σ x = 0 . , σ y = 2 . η ( x, t ) = Re[ u ( x, t ) e i ( x − t ) ]. The local maxima of the surface elevation is described by themodulus of the envelope | u | , which is the quantity of interest in this example.We consider a wave field with with initial characteristics given by u ( x,
0) = N/ X k = − N/ p k F ( k ∆ k ) e i ( ω k x + ξ k ) , F ( k ) = (cid:15) σ √ π e − k σ (39)where ξ k are independent, uniformly distributed random phases between 0 and 2 π .If the Benjamin-Feir Index (BFI), the ratio of steepness to bandwidth (cid:15)/σ is large enough,we have important probability for the occurrence of rare events due to nonlinear focusing. Suchnonlinear focusing is triggered by the energy localization over a specific region, which is theresult of phase differences between the various harmonics [7] in the stochastic background u b .These relative phases continuously change primarily due to the effect of linear dispersion and asthe BFI increases there is a higher probability for them to result in important energy localizationfor u b and a subsequent nonlinear focusing event u r .17t has been shown [8] that for unidirectional waves the occurrence of nonlinear focusing ofan arbitrary wavegroup (formed due to linear dispersion) in u b is controlled by the wavegroupamplitude A and its lengthscale L . This fact leads us to the adoption of the wavegroupcharacteristics A and L as a way to parameterize the instability region, i.e. the region wherenonlinear focusing occurs. R e We describe the dynamics of focusing wavegroups through a reduced order subspace, V s , obtainedusing a proper orthogonal decomposition [23] of the focusing wavegroups under MNLS dynamics.The proper orthogonal decomposition is appropriate in this case as the order reduction schemefor focusing wavegroups, since we do not have important spatial translations of the focusingwaves and this allows us to capture the dynamics with just a few modes.To capture the variations in the dynamics for different wavegroup lengthscales and amplitudes,we choose n = 8 sets of simulations of the MNLS for various ( A, L ) that undergo nonlinearfocusing. For each simulation we take snapshots u s ( x, t k ; A i , L i ) in time, where for each simulationwe ensure that the snapshots are capturing the dominant focusing action. Stacking the snapshotsfor all the simulations gives the matrix XX = (cid:2) u s ( x, t ; A , L ) u s ( x, t ; A , L ) . . . u s ( x, t m ; A n , L n ) (cid:3) , (40)where m is total number of snapshots. We use the method of snapshots [53] to determine theorthonormal POD modes, by solving the m × m eigenvalue problem X T XU = U Λ, for themodes V = XU Λ − / .By this procedure we obtain the local basis ˆ v n ( x, t ) and represent u r as u r ( x, t ) = s X n =1 a n ( t )ˆ v n ( x − x c , t − t c ) , (41)which is conditional upon a background state at ( x c , t c ). We used the projection of the fullMNLS around a zero background state as described in section 4.1.1. We found that a projectionupon s = 14 modes approximated the dynamics of V s well across a range of L and A (fig. 7).With this reduced order description of the dynamics for the intermittent component u r , weinvestigate the evolution of various wavegroups ( A, L ) in order to quantify the conditions forthe occurrence of rare transitions, which will determine R e . More specifically, we consider initialwavegroups of the form u r ( x, t c ) = Π V s [ A sech( x/L )] . We compute the value of the first spatiotemporal local maximum of u r ( x, t ) for a range of( A, L ) to investigate group dynamics under different lengthscales L and amplitudes A . In the space S u r f a c e E l e v a t i o n Figure 6: Example large amplitude freak wave formed due to nonlinear interactions.18igure 7: Simulation of an extreme wave group for A = 0 . , L = 11, comparing the exactMNLS (left) and the reduced order model (right) with 14 modes.left pane of figure 8 we display the value of this map divided by the initial amplitude A . Thiswave group amplification factor is 1 for defocusing groups and greater than 1 for groups thatundergo nonlinear focusing. Importantly, by this map we can partition the space (right paneof figure 8) and subscribe the region where u max ( A, L ) /A > λ T > u r ( x, t ) (conditional on the backgroundstate ( A, L )) is a focusing wave packet. This gives the set R e = (cid:26) ( A, L ) (cid:12)(cid:12)(cid:12)(cid:12) A u max ( A, L ) > (cid:27) , (42)which parameterizes the instability region in terms of just two parameters ( A, L ).Figure 8:
Left:
The wave group amplification factor 1
A u max ( A, L ) computed using the reduced-order dynamical system for the rare events u r ( x, t ). Right:
Partition of (
A, L ) plane into focusingand defocusing regions according to the reduced order model.19igure 9: Procedure for computing the conditional statistics for rare events. First we determinethe probability of wavegroups, ρ AL ( A, L ), in R e for a given spectrum. Then we compute thepdf for the quantity of interest for each focusing group, using the reduced order model. To determine ρ | u | ( q | k u k > ζ, u b ∈ R e ) we express it in terms of the wavegroup parameters, i.e.in the form ρ | u | ( q | k u k > ζ, A, L ∈ R e ), where A and L denote the parameters of an arbitrarywavegroup formed in the stochastic background u b and lead to a rare event defined by ζ = σ ,the standard deviation of the initial spectrum. In the rare event regime we approximate u ’ u r .The estimation of this quantity involves simulations of u r through the reduced-order model. Weperform the simulation until the wavegroup relaxes back to the background state, necessarilya short time simulation. We also sample spatial points for which the response has importantmagnitude, in practice this means that we consider points x ∈ [ − L, L ]. This is done forall unstable wavegroup parameters ( A, L ) for which the probability of occurrence ρ AL ( A, L )within the background dynamics is finite. A visual demonstration of this procedure is displayedin figure 9.We explicitly compute ρ AL ( A, L ) by generating many random fields according to (39). Thisrandom field represents typical realizations of the background state u b , where the dominanteffect is linear dispersion, which continuously mixes the phases between harmonics. We notethat we ignore the evolution of the spectrum due to weak nonlinearities (typically the spectrumtends to broaden [12, 69]) and we use the initial spectrum as the spectrum of the background u b . This simplification, which could be avoided by sampling the spectrum of the steady stateusing a relatively short simulation, does not cause any serious discrepancies to our final results.However, for more complex cases, such as two-dimensional nonlinear waves, where the spectrumis continuously evolving, the transient character of the background statistics may have to betaken into consideration. Our framework can support such situation without any modifications.20igure 10: Group density ρ AL ( A, L ) for Gaussian spectra with (cid:15) = 0 . σ = 0 . (cid:15) = 0 . σ = 0 . R e boundary.For each random field realization, we apply a group detection procedure (described in detailin [7]), which returns a set of the groups in the field along with the amplitude and length scaleof each group. This is a fast computation since it only requires generating and analyzing randomrealizations out of a given spectrum. After generating many random fields and computing thegroups, we have a set of samples of ( A, L ); subsequently we estimate the joint density ρ AL ( A, L ),scaling the amplitudes in the resulting pdf by their corresponding lengthscale.Example group densities for two different spectra are displayed in figure 10. In each figure weoverlay the R e boundary. Notice that in the case displayed on the right panel, the spectral widthis larger, meaning that groups are narrower. The effect of this change is to reduce the numberof focusing groups, and thus reduce the number of extreme waves. This is consistent with theBenjamin-Feir Index (2 √ (cid:15)/σ ), used in the nonlinear water-waves community to evaluate thelikelihood for rare events, which is reduced by increasing σ . The next ingredient is to estimate the conditional statistics for the background dynamics ρ | u | ( q | u = u b ). We will have only background components if and only if the occurringwavegroups belong in R ce . As a result, we have that ρ | u | ( q | u = u b ) = ρ | u | ( q | A, L ∈ R ce ) . Moreover, in this regime we have predominantly linear dynamics and therefore the statistics of thebackground state can be approximated by a Gaussian distribution. More accurate representationsthat take into account weak nonlinearities [62] may also be utilized to improve the accuracy ofthe main core of the distribution.For linear waves η ( t ) with a Gaussian distribution their envelope is Rayleigh distributed.Therefore, the conditional distribution for the quantity of interest in the background state isgiven by ρ | u | ( q | R ce ) = qσ s e − q / σ s , (43)where σ s is the variance of non-focusing wave groups, computed by taking the variance of21ave groups in R ce . This is estimated along side the computation described in section 6.4 fordetermining ρ AL ( A, L ), by using the scale selection algorithm to identify the wave groups in R ce and then by taking the variance of this set of wave packets. The final component is the computation of the total probability to have a rare event in anarbitrary spatial location. We assume that the system is ergodic, and we use (20) in the followingspatial form P ( k u k > ζ, u b ∈ R e ) = 1 X D Z u X e ( u ) ρ u b ( u ) du, (44)where X e ( u ) is the spatial extend of each wavegroup associated with nonlinear focusing, and X D is the spatial size of the domain. Taking into account that the pdf ρ AL already has beenweighted with respect to the spatial extend of the wavegroups the last formula takes the form P r = P ( k u k > ζ, u b ∈ R e ) = Z Z R e ρ AL ( A, L | R e ) dAdL. (45)This probability, of course, depends on the particular choice of the spectral properties of thefield. Increasing the spectral width, for example, will shift the distribution of ρ AL to the left,reducing in this way the total likelihood for the occurrence of a rare event (figure 10). We have now determined all the components required to compute the heavy tailed distributionfor | u | according to the decomposition-synthesis procedure. We have the following final result,by the total probability probability argument: ρ | u | ( q ) = ρ | u | ( q | R ce )(1 − P r ) + P r Z Z R e ρ | u | ( q | A, L ) ρ AL ( A, L ) dAdL. (46)In figure 11 we demonstrate the decomposition-synthesis method on four different spectraalongside the ‘true’ density from Monte-Carlo simulations. To compute the Monte-Carlo statisticswe use 10 realizations of the MNLS equation on a spatial domain 256 π and a total sampledduration of 1000 time units, discarding the first 500 time units for cases with σ = 0 .
10 case and500 time units for the cases with σ = 0 .
20 case. Overall, we see that the decomposition-synthesisapproach gives good approximations to the true density. For the most heavy-tailed cases, weobserve that large waves occur many orders of magnitude more frequently than predicted bylinear dynamics (Rayleigh distribution). Our method compares favorably with the Monte-Carloresults at a fraction of the computational cost.We note that the computational savings of the decomposition-synthesis method is due toseveral important advantages over the direct brute force Monte-Carlo method. First and veryimportantly, we do not have to wait for extreme waves to emerge from the background wavefield; moreover, we do not have to wait until a sufficient number of extreme waves occur in orderto obtain reliable tail statistics. This is because we simulate extreme wave groups directly bycarefully choosing the initial/background conditions that trigger them. The fact that we usea low dimensional reduced order model to evolve extreme wave groups, makes applying ourestimation procedure very inexpensive for computing the extreme components. And finally, thedistribution of the background dynamics was obtained analytically.22e also emphasize that once we performed the decomposition-synthesis procedure, computingthe distribution for different background spectra is extremely cheap, since the only quantitythat needs to be recomputed is the distribution of various groups ρ AL ( A, L ) and σ s , which iseasy and fast to determine since it does not require the simulation of the original SPDE or thesolution of the reduced-order model. = 0.025 σ = 0.200.05 0.1 0.15 0.2 |u| -5 P D F Monte-CarloLinear TheoryDecomposition = 0.025 σ = 0.100.05 0.1 0.15 0.2 |u| -4 -2 P D F Monte-CarloLinear TheoryDecomposition = 0.050 σ = 0.200.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 |u| -6 -4 -2 P D F Monte-CarloLinear TheoryDecomposition = 0.050 σ = 0.100.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 |u| -4 -2 P D F Monte-CarloLinear TheoryDecomposition
Figure 11: Probability distribution function approximation by the decomposition-synthesismethod for the MNLS equation (38). In each case, the blue curve denotes Monte Carlo simulationof the MNLS and the red line represents the approximation by the decomposition-synthesisprocedure. The vertical dashed line denotes 4 (cid:15) (i.e. 4 standard deviations) and x -axis plotted to8 (cid:15) . We also display the Rayleigh distribution of | u | (green dashed line). Ordered in increasingBFI regimes: top left BFI = 0 .
35, top right BFI = 0 .
71, bottom left BFI = 0 .
71, bottom rightBFI = 1 . We have considered the problem of quantifying rare event statistics due to internal and finite-timeinstabilities in general dynamical systems. Our analysis is based on the assumption that theconditional dynamics involving rare events are of low-dimensionality, although the total response23f the system can still be very high dimensional. Relying on this setup, we have formulated aconditional decomposition into low-dimensional extreme events caused by internal instabilities,and a high dimensional stochastic background. This decomposition allows for the study ofthe two components separately (but taking into account mutual interactions), using different uncertainty quantification methods that (i) take into consideration the possibly high-dimensional(broad spectrum) character of the stochastic background, and (ii) the nonlinear and unstablecharacter of the rare events.The adopted decomposition is in full correspondence with a partition of the phase spaceinto a stable region, where we have no internal instabilities, and a region where non-linearinstabilities lead to extreme transitions with high probability. We quantify the statistics in thestable region using a Gaussian approximation, while the non-Gaussian distributions associatedwith the intermittently unstable regions of the phase space, are inexpensively computed throughorder-reduction methods that take into account the strongly nonlinear character of the dynamics.The probabilistic information for the two domains is analytically synthesized through a totalprobability argument.The proposed approach allows for the derivation of the statistics for any quantity of interestin a semi-analytical form, where only a few carefully selected simulations through a reducedorder model are sufficient for the accurate determination of the heavy tail structure. For low-dimensional systems the developed framework allows for the derivation of fully analytical forms,while for more complex systems it provides an inexpensive computational method to determineextreme event statistics.To demonstrate the new method we considered two systems of increasing complexity wherenon-trivial energy exchanges occur due to internal instabilities, leading to extreme responses.The first application was a two-degree-of-freedom nonlinear system of coupled mechanicaloscillators encountered in a viariety of engineering settings. This setup leads to non-Gaussianstatistics with heavy tails characterized by qualitatively different regimes. Through numericalexperiments we demonstrated that our method is able to capture very accurately the core ofthe response distribution, the exponential like heavy-tails at extreme values, and the subsequentsub-exponential decay at very extreme values of distribution.The second application was a prototype nonlinear envelope equation, that describes theone-dimensional propagation of deep water waves, where extreme waves (known as rogue orfreak waves) randomly appear due to nonlinear focusing in a wavefield. This is an example thatis particularly challenging due to its nonlinear, dispersive, and infinite dimensional character.In such case our method is very advantageous as it allows to separately quantify the extremewavegroups from the background field. Comparisons with direct Monte-Carlo simulation demon-strated the effectiveness of our approach on semi-analytically and inexpensively capturing theheavy-tailed statistics for the distribution of the local wave field maxima, for four differentspectra of increasing heavy-tailed statistics . We also demonstrated the value of our approach oncomputing the pdf of interest for different spectrum parameters in comparison with Monte-Carloresults where the simulation would have to be run anew.
Acknowledgments
TPS has been supported through the Office of Naval Research grant ONR N00014-14-1-0520,the Army Research Office Young Investigator Award 66710-EG-YIP, and the DARPA grantHR0011-14-1-0060. MAM and WC have been partially supported by the first and second grants.24 eferences [1] L. Arnold, I. Chueshov, and G. Ochs. Stability and capsizing of ships in random sea - Asurvey.
Nonlinear Dynamics , 36:135–179, 2004.[2] V. L. Belencky and N. B. Sevastianov.
Stability and Safety of Ships: Risk of Capsizing .The Society of Naval Architects and Marine Engineers, 2007.[3] P. G. Bolhuis. Transition-path sampling of beta-hairpin folding.
Proceedings of the NationalAcademy of Sciences of the United States of America , 100(21):12129–34, oct 2003.[4] T. K. Caughey and F. Ma. The exact steady-state solution of a class of non-linear stochasticsystems.
Int. J. Non. Linear. Mech. , 17:137–142, 1982.[5] P.-L. Chow. Some parabolic Ito equations.
Comm. Pure Appl. Math. , 45:97, 1992.[6] W. Cousins and T. P. Sapsis. Quantification and prediction of extreme events in a one-dimensional nonlinear dispersive wave model.
Physica D , 280:48–58, 2014.[7] W. Cousins and T. P. Sapsis. Reduced order prediction of rare events in unidirectionalnonlinear water waves.
Submitted to Journal of Fluid Mechanics , 2015.[8] W. Cousins and T. P. Sapsis. The unsteady evolution of localized unidirectional deep waterwave groups.
Physical Review E , 91:063204, 2015.[9] A. Dembo and O. Zeitouni.
Large Deviations Techniques and Applications . Springer-Verlag,New York, 2000.[10] K. Dysthe, H. E. Krogstad, and P. Müller. Oceanic Rogue Waves.
Annu. Rev. Fluid Mech. ,40(1):287, 2008.[11] K. B. Dysthe. Note on a modification to the nonlinear Schrodinger equation for applicationto deep water waves.
Proceedings of the Royal Society of London. A. Mathematical andPhysical Sciences , 369(1736):105–114, 1979.[12] K. B. Dysthe, K. Trulsen, H. E. Krogstad, and H. Socquet-Juglard. Evolution of anarrow-band spectrum of random surface gravity waves.
J. Fluid Mech , 478:1–10, 2003.[13] W. E. and E. Vanden-Eijnden. Towards a Theory of Transition Paths.
Journal of StatisticalPhysics , 123(3):503–523, may 2006.[14] W. E and E. Vanden-Eijnden. Transition-path theory and path-finding algorithms for thestudy of rare events.
Annual review of physical chemistry , 61:391–420, jan 2010.[15] D. R. Easterling. Climate Extremes: Observations, Modeling, and Impacts.
Science ,289(5487):2068–2074, sep 2000.[16] P. Embrechts, C. Kluppelberg, and T. Mikosch.
Modeling Extremal Events . Springer, 2012.[17] E. S. Epstein. Stochastic Dynamic Predictions.
Tellus , 21:739–759, 1969.[18] H. Eyring. The Activated Complex and the Absolute Rate of Chemical Reactions.
ChemicalReviews , 17(1):65–77, aug 1935.[19] C. Frei and C. Schär. Detection Probability of Trends in Rare Events: Theory andApplication to Heavy Precipitation in the Alpine Region.
Journal of Climate , 14(7):1568–1584, apr 2001. 2520] M. I. Freidlin and A. D. Wentzell.
Random Perturbations of Dynamical Systems . BerlinSpringer Verlag, 2nd editio edition, 1998.[21] J. Galambos.
The asymptotic theory of extreme order statistics . Wiley Series in Probabilityand Statistics, 1978.[22] A. Hense and P. Friederichs. Wind and Precipitation Extremes in the Earth’s Atmosphere.In
Extreme Events in Nature and Society , pages 169–187. Springer-Verlag, 2006.[23] P. Holmes, J. Lumley, and G. Berkooz.
Turbulence, Coherent Structures, Dynamical Systemsand Symmetry . Cambridge University Press, 1996.[24] C. Kharif and E. Pelinovsky. Physical mechanisms of the rogue wave phenomenon.
EuropeanJournal of Mechanics, B/Fluids , 22(6):603–634, 2003.[25] V. Kishore, M. S. Santhanam, and R. E. Amritkar. Extreme events and event sizefluctuations in biased random walks on networks.
Physical Review E , 85:56120, 2012.[26] E. Kreuzer and W. Sichermann. The effect of sea irregularities on ship rolling.
Computingin Science and Engineering , May/June:26–34, 2006.[27] P. Kundur.
Power System Stability and Control . McGraw-Hill Education, 1994.[28] M. R. Leadbetter, G. Lindgren, and H. Rootzen.
Extremes and related properties of randomsequences and processes . Springer, New York, 1983.[29] A. J. Majda and M. Branicki. Lessons in Uncertainty Quantification for Turbulent DynamicalSystems.
Discrete and Continuous Dynamical Systems , 32:3133–3221, 2012.[30] A. J. Majda and J. Harlim.
Filtering Complex Turbulent Systems . Cambridge UniversityPress, 2012.[31] A. J. Majda and Y. Lee. Conceptual dynamical models for turbulence.
Proceedings of theNational Academy of Sciences of the United States of America , 111(18):6548–53, may 2014.[32] A. J. Majda, D. W. McLaughlin, and E. G. Tabak. A one-dimensional model for dispersivewave turbulence.
Journal of Nonlinear Science , 7(1):9–44, 1997.[33] A. J. Majda, D. Qi, and T. P. Sapsis. Blended particle filters for large-dimensional chaoticdynamical systems.
Proceedings of the National Academy of Sciences , 111(21):7511–7516,2014.[34] P. Metzner, C. Schütte, and E. Vanden-Eijnden. Transition Path Theory for Markov JumpProcesses.
Multiscale Modeling & Simulation , 7(3):1192–1219, jan 2009.[35] Philipp Metzner, Christof Schütte, and Eric Vanden-Eijnden. Illustration of transition paththeory on a collection of simple examples.
The Journal of chemical physics , 125(8):084110,aug 2006.[36] M. A. Mohamad and T. P. Sapsis. Probabilistic description of extreme events in inter-mittently unstable systems excited by correlated stochastic processes.
SIAM ASA J. ofUncertainty Quantification , 3:709–736, 2015.[37] A Naess. Extreme value estimates based on the envelope concept.
Applied Ocean Research ,4(3):181, 1982. 2638] A. Naess and T. Moan.
Stochastic Dynamics of Marine Structures . Cambridge UniversityPress, 2012.[39] M. Nicodemi. Extreme value statistics. in Encyclopedia of complexity and systems science ,E:3317, 2009.[40] M. Olagnon and M. Prevosto.
Rogue Waves 2004: Proceedings of a Workshop Organizedby Ifremer and Held in Brest, France, 20-21-22 October 2004, Within the Brest Sea TechWeek 2004 . Editions Quae, 2005.[41] J. K. Paik and A. K. Thayamballi.
Ultimate Limit State Design of Steel-Plated Structures .John Wiley and Sons, 2003.[42] P. Pourbeik, P. Kundur, and C. Taylor. The anatomy of a power grid blackout - Root causesand dynamics of recent major blackouts,.
IEEE Power and Energy Magazine , 4(5):22–29,2006.[43] L. R. Pratt. A statistical method for identifying transition states in high dimensionalproblems.
The Journal of Chemical Physics , 85(9):5045, nov 1986.[44] D. Qi and A. J. Majda. Predicting Fat-Tailed Intermittent Probability Distributions inPassive Scalar Turbulence with Imperfect Models through Empirical Information Theory.
Submitted to Physica D , 2015.[45] R.-D. Reiss and M. Thomas.
Statistical analysis of extreme values . Birkh{ä}user; 3rdedition, 2007.[46] T. P. Sapsis. Attractor local dimensionality, nonlinear energy transfers, and finite-timeinstabilities in unstable dynamical systems with applications to 2D fluid flows.
Proceedingsof the Royal Society A , 469(2153):20120550, 2013.[47] T. P. Sapsis and P. F. J. Lermusiaux. Dynamically Orthogonal field equations for continuousstochastic dynamical systems.
Physica D , 238:2347–2360, 2009.[48] T. P. Sapsis and A. J. Majda. A statistically accurate modified quasilinear Gaussian closurefor uncertainty quantification in turbulent dynamical systems.
Physica D , 252:34–45, 2013.[49] T. P. Sapsis and A. J. Majda. Blended reduced subspace algorithms for uncertaintyquantification of quadratic systems with a stable mean state.
Physica D , 258:61, 2013.[50] T. P. Sapsis and A. J. Majda. Blending Modified Gaussian Closure and Non-GaussianReduced Subspace methods for Turbulent Dynamical Systems.
Journal of Nonlinear Science ,23:1039, 2013.[51] T. P. Sapsis and A. J. Majda. Statistically Accurate Low Order Models for UncertaintyQuantification in Turbulent.
Proceedings of the National Academy of Sciences , 110:13705–13710, 2013.[52] D. Shadman and B. Mehri. A non-homogeneous Hill’s equation.
Applied Mathematics andComputation , 167:68–75, 2005.[53] L. Sirovich. Turbulence and the dynamics of coherent structures, Parts I, II and III.
Quart.Appl. Math. , XLV:561–590, 1987.[54] K. Sobczyk.
Stochastic Differential Equations . Kluwer Academic Publishers, Dordrecht,The Netherlands, 1991. 2755] C. Soize. Steady-state solution of Fokker-Planck equation in higher dimension.
Probab.Eng. Mech. , 3:196–206, 1988.[56] C. Soize.
The Fokker-Planck equation for stochastic dynamical systems and its explicitsteady state solutions . World {S}cientific, 1994.[57] T. Soong and M. Grigoriu.
Random Vibration of Mechanical and Structural Systems . PTRPrentice Hall, 1993.[58] R. Sowers. Large Deviations for a reaction diffusion equation with non-Gaussian perturba-tions.
Ann. Probab. , 20:504, 1992.[59] S. S. Sritharan and P. Sundar. Large deviations for the two-dimensional Navier–Stokesequations with multiplicative noise.
Stochastic Processes and their Applications , 116:1636,2006.[60] D. Stroock.
An Introduction to the Theory of Large Deviations . Springer-Verlag, New York,1984.[61] Y. Susuki and I. Mezic. Nonlinear Koopman Modes and a Precursor to Power SystemSwing Instabilities.
IEEE Transactions on Power Systems , 27(3):1182–1191, aug 2012.[62] M. A. Tayfun. Narrow-band nonlinear sea waves.
Journal of Geophysical Research ,85(C3):1548, 1980.[63] X. Tong and A. J. Majda. Intermittency in Turbulent Diffusion Models with a MeanGradient.
Submitted to Nonlinearity , 2015.[64] K. Trulsen and K. B. Dysthe. A modified nonlinear Schrödinger equation for broaderbandwidth gravity waves on deep water.
Wave motion , 24(3):281–289, 1996.[65] S. R. S. Varadhan.
Large Deviations and Applications . SIAM, 1984.[66] S. R. S. Varadhan. Special invited paper: Large deviations.
The Annals of Probability ,36:397, 2008.[67] R. Wang and Z. Zhang. Exact stationary solutions of the Fokker-Planck equation fornonlinear oscillators under stochastic parametric and external excitations.
Nonlinearity ,13:907–920, 2000.[68] E. Wigner. The transition state method.
Transactions of the Faraday Society , 34:29, jan1938.[69] W. Xiao, Y. Liu, G. Wu, and D. K. P. Yue. Rogue wave occurrence and dynamics by directsimulations of nonlinear wave-field evolution.
Journal of Fluid Mechanics , 720:357–392,2013.[70] W. Q. Zhu. Exact Solutions for Stationary Responses of Several Classes of NonlinearSystems to Parametric and/or External White Noise Excitations.