Abstract
A technique based on the Wiener path integral (WPI) is developed for determining the stochastic response of diverse nonlinear systems with fractional derivative elements. Specifically, a reduced-order WPI formulation is proposed, which can be construed as an approximation-free dimension reduction approach that renders the associated computational cost independent of the total number of stochastic dimensions of the problem. In fact, the herein developed technique can determine, directly, any lower-dimensional joint response probability density function corresponding to a subset only of the response vector components. This is done by utilizing an appropriate combination of fixed and free boundary conditions in the related variational, functional minimization, problem. Notably, the reduced-order WPI formulation is particularly advantageous for problems where the interest lies in few only specific degrees-of-freedom whose stochastic response is critical for the design and optimization of the overall system. An indicative numerical example is considered pertaining to a stochastically excited tuned mass-damper-inerter nonlinear system with a fractional derivative element. Comparisons with relevant Monte Carlo simulation data demonstrate the accuracy and computational efficiency of the technique.
1 Introduction
Various methodologies have been developed over the last six decades in the field of stochastic engineering dynamics for solving the governing stochastic differential equations of motion and for determining response statistics of diverse structural and mechanical systems; see Refs. [1–4] for a broad perspective. Nevertheless, in many cases these methodologies are limited to treating only conventional forms of the equations of motion based on classical mechanics theories. In fact, generalizing the aforementioned solution approaches to account for systems exhibiting time- and space-localized behaviors described by operators based on wavelets and/or noninteger order (fractional) derivatives can be a rather challenging task (e.g., Refs. [5–7]).
Specifically, fractional calculus can be construed as a generalization of ordinary calculus (e.g., Ref. [8]) that has been successfully employed, indicatively, for developing nonlocal continuum mechanics theories (e.g., Refs. [9,10]), and for modeling viscoelastic materials (e.g., Ref. [11]). Further, a wide range of methodologies, including extensions of the standard statistical linearization solution scheme [2], have been developed over the past few years for determining the stochastic response of nonlinear systems endowed with fractional derivative terms (e.g., Refs. [12–16]). Note, however, that in many cases the performance of the aforementioned methodologies is hindered by the associated significant computational cost and/or the quite high degree of the involved approximations. For example, a statistical linearization solution treatment was proposed in Ref. [16], where the original single-degree-of-freedom (single-DOF) oscillator with a fractional derivative element was cast, equivalently, into a multi-DOF system with integer-order derivatives only. This was done by augmenting the response vector and by considering additional auxiliary state variables and equations; see also Ref. [17]. Inevitably, the increased dimensionality of the problem due to the augmentation of the response vector translated into nontrivial additional computational cost. Further, due to the Gaussian response assumption [2], a statistical linearization solution treatment exhibits satisfactory accuracy in determining only first- and second-order response statistics (i.e., mean vector and covariance matrix), and not the complete joint probability density function (PDF).
Over the past few years, Kougioumtzoglou and coworkers have developed a Wiener path integral (WPI) technique for determining the stochastic response of diverse nonlinear dynamical systems (e.g., Refs. [18–20]). The fundamental concept of the WPI technique, which exhibits both low computational cost and high accuracy, relates to treating the system response joint transition PDF as a functional integral over the space of all possible paths connecting the initial and the final states of the response vector. Further, the functional integral is evaluated, ordinarily, by resorting to an approximate approach that considers the contribution only of the most probable path. This corresponds to an extremum of the functional integrand and is determined by solving a functional minimization problem that takes the form of a deterministic boundary value problem (BVP).
Notably, the WPI technique has been extended recently to treat systems endowed with fractional derivative elements; see Refs. [21–23]. Specifically, a single-DOF nonlinear oscillator with a fractional derivative element was considered in Ref. [21] that was cast, equivalently, into a multi-DOF nonlinear system with integer-order derivatives only. Next, to circumvent the challenge of increased dimensionality of the problem due to the augmentation of the response vector, an efficient reduced-order WPI variational formulation was employed that was originally developed in Ref. [20]. This reduced-order formulation can be construed as an approximation-free dimension reduction approach that renders the associated computational cost independent of the total number of stochastic dimensions of the problem. Concisely stated, the technique in Ref. [20] can determine directly any lower-dimensional joint response PDF corresponding to a subset only of the response vector components by utilizing an appropriate combination of fixed and free boundary conditions. Thus, the original single-DOF oscillator joint PDF corresponding to the response displacement and velocity was determined efficiently, while circumventing the rather challenging task of treating directly equations of motion involving fractional derivatives. Further, a conceptually different and more direct approach was proposed in Ref. [22], where the derived closed-form expression for the WPI-based joint response PDF involved the fractional derivative terms of the original equation of motion. This led to a fractional variational problem that was solved numerically based on a Rayleigh-Ritz scheme to determine the most probable path, and ultimately, to obtain the system response joint PDF. Furthermore, based on a functional change of variables treatment, it was proved in Ref. [23] that the mathematical formulation of Ref. [22] can be generalized to treat a broader class of systems with non-Markovian response processes. In fact, it was shown that the WPI technique can treat nonlinear systems with a history-dependent state, such as hysteretic structures or oscillators endowed with fractional derivative elements, in a straightforward manner; that is, without resorting to ad hoc modifications pertaining, typically, to employing additional auxiliary filter equations and state variables.
In this paper, the reduced-order WPI formulation of Ref. [20] is extended to account for systems endowed with fractional derivative elements. Alternatively, the herein developed technique can be construed as a generalization of the results in Ref. [22] to account for mixed boundary conditions in the variational formulation. An indicative numerical example is considered pertaining to a stochastically excited tuned mass-damper-inerter nonlinear system with a fractional derivative element. Comparisons with relevant Monte Carlo simulation (MCS) data are included as well for demonstrating the accuracy and computational efficiency of the reduced-order formulation.
2 Preliminaries
2.1 Stochastic Differential Equations of Motion.
In this section, some basic aspects of the theory of stochastic differential equations (SDEs) are presented for completeness. The interested reader is also directed to Refs. [24] and [25] for a broader perspective.
where the dot above a variable denotes differentiation with respect to time t, and denotes a zero-mean and delta-correlated Gaussian white noise stochastic process of intensity one. That is, and for any , where δjk is the Kronecker delta and is the Dirac delta function. Regarding the relation between the Wiener and the white noise processes, can be defined as an infinitesimal jump of the Wiener process, i.e.,, . Thus, it is often, informally, written as the time derivative of the Wiener process in the form ; see also Refs. [24] and [26] for a more detailed discussion on the topic.
where the square root of matrix is given by .
2.2 Wiener Path Integral Formalism.
In this section, the salient aspects of a recently developed novel formalism of the WPI technique for determining the stochastic response of diverse dynamical systems are presented concisely for completeness; see Ref. [23] for more details. In comparison to alternative derivations in the literature, which resort to the Chapman–Kolmogorov equation as the starting point, the novel formalism circumvents the Markovian assumption for the system response process. In this regard, nonlinear systems with a history-dependent state, such as hysteretic structures or oscillators endowed with fractional derivative elements, can be treated in a direct manner; that is, without resorting to ad hoc modifications of the WPI technique pertaining, typically, to employing additional auxiliary filter equations and state variables (e.g., Ref. [21]).
Next, taking into account the relationship between the system response process and the Wiener process W described by Eq. (1), the transition PDF is obtained by applying a functional change of variables to Eq. (8); see Ref. [23] for more details.
2.3 Most Probable Path Approximation.
It is remarked that the analytical evaluation of the WPI of Eq. (13) is, in general, an impossible task. Thus, alternative approaches are typically pursued in the literature for evaluating approximately Eq. (13), such as the most probable path approach (e.g., Ref. [27]). Note that the most probable path approximation has exhibited a quite high degree of accuracy in various diverse engineering mechanics applications (e.g., Refs. [28] and [29]).
In passing, it is remarked that in cases where an even higher degree of accuracy is required (e.g., reliability assessment applications), fluctuations around the most probable path can be also accounted for (e.g., Ref. [19,30]). Further, various methodologies can be employed for treating the optimization problem of Eq. (17) and for determining . These range from standard Rayleigh-Ritz type numerical solution schemes (e.g., Refs. [28,31]) to more recently developed techniques relying on computational algebraic geometry tools [32]. Alternatively, considering Eq. (17) and resorting to calculus of variations yields the corresponding Euler–Lagrange equations, which take the form of a BVP to be solved for obtaining the most probable path (e.g., Ref. [33]).
3 Mathematical Formulation
In this section, the WPI technique is generalized by developing a variational formulation with mixed fixed/free boundary conditions for determining the stochastic response of diverse nonlinear systems endowed with fractional derivative elements. The herein developed technique can be construed as an extension of the methodology proposed in Ref. [20] to account for fractional derivatives in the system equations of motion.
In Eq. (20), U0 and U1 are sets containing the indices of the coordinates of x and , respectively, that participate in the lower-dimensional joint response PDF. In this regard, the boundary conditions at tf are considered fixed for these indices, while the rest are treated as free. Notably, the associated variational formulation of the WPI technique allows for determining the dimensional, marginalized, joint response PDF in a direct manner at a reduced computational cost. In fact, the cost depends only on the reduced dimension p of the target PDF and not on the dimension 2 m of the original system.
with denoting the Gamma function.
3.1 Variational formulation and Euler–Lagrange Equations With Mixed Fixed/Free Boundary Conditions.
It is remarked that the herein developed technique is an extension of the reduced-order WPI formulation in Ref. [20] to account for systems endowed with fractional derivative elements. Alternatively, the technique can be construed as a generalization of the results in Ref. [22] to account for mixed boundary conditions in the variational formulation. Indeed, it can be readily seen that for the special case of all boundary conditions being fixed, i.e., , the technique reduces to the formulation in Ref. [22].
3.2 Computational Aspects.
Truncating the expression in Eq. (40) after a number of terms and substituting into the Lagrangian of Eq. (24) yields a BVP to be solved for the most probable path that involves integer order derivatives only. This can be readily solved by a plethora of well-established numerical schemes (e.g., Ref. [43]).
Further, as noted in Sec. 3.1, for the special case where , all final conditions at tf in Eq. (38) become fixed, and Eq. (39) corresponds to the complete dimensional joint PDF. In this regard, for a given time instant tf, a standard brute-force numerical implementation of the technique entails the discretization of the PDF effective domain into points, where N is the number of points in each dimension, and the evaluation of the PDF is performed point-wise on the discretized lattice. In other words, BVPs in the form of Eqs. (37) and (38) need to be solved numerically, yielding an exponential increase of the computational cost with the dimension m of the system.
Remarkably, the variational formulation with mixed fixed/free boundaries developed in Sec. 3.1, which extends the technique in Ref. [20] to account for fractional derivative elements, can reduce the associated computational cost drastically by determining directly any lower-dimensional () joint response PDF. This capability is particularly advantageous for problems where the interest lies in few only specific DOF whose stochastic response is critical for the design and optimization of the overall DOF system. In this regard, the BVPs required to be solved by the standard formulation of the technique decrease to Np problems only, where is a reasonable range of values for various diverse engineering dynamics applications (e.g., Refs. [22,23,29]). Note that for small values of p relative to 2 m, it has been shown in Ref. [20] that the reduced-order WPI technique becomes orders of magnitude more efficient than an alternative MCS solution treatment.
3.3 Mechanization of the Numerical Implementation.
Concisely stated, the mechanization of the numerical implementation of the herein developed technique comprises the following steps:
For a given time instant tf, select sets U0 and U1 corresponding to the system response variables of interest.
Discretize the dimensional domain, where , with N points in each dimension.
For each point of the discretized domain, determine the most probable path by solving numerically Eq. (37) in conjunction with the boundary conditions of Eq. (38).
Obtain the dimensional joint response PDF by employing Eq. (39).
4 Numerical Example
An indicative numerical example pertaining to a stochastically excited tuned mass-damper-inerter (TMDI) nonlinear system with a fractional derivative element is considered next for demonstrating the efficiency and accuracy of the herein developed reduced-order WPI technique.
Specifically, TMDI systems can be construed as generalizations of the classical tuned mass-damper (TMD) concept for passive vibration control applications. In fact, TMDI configurations have demonstrated enhanced performance compared to standard TMD systems by exploiting the mass amplification effect of the inerter, a two-terminal flywheel device developing forces proportional to the relative acceleration of its terminals. The interested reader is directed to the recent review papers [44] and [45] for a broad perspective.
Notably, the majority of research efforts for the optimization of stochastically excited TMDI structures refer to linear systems for which convenient closed-form input–output relationships exist in the time or frequency domains (e.g., Refs. [46,47]). However, in many structural engineering applications, the primary structure (or the combined system) exhibits nonlinear/hysteretic response behavior (e.g., Refs. [45,48,49]). In such cases, employing a standard MCS-based solution treatment of the equations of motion translates into performing a relatively large number of computationally demanding time-domain response time-history analyses. In fact, the computational cost as a function of the number of MCS-based response analyses increases significantly and can become prohibitive when constraints pertaining to low probability events (e.g., failures) are considered in the TMDI optimization problem. Further, alternative, computationally efficient approximate techniques, such as statistical linearization [2], exhibit typically a low degree of accuracy and are incapable of treating optimization problems with low probability constraints. Indeed, as shown in Refs. [29] and [50] where a class of nonlinear electromechanical harvesters was considered, statistical linearization was strikingly outperformed in terms of accuracy by a WPI-based solution treatment in determining the system response PDF.
where the values , a = 0.5 and are used for the parameters in Eqs. (42)–(45). The restoring force of the primary structure exhibits cubic stiffness nonlinearities, whereas the secondary system is connected to the primary structure with a viscoelastic element whose force is modeled by a fractional derivative term as in Ref. [51].

Stochastically excited tuned mass-damper-inerter (TMDI) nonlinear system with a fractional derivative element
Next, in the ensuing implementation of the reduced-order WPI technique, only the response displacement of the primary structure, i.e., x1, is considered fixed at the final time instant tf. Thus, according to Eq. (39), the one-dimensional response PDF with fixed zero initial conditions at ti is determined in a direct manner by solving only N = 41 BVPs of the form of Eqs. (37) and (38). For the numerical solution of the BVPs, the approximation of Eq. (40) is adopted and the expansion is truncated at n = 2. Also, the initial time instant is set at , and results are obtained and plotted for the relative final time instant . In particular, Fig. 2 shows the evolution in time of the marginal response displacement PDF . Comparisons with MCS-based estimates (50,000 realizations) demonstrate a satisfactory degree of accuracy. For the MCS analyses, the L1 algorithm (e.g., Ref. [52]) has been used for integrating numerically Eq. (21) and for determining response realizations. Further, the quite high accuracy exhibited by the WPI technique is also shown in Fig. 3, where the response displacement PDF at various time instants is compared with pertinent MCS-based results.

Evolution in time of the response displacement PDF of the primary structure of a stochastically excited TMDI nonlinear system with a fractional derivative element: (a) WPI-based estimates and (b) MCS data—50,000 realizations

Wiener path integral (WPI)-based response displacement PDF of the primary structure of a stochastically excited TMDI nonlinear system with a fractional derivative element at various time instants; comparisons with MCS data—50,000 realizations
Remarkably, compared to a standard implementation of the WPI technique with all boundary conditions fixed, which requires BVPs to be solved and yields the four-dimensional joint PDF , the proposed reduced-order implementation minimizes the associated computational cost by obtaining the marginal PDF directly by solving only N = 41 BVPs. Notably, for this particular example, MCS based on 50,000 realizations requires approximately 15 min of computation time, whereas the marginal PDF of x1 at a specific time instant is determined via the proposed reduced-order WPI technique in approximately 4 s on the same computer.
Lastly, it can be argued that the reduced-order WPI technique is ideally suited for the considered engineering dynamics problem of vibration control, where the fundamental objective relates to determining and controlling, specifically, the DOF corresponding to the response displacement of the primary structure.
5 Concluding Remarks
In this paper, a reduced-order WPI technique has been developed for determining the stochastic response of nonlinear systems with fractional derivative elements. The technique can be construed as an extension of the formulation in Ref. [20] to account for fractional derivative elements in the system equation of motion. Alternatively, it can be viewed as a generalization of the formulation in Ref. [22] to yield enhanced computational efficiency.
Specifically, the developed technique constitutes an approximation-free dimension reduction approach that renders the associated computational cost independent of the total number of stochastic dimensions of the problem. In this regard, it can determine directly any lower-dimensional joint response PDF corresponding to a subset only of the response vector components. This is done by utilizing an appropriate combination of fixed and free boundary conditions in the related functional minimization problem. An indicative numerical example pertaining to a stochastically excited TMDI nonlinear system with a fractional derivative element has been considered for demonstrating the accuracy and efficiency of the technique compared to relevant MCS data.
Overall, it has been shown that the reduced-order WPI technique is ideally suited for problems where the objective relates to determining the response of few only specific DOF that are critical for the design and optimization of the overall system.
Acknowledgment
I. A. Kougioumtzoglou gratefully acknowledges the support through his CAREER award by the CMMI Division of the National Science Foundation (Award No. 1748537).
Funding Data
Directorate for Engineering (Award No. 1748537; Funder ID: 10.13039/100000084).
Data Availability Statement
The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.