## Abstract

A probabilistic approach is needed to address systems with uncertainties arising in natural processes and engineering applications. For computational convenience, however, the stochastic effects are often ignored. Thus, numerical integration routines for stochastic dynamical systems are rudimentary compared to those for the deterministic case. In this work, the authors present a method to carry out stochastic simulations by using methods developed for the deterministic case. Thereby, the well-developed numerical integration routines developed for deterministic systems become available for studies of stochastic systems. The convergence of the developed method is shown and the method's performance is demonstrated through illustrative examples.

## 1 Introduction

Natural processes are inevitably uncertain, and systems in engineering commonly have uncertainties. To capture this stochasticity, one needs a probabilistic approach. Thus, stochastic models are widely used in physics (e.g., Gardiner [1]), engineering (e.g., Wirsching et al. [2]), and biology (e.g., Wilkinson [3]). Uncertain parameters are classically related to manufacturing imperfections, finite measurement resolution, or incomplete data, whereas random loads arise in complex environments due to wind, ocean waves, seismic excitations, or road roughness. However, for computational convenience, stochasticity is often ignored in studies of associated systems. Thus, a variety of methods exist to simulate deterministic systems, while the same is not true for carrying out simulations of stochastic systems. In particular, for high dimensional systems, methods are still in their infancy.

The response of stochastic dynamical systems is fully prescribed by the time evolution of the probability density function (PDF). If the stochasticity is generated by a Wiener process, then the PDF is governed by the Fokker-Planck equation, a partial differential equation (e.g., Risken [4]). Exact solutions have been obtained by, for example, Soize [5], Caughey [6], and Lin and Cai [7]. To obtain these results one requires a balance between dissipative terms and white noise perturbations. This balance, however, is hardly met in realistic systems (Lin and Cai [7]).

In the absence of exact solutions to the Fokker–Planck equation, it can be solved with numerical methods. In principle, space and/or time can be discretized to obtain a finite-dimensional dynamical system. Examples include methods based on finite differences (e.g., Pichler et al. [8]), the path integral method (e.g., Wehner and Wolfer [9] or Yu et al. [10]) and a finite element discretization formulated by Spencer and Bergman [11]. However, Masud and Bergman [12] point out that the computational burden is so significant, that the applications of those methods are limited to two or three-dimensional systems. In fact, the memory requirements of such methods generally grow exponentially with the dimensions of the considered dynamical system.

The most common and effective approach to approximate solutions to stochastic differential equations are discrete-time approximations (cf. Kloeden and Platen [13]). With this approach, the computational costs and memory requirements grow only polynomially with respect to number of dimensions. However, for example, as remarked by Rüemelin [14], numerical time integration schemes derived to approximate solutions of deterministic ordinary differential equations do not carry over to the stochastic setting in a straightforward manner. An underlying reason is that additional terms arise when solutions to stochastic differential equations are expanded in a Taylor series. These terms lead to the so-called Wagner-Platen formula (cf. Wagner and Platen [15] or Milstein [16]). Therefore, numerical integration of stochastic differential equations requires fundamentally different algorithms.

Schurz through his chapter in Ref. [17], Kloeden and Platen [13], and the introductory treatment by Platen [18] offer an extensive list of available integration schemes. The most prominent and basic ones are the Euler–Maruyama scheme and Milstein's method (e.g., Milstein Ref. [16]). Higher order schemes include the stochastic Runge–Kutta in which a predictor-corrector scheme is employed similar to the deterministic equivalent to enhance numerical stability (Schurz in Ref. [17]). Moreover, Boyce [19] developed a scheme with adaptive step size control. Milstein and Tret'yakov [20,21] present integration methods for stochastic differential equations with small noise terms. Such a setting is especially relevant for engineering applications, wherein the stochastic part can often be considered as small or weak. However, the computational expense of such higher-order schemes can be significant (e.g., Schurz in Ref. [17]), so much so, that these costs outweigh their benefits (cf. Mannella [22]).

In terms of the current state of the art, the available numerical methods to solve the Fokker–Plank equation suffer from the curse of dimensionality, whereas the discrete-time approximations of stochastic differential equations are often limited to the stochastic equivalent of the forward Euler scheme. However, in many applications, accurate and effective approximations of solutions to the deterministic part already require intricate numerical integration routines employing, for example, adaptive step-size control (cf. Gear [23]), predictor-corrector schemes (cf. de-Jalon and Bayo [24]), and customized methods such as Newmark's method for structural dynamics (cf. Géradin and Rixen [25]). Unfortunately, these well-developed routines are not applicable in a stochastic setting.

Here, the authors propose an algorithm to extend any deterministic numerical integrator to a stochastic setting. To that end, they proceed as follows. In many applications, especially in engineering, the noise intensity is small and one can perform a parameter expansion (cf. Sec. 3). In this setting, the authors formulate a small noise integrator (SNI), which relies on an appropriate resampling of distributions along sample paths (cf. Sec. 4.1). Within this algorithm, any deterministic integrator can be used to approximate the sample paths of stochastic dynamical systems. After proving the strong convergence of the SNI, the algorithm performance is demonstrated by considering various nonlinear mechanical systems with up to one hundred degrees-of-freedom in Sec. 4.2. Furthermore, in Sec. 5.1, a Gaussian kernel is presented to rigorously deduce a smooth distribution from a finite number of samples and the corresponding performance is shown for nonlinear mechanical systems (cf. Sec. 5.2).

## 2 System Setting

**x**denotes the states and the integer

*N*is the numbers of degrees-of-freedom. The parameter

*σ*in Eq. (1) is small and scales the stochastic excitation defined in differential form as

where $W1,W2,\u2026,WM$ denote uncorrelated one-dimensional standard Wiener processes (e.g., Kloeden and Platen [13]). The *m*-th column of the matrix $B(x,t)$ prescribes the direction in which the *m*-th Wiener process acts. The directions can generally depend on the states as indicated in Eq. (2). Such multiplicative noise commonly arises, for example,, in micro-electromechanical devices (MEMS, cf. Vig and Yonkee [26]). The parameter *σ* can be considered to be representative of small uncertainities that arise in various settings, including high precision manufacturing, measurements, and environmental modeling.

The time dependence of the deterministic part $f(x,t)$ and the stochastic part $B(x,t)$ is general. Thus, system (1) includes autonomous (i.e., time-invariant) systems as well as systems with time-periodic, quasi-periodic, and even more general, possibly aperiodic, time dependence. Moreover, it is assumed that the matrix $B(x,t)$ and the vector $f(x,t)$ are at least twice continuously differentiable with respect to the spatial coordinate **x**.

**q**is the vector of coordinates and the integer $N*$ denotes the numbers of degrees-of-freedom. The mass $M\u2208\mathbb{R}N*\xd7N*$ is assumed to be invertible, whereas the damping $C\u2208\mathbb{R}N*\xd7N*$ and stiffness $K\u2208\mathbb{R}N*\xd7N*$ matrix are assumed to be general matrices. In many engineering applications the aforementioned matrices, however, are assumed to be positive definite (cf. Géradin and Rixen [25] or Balachandran and Magrab [27]). The vector $S(q,q\u02d9,t)$ includes nonlinear terms depending on positions and velocities as well as time varying terms such as external and/or parameteric excitation. Introducing the state-space coordinate $x:=[q\u22a4,q\u02d9\u22a4]\u22a4$, the notation

**B**as

where $x0$ denotes a deterministic initial condition.

## 3 Small Noise Expansion

Similar to approximate solutions obtained through perturbation analysis in the deterministic case (e.g., Verhulst [28] or Nayfeh [29]), approximate solutions to system (6) can be obtained through an expansion in the small parameter *σ*. The convergence of such a series can be guaranteed by the following result due to Blagoveshchenskii [30]:

*Assume that*$f(x,t)$

*and*$B(x,t)$

*and their partial derivatives with respect to the coordinate*

**x**

*up to order two are bounded and Lipschitz continuous for*$t0<t<t1$

*. Then, there exists*$\sigma 0>0$

*such that*

*holds, whereby the remainder*$r(t,\sigma )$

*is bounded in the mean square sense; that is*

*Proof.* This is a restatement of a Theorem by Blagoveshchenskii [30] for the current setting of the paper. Related versions have alo been stated by Freidlin and Wentzell [31] and Kaszás and Haller [32]. ◻

*Remark 3.1*. Polynomial nonlinearities appearing in common nonlinear oscillator prototypes such as Duffing's equation (cf. Kovacic and Brennan [33]) or the van der Pol oscillator (e.g., Guckenheimer and Holmes [34]) are neither bounded nor globally Lipschitz continuous. These nonlinearities generally grow unbounded for infinite states. Hence, Theorem 3.1 is not applicable in those cases. On the other hand, realistic systems, do not allow for arbitrarily large states. Thus, Eq. (1) is necessarily only valid for states inside a bounded domain; that is, the model's domain of validity. Thus, for example, one can follow an argument of Breunung [35] and modify Eq. (1) outside its domain of validity such that the requirements of Theorem 3.1 are met.

*σ*, the result is the following hierarchy of equations for the first two orders

*t*=

*t*

_{0}to their position at time

*t*is denoted by $Ft0t(x0)$. The gradient of $Ft0t(x0)$ with respect to the initial condition is denoted by $DFt0t(x0)$. The linearized flow map $DFt0t(x0)$ is also the fundamental solution to the deterministic part of the first order Eq. (10), i.e., it satisfies

which have been obtained by van Kampen [36]. Caughey [6] and Risken [4] present solutions for the special case of a constant Jacobian $\u2202xf(x0,t)=A1$ in Eq. (10).

**U**is orthogonal and the matrix $\Lambda $ is diagonal containing the $L\u2264N$ positive eigenvalues of $\Sigma (t;t0,x0)$. The pseudo-inverse of $\Lambda $ is defined elementwise by $\Lambda \u0303jl\u22121:=\delta jl/\Lambda ll$ for $1\u2264j,l\u2264L$ and $\Lambda \u0303jl\u22121:=0$ otherwise. With this notation, the pseudo-inverse of $\Sigma (t;t0,x0)$ is given by

which generally depends on the solution to the first order $x0$.

Thus, the PDF can be approximated by a Gaussian distribution whereby the mean is determined by solution to the deterministic limit $Ft0t(x0)$ and the variance is given in Eq. (12). The authors emphasize that both, the mean and variance, are defined through purely deterministic equations and hence, they can be computed without relying on stochastic integration methods.

is considered. First, the deterministic limit ($\sigma \u21920$) is analyzed. Varying the excitation frequency Ω in the vicinity of the natural frequency, the frequency response curve is computed with the numerical continuation package coco [37] and shown in Fig. 1(a). For the forcing frequency $\Omega =1.2$, two stable periodic responses exist for system (18). After selecting the initial condition $x0$ to be at the orbit with higher amplitude, the variances (12) are computed and depicted in Fig. 1(b). For comparison, 10^{4} approximations of solutions to the stochastic differential Eq. (18) are calculated with the Euler–Maruyama scheme (e.g., Kloeden and Platen [13]) with the noise intensity $\sigma =0.01$. Subsequently, the variances are calculated from the Euler–Maruyama samples.

From Fig. 1(b), it is evident that the approximation (16) is only accurate for about a quarter of the period *T*. After half a period, the computed variances (12) differ significantly from the results obtained through Monte Carlo sampling. This discrepancy is not an artifact of the obtained numerical approximation, as increasing the sample size or decreasing the step size in the Euler–Maruyama approximation does not alter Fig. 1(b). Furthermore, Fig. 1(b) is not in contradiction with Theorem 3.1, according to which, for any final time *t*_{1} there exists some *σ*_{0} such that the series (7) converges and the remainder is bounded in the mean square sense (cf. estimate (8)). It can be said that Fig. 1(b) solely indicates that for $t1=T$ and $\sigma =0.01$, the series (7) either does not converge or the remainder is unacceptably large. However, for $t1=T/4$, the approximation (16) is found to be acceptable. This observation is in line with experience from deterministic dynamical systems. In this setting, Nayfeh [29] has noted that the straightforward expansions in the form of Eq. (7) are often limited to small time intervals and more sophisticated methods need to be employed to obtain solutions valid over larger time intervals.

## 4 Approximation of Sample Paths

In the previous section, it was demonstrated that the explicit formula (16) works well for small enough times, but the formula's accuracy deteriorates significantly for longer time spans. To address this shortcoming, the small noise integrator (SNI) is proposed in Sec. 4.1. Subsequently, the integrator's performance is evaluated on nonlinear mechanical systems with up to one hundred degrees-of-freedom.

### 4.1 Small Noise Integrator.

To overcome the small time horizon of the approximation (16), the Algorithm 1 is proposed next. After fixing a time-step *τ* and providing an initial condition $x0$ at *t*_{0}, the distribution (15) for $t=\tau $ is computed and then this distribution is sampled to generate the new initial condition $x\tau $ according to Eq. (16). Then, the algorithm is restarted with the initial condition $x\tau $ to generate a new sample $x2\tau $. By repeating these two steps, samples are generated at later time instances $xk\tau $.

Result:$x(K\tau ;x0,t0)$ |

Set time step τ and provide initial condition $x0$; |

K = 1; |

whilek < Kdo |

Solve Eq. (9) to obtain the deterministic solution $F(k\u22121)\tau k\tau (xk\u22121)$; |

Compute the variance (12) ; |

Sample the distribution $p1(x,k\tau ;xk\u22121,(k\u22121)\tau )$ (cf. Eq. (15)) to obtain |

$x1(k\tau ;xk\u22121,(k\u22121)\tau )$; |

$xk=F(k\u22121)\tau k\tau (xk\u22121)+\sigma x1(k\tau ;xk\u22121,(k\u22121)\tau )$; |

$k=k+1$; |

end |

Result:$x(K\tau ;x0,t0)$ |

Set time step τ and provide initial condition $x0$; |

K = 1; |

whilek < Kdo |

Solve Eq. (9) to obtain the deterministic solution $F(k\u22121)\tau k\tau (xk\u22121)$; |

Compute the variance (12) ; |

Sample the distribution $p1(x,k\tau ;xk\u22121,(k\u22121)\tau )$ (cf. Eq. (15)) to obtain |

$x1(k\tau ;xk\u22121,(k\u22121)\tau )$; |

$xk=F(k\u22121)\tau k\tau (xk\u22121)+\sigma x1(k\tau ;xk\u22121,(k\u22121)\tau )$; |

$k=k+1$; |

end |

The SNI is motivated by the fact that the small noise expansion is accurate only for a limited time span. Thus, by selecting *τ* to be within this time span, one can ensure that the approximation (16) is accurate. Repeated sampling removes the limited time horizon of the small noise expansion and the following Lemma guarantees the accuracy of the SNI 1 for longer time spans:

*Assume that the flow map of the deterministic limit (9) is approximated with an accuracy not less than*$O(\tau )$

*. Then the SNI 1 strongly converges to the sample paths of the stochastic system (6) with order*$O(\tau 12)$

*; that is, the following holds*

*where*$x(k\tau )$*denotes a solution to system (6),*$xkSNI$*is an approximation obtained by the SNI 1, and C is a finite constant.*

*Proof.* The above Lemma 4.1 is proven by showing that for system (6) the SNI 1 converges to the Euler–Maruyama approximation for a small enough time-step *τ*. Then, the convergence follows from well-established convergence results (e.g., Kloeden and Platen [13] or Schurz's chapter in Ref. [17]). The details are presented in Appendix A. ▪

*Remark 4.1.* Notably, the convergence result of Lemma 4.1 does not depend on the size of the noise intensity *σ*. More specifically, in Appendix A only a parameter expansion in the time-step *τ* is performed and no restrictions on the noise intensity *σ* are imposed. For small enough time-step *τ*, the SNI 1 resembles the Euler–Maruyama scheme, which is valid for arbitrarily large *σ*. Hence, the SNI 1 will also converge for large *σ*.

*Remark 4.2.* As utilized in the Euler–Maruyama scheme and also shown by Risken [4], solutions to Eq. (6) can be approximated by a Gaussian process on small time scales regardless of the size of the noise intensity. With the SNI 1 one arrives at a similar conclusion and exploits the convergence of the expansion (7) for small noise intensities. Within the SNI 1 a non-Gaussian distribution arises on longer time scales when multiple sample paths are computed.

*N*integrations of the equations of variation, the new variable $V(t):=DFt0t(x0)B$ is introduced. Then, the variance (12) can be reformulated as

Thus, it suffices to compute the matrix **V** which requires only *M* integrations of the linearized flow.

The fact, that for small times solutions of system (6) can be approximated by a Gaussian process is well-known (see, e.g., Risken [4] or Freidlin and Wentzell [31]) and has also been exploited in the construction of solution schemes that repeatedly resample an arising Gaussian distribution in a similar manner as the SNI 1. For example, Sun and Hsu [38] introduce a short-time Gaussian approximation for the cell mapping method or Yu et al. [10] and Wehner and Wolfer [39] use a Gaussian transition probability density function in the formulation of a path integral solution to the Fokker–Planck equation. In these schemes, however, the mean and variance of the Gaussian approximation are obtained differently from the SNI 1. Wehner and Wolfer [39] utilize the short time propagator also used in the Euler–Maruyama scheme and Sun and Hsu [38] and Yu et al. [10] rely on moment equations with heuristic closure schemes. Contrarily, the SNI 1 utilizes the rigorously deduced mean and variance (16) for the arising Gaussian processes.

where $\Delta W$ denotes the increments of the *m*-dimensional of the normalized Wiener process. In this setting, the SNI 1 replaces the forward Euler scheme of the deterministic part in the Euler–Maruyama scheme ($xk\u22121e+\tau f(xk\u22121e,(k\u22121)\tau )$) by the deterministic flow map $F(k\u22121)\tau k\tau (xk\u22121)$ at each time-step. Moreover, instead of sampling the normal distribution with zero mean and variance $\tau B(xk\u22121e,(k\u22121)\tau )[B(xk\u22121e,(k\u22121)\tau )]\u22a4$ the variance (12) is sampled for the stochastic increment in the SNI 1 at each time-step.

Although the order of convergence of the SNI is only one-half, which is the order of the classic Euler–Maruyama scheme, the SNI 1 has inherent advantages making it an appealing alternative to existing techniques to solve Eq. (6) numerically. First, the SNI requires only sampling of normal distributions and hence it does not require sampling of intricate distributions as some competing method (e.g., Schurz's chapter in Ref. [17]). Second, the SNI allows for the use of deterministic integration routines. More specifically, the flow-map $Ft0t(x0)$ can be approximated by any deterministic numerical integrator. This is especially appealing in applications, where accurate and effective approximations of solutions to the deterministic limit of Eq. (6) already require sophisticated numerical integration routines. For example, the forward Euler approximation of the deterministic flow map in the Euler–Maruyama scheme or Milstein's scheme (e.g., Kloeden and Platen [13]) can be unstable for simple linear equations (see, e.g., Butcher [40]).

Most importantly, the SNI's performance is drawn from not only the limit of an asymptotically small time-step, but this algorithm is also accurate for large time steps, given that the noise terms are sufficiently small. Thus, it is anticipated that the time-step *τ* can be chosen considerably larger than in the available numerical integration routines for stochastic differential equations.

Overall, a possible increase in speed and accuracy gained by using a deterministic integrator for system (1) is counteracted by the necessity of computing the linearized flow map $DFt0t(x0)$ and the variance (12). In special cases, this burden can be reduced by computing the matrix **V** instead (cf. Eq. (22)). These observations lead to the conclusion that the SNI 1 is superior to other stochastic integration methods if the deterministic limit is challenging enough such that the speed up gained by the deterministic numerical integration routines is larger than the computational effort of computing the linearized flow map $DFt0t(x0)$ and the variance (12). In the next section, it is shown that this can be the case, through many examples.

### 4.2 Numerical Investigations.

As a benchmarking study, the SNI 1 is compared to Euler–Maruyama (EM) approximations of solutions to Eq. (6) (e.g., Kloeden and Platen [13]). Although many higher-order schemes are available (cf. Schurz's overview in Ref. [17]), as reported by Mannella [22], the performance gain can be insignificant at the expense of additional computational expense. Thus, for a first comparison, the widely-used Euler–Maruyama method is used^{2}^{,}^{3}. Within the SNI 1 the following numerical routines are used. Equations (9) and (10), and in turn Eq. (22), are solved using Matlab's ode45 routine. Moreover, the integral to compute the variance (12), respectively Eq. (21), is approximated by using the trapezoidal rule.

For an initial demonstration of the accuracy of the SNI, the Duffing system (18) with parameter values (19) and excitation frequency $\Omega =1.2$ is investigated for two different choices of the time-step *τ*. For the first choice, $\tau =T/5$ is selected, and subsequently, for the second choice, it is increased to *T* for comparison. The results of the SNI are shown in Fig. 2. While the distribution obtained with the SNI 1 with time-step $\tau =T/5$ shown in Fig. 2(a) matches the distribution from the EM-approximations, the SNI does not yield an accurate approximation for the time-step $\tau =T$ (cf. Fig. 2(b)). The accuracy can be quantified by computing the first two statistical moments for all three distributions. In Fig. 2(a), the mean and standard deviation of the EM-approximation and the SNI agree up to an accuracy of $10\u22123$, whereas the standard deviations differ by more than 0.01 in Fig. (2(b)). With the SNI 1, one computes the 10^{3} samples shown in Fig. 2 about 20-times faster than the Euler–Maruyama approximation.

In Fig. 3(a), the authors show the convergence of the SNI relative to the step size *τ*. The errors in the first two statistical moments, that is the mean and variance, computed with the SNI 1 are compared to those associated with from an Euler–Maruyama approximation with step size of $\tau =T/(2\xb7105)$. For the selected step sizes, the mean computed with the SNI 1 remains accurate with an error of the order of about $10\u22123$. This accuracy is also recognizable in Fig. 2. For both step sizes $\tau =T$ and $\tau =T/5$, the sample populations are centered at the same phase space location indicating that both distributions have a similar mean. However, the accuracy of sample variance computed with the SNI 1 increases significantly between $\tau =T/2$ and $\tau =T/3$. This observation is also discernible in Fig. 2. The sample population computed with the SNI 1 for the step size $\tau =T$ spans an evidently larger area than the converged sample distribution (cf. Fig. 2(b)). The same numerical convergence analysis is also performed for the Euler–Maruyama scheme and shown in Fig. 3(a). A converged result from the Euler–Maruyama scheme requires a significantly higher number of time steps per period compared to that for SNI 1.

*n*th mass is

*s*induce coupling between the individual degrees-of-freedom of system (24). The quantity

_{n}*a*is the amplitude of the deterministic harmonic component in the

_{n}*n*th forcing, and the quantity

*g*is the level of the noise term

_{n}*W*

_{1}in the

*n*th coordinate. It is noted that considering the oscillator array (24) with a single mass ($N*=1$) yields the Duffing Eq. (18) investigated in Figs. (2) and (3). Next, identical parameters for each oscillator are chosen with the following dimensionless values

and the noise intensity is set to $\sigma =0.01$.

The excitation frequency Ω is set to 1.2 and for each oscillator, the same initial condition $x0$ shown in Fig. 1(a) is selected. Then, the sample paths of the oscillator array (24) are approximated for 1, 10, and 100 forcing periods. To ensure that the Euler approximation yields the correct deterministic limit, the step size of the EM-approximation is decreased until it matches the solution of Matlab's ode45 with an error less than $10\u22123$. For the SNI 1, $\tau =T/5$ is selected and Matlab's ode45 algorithm is used to solve Eqs. (6) and (22). For both approximations, 10^{3} sample paths are obtained and the first two statistical moments are compared. It is observed that in all cases there is an agreement with an accuracy of $10\u22123$. The run-time comparison of the EM-approximation and the SNI 1 is shown in Table 1.

Final time t = T | Final time $t=10T$ | Final time $t=100T$ | ||||
---|---|---|---|---|---|---|

# dof ($N*$) | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up |

1 | $4\xd7104$ | 21 | $6\xd7104$ | 32 | $1\xd7105$ | 47 |

2 | $6\xd7104$ | 32 | $8\xd7104$ | 45 | $2\xd7105$ | 95 |

5 | $8\xd7104$ | 37 | $1\xd7105$ | 52 | $4\xd7105$ | 194 |

10 | $1\xd7105$ | 46 | $2\xd7105$ | 195 | $6\xd7105$ | 259 |

20 | $2\xd7105$ | 72 | $3\xd7105$ | 115 | $6\xd7105$ | 175^{a} |

50 | $3\xd7105$ | 46 | $4\xd7105$ | 45 | $8\xd7105$ | 102^{a} |

100 | $4\xd7105$ | 21 | $6\xd7105$ | 32 | $1\xd7106$ | 62^{a} |

Final time t = T | Final time $t=10T$ | Final time $t=100T$ | ||||
---|---|---|---|---|---|---|

# dof ($N*$) | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up |

1 | $4\xd7104$ | 21 | $6\xd7104$ | 32 | $1\xd7105$ | 47 |

2 | $6\xd7104$ | 32 | $8\xd7104$ | 45 | $2\xd7105$ | 95 |

5 | $8\xd7104$ | 37 | $1\xd7105$ | 52 | $4\xd7105$ | 194 |

10 | $1\xd7105$ | 46 | $2\xd7105$ | 195 | $6\xd7105$ | 259 |

20 | $2\xd7105$ | 72 | $3\xd7105$ | 115 | $6\xd7105$ | 175^{a} |

50 | $3\xd7105$ | 46 | $4\xd7105$ | 45 | $8\xd7105$ | 102^{a} |

100 | $4\xd7105$ | 21 | $6\xd7105$ | 32 | $1\xd7106$ | 62^{a} |

The computations have been performed using matlab 2020a installed on a Windows PC with Intel Xeon CPU E5-2687 W @ 3.1 GHz and 64 GB RAM.

The sample size is 10^{3} and the convergence of the two methods is ensured by computing the first two statistical moments. The excitation frequency Ω is 1.2.

Due to the excessive computation time of the Euler–Maruyama scheme (more than one day), the run-time for the full sample size is estimated by extrapolating the run-time for one sample.

Overall, from Table 1, it can be discerned that a significant speed-up of the SNI 1 is obtained in comparison to the classical Euler–Maruyama scheme. In all cases, the speed-up is at least one order of magnitude. Especially large is the speed-up for oscillator arrays with 5 to 50 oscillators. Moreover, the efficiency of the SNI 1 increases with the time span, which makes it an appealing choice for long time horizons.

To further demonstrate the versatility of the SNI 1, the excitation frequency is increased to $\Omega =1.9$ and for each oscillator, the initial condition is selected to be on the upper stable branch depicted in Fig. 1(a). At this frequency, Matlab's ode45 algorithm can be used to compute the stable periodic orbit effectively, whereas the step size of the primitive EM-approximation needs to be decreased significantly to converge to a stable periodic orbit of the deterministic limit. Proceeding as previously described for the stochastic simulations, the SNI 1 yields a significant computational gain compared to the Euler–Maruyama scheme (cf. Table 2).

$\Omega =1.9$ | Random parameters | |||||
---|---|---|---|---|---|---|

Final time t = T | Final time t = T | Final time $t=10T$ | ||||

# dof ($N*$) | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up |

1 | $2\xd7105$ | 87 | $6\xd7104$ | 21 | $3\xd7105$ | 89 |

2 | $3\xd7105$ | 130 | $2\xd7104$ | 16 | $2\xd7105$ | 161 |

5 | $5\xd7105$ | 156 | $2\xd7104$ | 13 | $2\xd7105$ | 163 |

10 | $8\xd7105$ | 303 | $5\xd7104$ | 27 | $5\xd7105$ | 270 |

20 | $1\xd7106$ | 257 | $5\xd7105$ | 60 | $4\xd7106$ | 488^{a} |

50 | $2\xd7106$ | 270 | $5\xd7106$ | 22 | $1\xd7107$ | 250^{a} |

100 | $3\xd7106$ | 153 | $5\xd7105$ | 7 | $5\xd7106$ | 201^{a} |

$\Omega =1.9$ | Random parameters | |||||
---|---|---|---|---|---|---|

Final time t = T | Final time t = T | Final time $t=10T$ | ||||

# dof ($N*$) | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up | EM: Steps per T | Speed-up |

1 | $2\xd7105$ | 87 | $6\xd7104$ | 21 | $3\xd7105$ | 89 |

2 | $3\xd7105$ | 130 | $2\xd7104$ | 16 | $2\xd7105$ | 161 |

5 | $5\xd7105$ | 156 | $2\xd7104$ | 13 | $2\xd7105$ | 163 |

10 | $8\xd7105$ | 303 | $5\xd7104$ | 27 | $5\xd7105$ | 270 |

20 | $1\xd7106$ | 257 | $5\xd7105$ | 60 | $4\xd7106$ | 488^{a} |

50 | $2\xd7106$ | 270 | $5\xd7106$ | 22 | $1\xd7107$ | 250^{a} |

100 | $3\xd7106$ | 153 | $5\xd7105$ | 7 | $5\xd7106$ | 201^{a} |

The computations have been performed using MATLAB 2020a installed on a Windows PC with Intel Xeon CPU E5-2687 W @ 3.1 GHz and 64 GB RAM.

(i) $\Omega =1.9$ and uniform parameters (25) and (ii) random parameters (cf. also Appendix B). The sample size is 10^{3} and the convergence of the two methods is ensured by computing the first two statistical moments.

Due to the excessive computation time of the Euler–Maruyama scheme (more than one day), the run-time for the full sample size is estimated by extrapolating the run-time for one sample.

To verify that the performance of the SNI 1 is independent of specific parameter values (25), the parameter values of system (24) are assigned randomly as described in Appendix B. Also for this numerical experiment the SNI 1 is found to outperform the Euler–Maruyama approximation (cf. Table 2). Overall, from Table 2, one can discern a speed-up of about two orders in magnitude. Moreover, for longer simulation times an even higher speed-up of the SNI 1 is expected, since the small step size to be used for the EM-approximation will increase the computational burden excessively.

Although the SNI 1 relies on a repeated sampling of Gaussian distributions, the resulting sample distribution is not necessarily Gaussian. To demonstrate this, the oscillator array (24) with two masses $N*=2$ and a coupling spring stiffness $kc=0.025$ is considered. For these parameter values, four stable periodic orbits exist. For one of there orbits, both masses oscillate with a high amplitude (i.e., the high amplitude orbit), whereas along the low amplitude orbit, the displacements of both masses are low. Moreover, two *localized modes* exist. Along one of the localized modes, the first mass oscillates considerably, whereas the amplitude of the second mass is low. For the other localized mode, the energy distribution is reversed; that is, the amplitude of the first mass is low and the second mass oscillates with a high amplitude. Such energy localization phenomena have been previously investigated, for example, Sievers and Takeno [41], Vakakis and Cetinkaya [42] and Dick et al. [43]. Recently, there has been an interest to study the effects of noise on such energy localization (e.g., Perkins et al. [44] and Balachandran et al. [45]). Selecting the high amplitude orbit as an initial condition, 10^{3} sample paths for 100 periods have been computed with the SNI 1. The sample locations are shown in Fig. 5. The projections show that some samples remained in the neighborhood of the high amplitude orbit, whereas other realizations have escaped toward the low amplitude periodic orbit. Overall, the distribution shown in Fig. 5 is clearly non-Gaussian.

## 5 Approximation of the Probability Density Function

A straightforward method to compute probability densities is to discretize the state space into volumes, compute sample paths (e.g., with the SNI 1 or the EM-approximation) and then count the number of realizations within each volume. The result of this Monte Carlo approach is an invariably nonsmooth approximation of the probability density function. This requires a high number of samples to accurately represent the PDF. To overcome these shortcomings a more effective method is proposed in Sec. 5.1. This approach relies on the rigorously deduced Gaussian kernel (12) to smoothly approximate the probability density function. Next, the performance of the proposed method is examined through examples in Sec. 5.2.

### 5.1 Gaussian Kernel Approximation.

*R*realizations $xr(K\tau )=xr(K\tau ;x0,t0)$ of the stochastic process (6) have been computed. At time $t=(K\u22121)\tau $, the sample locations are given by $xr((K\u22121)\tau )=xr((K\u22121)\tau ;x0,t0)$. Now advancing the sample population from time $t=(K\u22121)\tau $ to the final time $t=K\tau $, each individual sample results in the generation of the Gaussian distribution

that is, a Gaussian kernel approximation (GKA) to the probability density at the final time. Equation (27) is an approximation of the probability density $p(x,K\tau ;x0,t0)$ obtained by a sum of Gaussian distributions centered at $F(N\u22121)\tau N\tau (xr((K\u22121)\tau ))$ with variance (12), which generally differ for each realization *r*. Hence, it can also be viewed as kernel smoothening method (e.g.,, Friedman et al. [46]). Compared to the usual adhoc chosen kernel distributions, the kernel (27) is rigorously chosen based on the Gaussian kernel (12). Moreover, since (27) is the standard Monte Carlo estimator, the corresponding convergence can be guaranteed under appropriate assumptions; that is, small enough *τ* and small enough *σ*.

### 5.2 Numerical Investigations.

In the following, the time-varying probability density function of the stochastic process (6) is computed with the Gaussian kernel approximation (27) and compared with the probability density obtained by using the straight forward Monte Carlo simulations, wherein the state space is discretized into volumes and the number of realizations in each volume is counted.

For a first comparison, the single degree-of-freedom oscillator (18) with the parameters (19) is reinvestigated. The initial condition is depicted in Fig. 1(a) and the solutions of the stochastic process (18) are approximated via the SNI 1 for one period. In Fig. 6, the authors show the obtained probability density functions for various sample sizes. For 10^{5} samples, the probability density function of both methods, Monte Carlo and GKA (27) match very well. Reducing the sample size in the GKA to include only 10^{3} samples, still yields a PDF, which matches the PDF obtained at convergence with 10^{5} samples. However, with a reduced number of samples in the Monte Carlo method, there is a significant deviation from the converged results.

To quantify the superior convergence of the GKA, the *L*_{2}-error relative to Monte Carlo simulations with 10^{6} samples is visualized for both approximations and various sample sizes. The relative error of the GKA shown in Fig. 7 is about two orders of magnitude smaller than that of the Monte Carlo approximation. Only for the largest sample size ($Ns=105$), the error of the GKA is one order of magnitude less than that of the Monte Carlo simulations. However, the error with the Monte Carlo simulations with 10^{5} samples is still larger than the error with the GKA with only 10^{3} samples. Thus, in this example, the sample size of the GKA-approximation can be reduced by two orders of magnitude compared to the Monte Carlo method.

Before proceeding with a high dimensional oscillator array (24), it is insightful to emphasize the enormous sample size needed for the Monte Carlo simulations to obtain a converged PDF. For the one degree-of-freedom oscillator considered in Figs. 6 and 7, about 10^{5} samples were necessary to compute a converged PDF with the Monte Carlo method. Thus, an estimated number of $105N*$ samples would be required to obtain a converged PDF for an $N*$-oscillator array (24). Even for three oscillators in system (24), the computational effort necessary to generate such a high number of samples is overwhelming^{4}.

To verify whether one obtains the correct results with the GKA-approximation (27) in higher dimensions, an perturbation scheme is employed. For the weak coupling spring stiffness *s _{n}*, the time evolution the oscillator array (24) can be estimated with

*N*uncoupled oscillators; that is, the

*uncoupled*limit $sn\u21920$. More precisely, it is assumed that the coupling spring stiffness

*s*is of order $O(\sigma 2)$. Then, for the uniform parameters (25) each oscillator

_{n}*q*and $q\u02d9j$ are identical up to order $O(\sigma 2)$. Thus, in this limit, it is sufficient to simulate a single degree-of-freedom oscillator to infer about the full oscillator array (24).

_{j}In Fig. 8, the authors depict the PDF of the oscillator array (24) with three masses *N *=* *3 along the two-dimensional plane where the positions and velocities of the second and third mass are at their sample mean; that is, $q2=q\xaf2,\u2009q3=q\xaf3,\u2009q\u02d92=q\u02d9\xaf2$ and $q\u02d92=q\u02d9\xaf2$, where the bar indicates the mean values. Both approximations match with an relative *L*_{2}-error of less than 0.01. It is emphasized, that for the Monte Carlo simulations approximate symmetries of system (24) are explicitly exploited, whereas no such reduction technique has been employed for the GKA^{5}. Similar to the single degree-of-freedom case, with the GKA (27), one converges to a significantly lower number of samples than with the Monte Carlo sampling.

Finally, the two-degree-of-freedom system with localized modes from Sec. 4.2 is reinvestigated to verify whether samples escaped from the high amplitude orbit to the other stable periodic orbits. With the two-dimensional projection shown in Fig. 5, one invariably ignores the other directions and hence such projections are of limited use to infer about an escape from the high amplitude orbit. Accordingly, measuring distances of samples to the periodic orbits is more appropriate. For each periodic orbit, the distance is a random variable, whose PDF can be approximated with the GKA. The arising four PDFs are shown in Fig. 9 for 10^{5} samples. The peak in Fig. 9 is due to the initial condition, whose distance to the localized modes and the low amplitude orbit is about 2.7. The distances to the high amplitude orbit and the low amplitude orbit (blue and red lines in Fig. 9) confirm the impression from Fig. 5 that realizations do not stay close to the high amplitude orbit and escape toward the low energy orbit.

## 6 Conclusions

In this work, the small noise expansion (cf. Eq. (7)) is exploited to propose a method to obtain approximations for the following: (i) the sample paths of the stochastic dynamical system (1) and (ii) the associated time-varying probability density function. With the formulated small noise integrator 1, one removes the limited time horizon of the small noise expansion, by using an appropriate resampling of the Gaussian distribution (12) along sample paths. For the SNI 1, one only requires the solution to deterministic differential equations, which means that deterministic integration routines can be used in the stochastic setting. Additionally, the authors have proven convergence of the proposed algorithm, which notably also holds for arbitrarily large noise intensities.

The computational benefit of the proposed SNI is examined with a series of coupled oscillator arrays with up to one hundred degrees-of-freedom, including a randomly parameterized array. Compared to the standard Euler–Maruyama scheme, a speed-up of two orders of magnitude is observed for a wide range of investigated systems. Especially for longer simulation times, the computational gain of the SNI is significant.

Moreover, the Gaussian kernel approximation (27) yields a justified approach to recover a smooth probability density function from samples, without relying on adhoc kernel choices or interpolation. With this method, the necessary sample size to accurately approximate the probability density function can be reduced drastically. For a single degree-of-freedom oscillator, the number of samples can be reduced by three orders of magnitude compared to established Monte Carlo methods. For a three degree-of-freedom system, the GKA (27) is found to yield an accurate PDF, while a computation with Monte Carlo methods is simply infeasible with reasonable computational resources. Thus, the GKA (27) opens up a new horizon to compute PDFs for higher dimensional systems beyond the current limitation to two or three dimensions.

While the SNI 1 can be used to extend the validity of the straightforward expansion (7) to longer time intervals, it does not yield a steady-state distribution, which can be of interest in applications. The underlying theory 3.1, unfortunately, is fundamentally restricted to finite time intervals. Thus, an extension of the underlying theory as well as computational algorithms is desirable to compute a statistical steady-state, if it exists.

The probability density function, a time-varying scalar quantity in a usually high dimensional space, is often difficult to understand, visualize or access. In many applications, the quantities of interest, for example, the exceedance probability, escape times, or frequency of occurrence, can, in principle be answered by computing the PDF, but this might not always be the most efficient approach. Thus, it would be of interest to tailor the method presented in such specific contexts.

## Acknowledgment

The authors gratefully acknowledge the support of the National Science Foundation through Grant No. CMMI 1760366 and the associated data science supplements. Additionally, T.B. is thankful to Lautaro Clienti and Abdulrahman Alofi for encouraging discussions throughout this work.

## Funding Data

National Science Foundation (Division of Civil, Mechanical and Manufacturing Innovation) (Grant No.: CMMI 1760366; Funder ID: 10.13039/100000147).

## Data Availability Statement

The data that supports the findings of this study is openly available online.^{6}

## Conflict of Interest

The authors have no conflicts to disclose.

## Footnotes

The Matlab provided Euler-Maruyama approximation is not utilized. Instead, a self-written code, available at github.com/tbreunung/SNI, is employed. For the cases tested the self-written code outperforms the Matlab provided code. Moreover, various optimizations of the Euler-Maruyama scheme summarized by Higham [47] (especially vectorization) have been tested and a version was selected that seemed optimal in terms of speed and memory requirements for the purposes of this study.

It is noted that for the investigated systems, the first order equivalent of systems (18) and (24), the matrix B (cf. Equation (6)) does not depend on the coordinates. Hence, the Euler-Maruyama scheme is equivalent to Milstein's scheme (e.g., Kloeden and Platen [13]). Thus, it convergences to the sample path of system (24) with order one.

To give an idea of the computational effort, the following is mentioned: the storage required to store 10^{15} samples with double precision floating point numbers is 48 petabytes ($=1015\xb73\xb72\xb78$ bytes). With the current advances in high-performance computing and data storage systems, such simulations may not be impossible but one requires a serious commitment in computing power and data storing capabilities, which is beyond the scope of this paper.

Without using the approximate symmetries, not a single realization out of the 10^{5} samples ended up in the discretized state space volume used to compute the results for Fig. 8. Hence, the Monte Carlo approximation for the PDF is identically zero.

### Appendix A: Convergence Proof of the SNI 1

*τ*is given by

*m*-dimensional of the normalized Wiener process; that is, they are independent and identically normal-distributed random variables with zero mean and variance one. First, the variance (12) is expanded for small

*τ*yielding

*M*columns of $\tau B(xk\u22121,(k\u22121)\tau )$ by samples drawn from the standard distribution. Thus, one time-step of the SNI can be written as

where the following is used: (i) the convergence result of the stochastic Euler–Maruyama approximation (A1) to the solutions of Eq. (10) (e.g., Kloeden and Platen [13]), (ii) the assumption that the approximated flow map $F\u0303k\tau (n\u22121)\tau (xkSNI)$ is $O(\tau )$ close to the true flow map $F\u0303k\tau (n\u22121)\tau (xkSNI)$ (cf. assumption in Lemma 4.1), and iii) the convergence of the deterministic forward Euler scheme. The required estimate for Lemma 4.1 is given in Eq. (A4).

### Appendix B: Random Parameter Selection for Oscillator Array (24)

*μ*and variance $\sigma N$ or the uniform distribution $U[a,b]$ with

*a*denoting the minimum and

*b*the maximum value are sampled. The parameter values are drawn from the following distributions:

The distributions (A6) are selected such that the linear unforced limit of the oscillator array (24) ($\kappa \u21920$ and $\sigma \u21920$) is a weakly damped oscillator customary in the structural dynamics literature (cf. Géradin and Rixen [25]). The parameter values for the numerical experiments presented in Table 2 are available online.^{6}