## Abstract

We present the analytical solutions of the second-, third-, and fourth-order response moments of a single-degree-of-freedom linear system subjected to a class of non-Gaussian random excitation. The non-Gaussian excitation is a zero-mean stationary stochastic process prescribed by an arbitrary probability density and a power spectrum whose peak is located at zero frequency. The excitation is described by an Itô stochastic differential equation in which the drift and diffusion coefficients are determined from the probability density and spectral density of the excitation. In order to obtain the analytical solutions of the response moments, first, we derive the third- and fourth-order autocorrelation functions of the non-Gaussian excitation using its Markov and detailed balance properties. The third-order correlation function is given by the same expression regardless of the difference in the probability density function of the excitation. On the other hand, the fourth-order correlation function is derived under the assumption that the excitation probability density belongs to the Pearson distribution family, which is one of the widest classes of probability distributions. Then, combining the autocorrelation functions of the excitation and the convolution representation of the response, we obtain the exact solutions of the response moments, and it is shown what kind of components the response moments are composed of. Finally, we investigate the dominant time-varying components of the response moments for several different excitation bandwidths.

## 1 Introduction

It is important to determine and examine the stochastic response of dynamic systems subjected to random excitations from the viewpoint of enhancement of reliability, safety, and comfortability of machines and structures. When analyzing stochastic dynamic systems, in many cases, a random excitation has been modeled by a Gaussian process. This is due to the fact that real random excitations often have the probability densities similar to a Gaussian distribution, and the analytical advantage that the probabilistic properties of Gaussian processes can be completely described by the statistics up to the second order, namely, the mean function and the autocorrelation function (or equivalently the power spectral density for stationary processes). However, some engineering systems are subjected to highly non-Gaussian random excitations, for example, vehicles running on rough road surfaces [1,2], trains running on tracks with irregularities [1,3], marine structures excited by wave forces [4,5], and low-rise buildings under wind pressure [6,7]. The response of a system under non-Gaussian excitation is generally also non-Gaussian even though the system is linear. Thus, if we assume a Gaussianity for a random excitation that is actually non-Gaussian, the assumption may lead to large errors in response analysis results. Therefore, the stochastic response analysis taking into account the non-Gaussianity of random excitations properly and the elucidation of the effects of the excitation non-Gaussianity on system responses have become one of the key issues in the community of random vibration.

A mathematically rigorous description of a non-Gaussian process requires its finite-dimensional distribution of arbitrary order, which is unrealistic in modeling real non-Gaussian phenomena. For this reason, in many practical situations, two essential probabilistic properties, the first-order probability density and the power spectral density, are used. In the last few decades, considerable effort has been devoted to developing models for a non-Gaussian process prescribed by the probability density and the power spectrum. And now, there are various models such as translation process [8,9], spectral representation method with translation process theory [10–15], Karhunen–Loève expansion model [16–18], polynomial chaos expansion model [19–21], stochastic differential equation model [22–28], etc.

Among such non-Gaussian models, this paper focuses on a stochastic differential equation model presented by Cai and Lin [25]. In this model, a non-Gaussian process is assumed to be zero-mean and stationary, and described by a one-dimensional Itô stochastic differential equation. The drift and diffusion coefficients in the Itô equation are adjusted to match the spectral density and the probability density, respectively. The advantages of the Cai and Lin’s model are as follows: (i) this model is very versatile in terms that it can treat a non-Gaussian process with an arbitrary probability density, which is useful to find the general characteristics of system response that commonly appear for various non-Gaussian excitations. (ii) No approximations are involved in modeling, and a model that exactly fits the given non-Gaussian distribution and power spectrum can be obtained. (iii) Since the random process is represented as a Markov process, it is possible to perform an analysis using the well-known Markov property.

So far, some studies have been conducted for linear and nonlinear systems under non-Gaussian random excitation described by the Cai and Lin’s model. Tsuchida and Kimura carried out Monte Carlo simulation to examine the stationary response probability density of single-degree-of-freedom (SDOF) linear and nonlinear-stiffness systems under the non-Gaussian excitation [29]. It was found that the shape of the response distribution is highly dependent on the bandwidth of the excitation power spectrum. Especially for a linear system, the displacement response distribution has a shape close to the excitation probability density when the excitation bandwidth is small compared to the bandwidth of the frequency response function of the system. Conversely, if the excitation bandwidth is larger, then the response distribution becomes almost Gaussian. Such a relationship between the excitation bandwidth and the stationary response of the system was also investigated in terms of the moments [30–32].

In many engineering problems, it is important to understand not only the stationary response but also the transient response. A few studies have addressed the transient response of systems subjected to non-Gaussian excitation modeled by the Cai and Lin’s procedure. Wu and Cai investigated the effects of the excitation probability density on the transient response properties for various linear and nonlinear dynamical systems in the case of wide-band excitation [33]. It was shown that the excitation probability distribution has a great effect on the transient behavior of the response for a dynamical system, either linear or nonlinear, and the effect reduces as the system response approaches the stationary state. Recently, Fukushima and Tsuchida also studied the transient response of a SDOF linear system subjected to the non-Gaussian excitation with a variety of bandwidths using Monte Carlo simulation [34]. The simulation results demonstrated that as in the case of the stationary response [29], the bandwidth of the excitation spectral density influences on the transient response characteristics significantly. Immediately after the excitation starts, the effect of non-Gaussianity of the excitation appears in the response distribution clearly, but when the excitation bandwidth is large, the response distribution approaches a Gaussian distribution almost monotonically and rapidly. On the other hand, for the small excitation bandwidth, the effect of the excitation non-Gaussianity remains in the response until the response reaches the stationary state, and the width and degree of non-Gaussianity of the response distribution change periodically in the transient state. In the above studies, however, it was not sufficiently investigated what kind of time-varying components the probability distribution and moments of the transient response are composed of, and which of them are dominant. A better understanding of such time-varying characteristics of the transient response is expected to contribute to improving the reliability and safety of non-Gaussian randomly excited systems. For non-Gaussian response, its probability density function generally has a shape with asymmetry and/or heavy tails. In terms of statistical moments, the asymmetry and tail shape of the probability density are characterized by the third- and fourth-order moments, which are related to the skewness and kurtosis of the response, respectively. Thus, the non-Gaussian characteristics of the response can be found by analyzing these higher-order response moments in detail.

In the present paper, we derive the analytical solutions of the second-, third-, and fourth-order time-dependent response moments of a SDOF linear system subjected to non-Gaussian random excitation expressed by the Cai and Lin’s model. In this regard, the higher-order autocorrelation functions of the non-Gaussian excitation, which are necessary for obtaining the response statistics, are also derived using its Markov and detailed balance properties. The response moments are then decomposed into their components to see what time-varying components are included. Furthermore, we also observe the relationship between the dominant components of the response moments and the excitation bandwidth.

## 2 Models of System and Non-Gaussian Random Excitation

### 2.1 Equation of Motion.

*X*denotes the displacement,

*ζ*is the damping ratio, and

*U*(

*t*) is a non-Gaussian random excitation.

*t*is the time variable non-dimensionalized by the undamped natural angular frequency of the system. The dot symbol is the differentiation with respect to

*t*. We assume that the system is initially at rest and the excitation

*U*(

*t*) is added to the system at

*t*= 0.

### 2.2 Non-Gaussian Random Excitation.

*U*(

*t*) considered in this paper is a zero-mean stationary stochastic process described by a model developed by Cai and Lin [25]. In this model,

*U*(

*t*) is prescribed by the probability density function

*p*

_{U}(

*u*) and the following power spectrum

*S*

_{U}(

*ω*):

*α*> 0 is the bandwidth parameter and E[

*U*

^{2}] is the mean square value of

*U*(

*t*).

*S*

_{U}(

*ω*) for five types of

*α*and E[

*U*

^{2}] = 1 is shown in Fig. 1. This spectral density has one peak at

*ω*= 0. In the case of small

*α*, the zero frequency component is dominant; on the other hand,

*S*

_{U}(

*ω*) is wide-band for large

*α*. Although

*α*is not often referred to as bandwidth when the peak of the spectrum is located at

*ω*= 0, in this paper, we refer to

*α*as bandwidth for the sake of convenience.

*U*(

*t*) with an arbitrary probability density

*p*

_{U}(

*u*) and the power spectrum

*S*

_{U}(

*ω*) given by Eq. (2) is expressed by the following one-dimensional Itô stochastic differential equation:

*α*is the same parameter in Eq. (2),

*B*(

*t*) is a standard Brownian motion, and −

*αu*and

*D*

^{2}(

*u*) are the drift and diffusion coefficients, respectively. Multiplying Eq. (3) by

*U*(

*t*−

*τ*) and taking the ensemble average, we obtain

*R*

_{U}(

*τ*) = E[

*U*(

*t*)

*U*(

*t*−

*τ*)] is the autocorrelation function of

*U*(

*t*), and in deriving this equation, the independent increments property of the Brownian motion

*B*(

*t*) was used. Since

*R*

_{U}(0) corresponds to E[

*U*

^{2}], the solution of Eq. (4) is given by

*U*(

*t*) is expressed by Eq. (2).

*p*

_{U}(

*u*) of

*U*(

*t*) is governed by the reduced Fokker–Planck equation corresponding to Eq. (3)

*G*(

*u*) is known as the probability flow. In the present one-dimensional case,

*G*(

*u*) must vanish everywhere [35], and Eq. (6) reduces to

*D*

^{2}(

*u*) in Eq. (3) is determined according to the probability density function

*p*

_{U}(

*u*). It is also noted that

*U*(

*t*) belongs to the class of detailed balance because Eq. (7) holds [35]. This fact will be used to derive the higher-order autocorrelation functions of

*U*(

*t*) in the next section.

Thus, the non-Gaussian random excitation *U*(*t*) described by Eq. (3) with *D*(*u*) in Eq. (8) possesses a given probability density *p*_{U}(*u*) and the spectral density *S*_{U}(*ω*) given by Eq. (2). For this non-Gaussian excitation, we analyze the second-, third-, and fourth-order response moments of the SDOF linear system (1).

## 3 Derivation of Higher-Order Autocorrelation Functions of Non-Gaussian Excitation

In this section, the third- and fourth-order autocorrelation functions of the non-Gaussian random excitation *U*(*t*) are derived.

### 3.1 Third-Order Autocorrelation Function.

*t*

_{1},

*t*

_{2}, and

*t*

_{3}, let

*t*

_{1}≤

*t*

_{2}≤

*t*

_{3}. To simplify the notation, the excitation at each time is denoted as

*U*(

*t*

_{1}) =

*U*

_{1},

*U*(

*t*

_{2}) =

*U*

_{2}, and

*U*(

*t*

_{3}) =

*U*

_{3}. Since the random excitation

*U*(

*t*) is described by the stochastic differential equation (3),

*U*(

*t*) is a diffusive Markov process. Therefore, the third-order joint probability density $pU1U2U3(u1,u2,u3)$ of

*U*(

*t*) can be written as

*U*

_{i}under the condition that

*U*

_{j}=

*u*

_{j}. By using Eq. (9), the third-order autocorrelation function E[

*U*

_{1}

*U*

_{2}

*U*

_{3}] of

*U*(

*t*) can be calculated as follows:

*U*

_{2}=

*u*

_{2}] is the conditional expectation conditioned on

*U*

_{2}=

*u*

_{2}at time

*t*

_{2}.

*U*

_{3}|

*U*

_{2}=

*u*

_{2}] and E[

*U*

_{1}|

*U*

_{2}=

*u*

_{2}] in Eq. (10). On both sides of Eq. (3), taking the conditional average E[·|

*U*

_{2}=

*u*

_{2}] under

*t*≥

*t*

_{2}and dividing by

*dt*lead to

*t*

_{3}≥

*t*

_{2}, we obtain

*U*

_{1}|

*U*

_{2}=

*u*

_{2}] (

*t*

_{1}≤

*t*

_{2}), Eq. (12) is unable to use directly. Accordingly, we utilize the fact that because the excitation

*U*(

*t*) satisfies the detailed balance as mentioned in Sec. 2.2, the following equation for its joint probability density function $pU1U2(u1,u2)$ holds [36].

*t*

_{2}≥

*t*

_{1}, the conditional expectation E[

*U*

_{2}|

*U*

_{1}=

*u*

_{2}] in the right-hand side of Eq. (15) can be obtained by the same procedure as in Eqs. (11)–(13)

*U*

_{1}

*U*

_{2}

*U*

_{3}] of the excitation

*U*(

*t*) is derived as

*U*(

*t*). Although Eq. (18) is valid only for

*t*

_{1}≤

*t*

_{2}≤

*t*

_{3}, the third-order correlation function with different order of

*t*

_{1},

*t*

_{2}, and

*t*

_{3}can be obtained by replacing the three indices of time

*t*in Eq. (18). For example, E[

*U*

_{1}

*U*

_{2}

*U*

_{3}] with

*t*

_{3}≤

*t*

_{1}≤

*t*

_{2}is given by

*p*

_{U}(

*u*) of the excitation

*U*(

*t*).

### 3.2 Fourth-Order Autocorrelation Function.

*t*

_{1},

*t*

_{2},

*t*

_{3}, and

*t*

_{4}, let

*t*

_{1}≤

*t*

_{2}≤

*t*

_{3}≤

*t*

_{4}. The excitation at each time is denoted as

*U*(

*t*

_{1}) =

*U*

_{1},

*U*(

*t*

_{2}) =

*U*

_{2},

*U*(

*t*

_{3}) =

*U*

_{3}, and

*U*(

*t*

_{4}) =

*U*

_{4}. As in Sec. 3.1, using the Markov property of the excitation

*U*(

*t*), the fourth-order joint probability density $pU1U2U3U4(u1,u2,u3,u4)$ of

*U*(

*t*) can be expressed as

*U*

_{1}

*U*

_{2}

*U*

_{3}

*U*

_{4}] of

*U*(

*t*) is calculated as follows:

*U*

_{1}|

*U*

_{2}=

*u*

_{2}] is given by Eq. (17), and E[

*U*

_{4}|

*U*

_{3}=

*u*

_{3}] can also be obtained through the same procedure as in Eqs. (11)–(13)

*U*

^{2}(

*t*)

*U*

_{2}=

*u*

_{2}] under

*t*≥

*t*

_{2}and dividing by

*dt*yield

*D*

^{2}(

*u*) is of the general form. However, if

*D*

^{2}(

*u*) is expressed by a quadratic polynomial in

*u*, we can find the solution of Eq. (25). Therefore, hereafter, consider the case where

*D*

^{2}(

*u*) is written in quadratic form as

*D*

^{2}(

*u*) =

*Au*

^{2}+

*Bu*+

*C*. This type of diffusion coefficient corresponds to the case in which the excitation probability distribution

*p*

_{U}(

*u*) belongs to the Pearson distribution family [25]. Substituting

*D*

^{2}(

*u*) =

*Au*

^{2}+

*Bu*+

*C*into Eq. (25), we have

*U*(

*t*)|

*U*

_{2}=

*u*

_{2}] is given by Eq. (12). The solution of Eq. (26) when

*A*≠

*α*, 2

*α*is

*U*

_{1}

*U*

_{2}

*U*

_{3}

*U*

_{4}] of

*U*(

*t*) as follows:

*U*

^{4}] are the second-, third-, and fourth-order moments of

*U*(

*t*), respectively. Similar to the third-order case in Sec. 3.1, E[

*U*

_{1}

*U*

_{2}

*U*

_{3}

*U*

_{4}] with different orders of

*t*

_{1},

*t*

_{2},

*t*

_{3}, and

*t*

_{4}can be obtained by replacing the four indices of time

*t*in Eq. (29). Incidentally, it is possible to get the autocorrelation function in the case of

*A*=

*α*and

*A*= 2

*α*by taking the limit of

*A*→

*α*and

*A*→ 2

*α*in Eq. (29), respectively. The solutions are provided in Appendix A.

## 4 Analysis of Transient Response Moments

*X*(

*t*) of the system (1) is written as a convolution integral of the excitation

*U*(

*t*) and the impulse response function

*h*(

*t*) of the system.

*h*(

*t*) is given by

*ω*

_{d}is the damped natural angular frequency $\omega d=1\u2212\zeta 2$. The velocity response $X\u02d9(t)$ is also described as

*n*th-order response moments E[

*X*

^{n}(

*t*)] and $E[X\u02d9n(t)]$ are

*U*(

*t*

_{1})

*U*(

*t*

_{2}) · · ·

*U*(

*t*

_{n})] is the

*n*th-order autocorrelation function of the excitation

*U*(

*t*).

In the following sections, we find the analytical solutions of the second-, third-, and fourth-order response moments by using the autocorrelation functions of *U*(*t*) derived in Sec. 3. Furthermore, based on the analytical solutions, we investigate what kind of time-varying components each moment consists of.

### 4.1 Second-Order Moments.

*X*

^{2}(

*t*)] and the second-order velocity moment $E[X\u02d92(t)]$ are expressed, respectively, as follows:

*U*(

*t*

_{1})

*U*(

*t*

_{2})] is the autocorrelation function of

*U*(

*t*), and by Eq. (5), E[

*U*(

*t*

_{1})

*U*(

*t*

_{2})] = E[

*U*

^{2}]exp ( −

*α*|

*t*

_{2}−

*t*

_{1}|). Substituting it into Eqs. (35) and (36), and performing integration, we can obtain the analytical solutions of E[

*X*

^{2}(

*t*)] and $E[X\u02d92(t)]$. The integral calculation was done using Matlab (R2021a) in this study. Since the analytical solutions are given by very long expressions, they are included in Appendix B. These solutions are not dependent on the shape of the probability density

*p*

_{U}(

*u*) of the excitation

*U*(

*t*) although they are proportional to E[

*U*

^{2}]. From the analytical solutions, it is found that both E[

*X*

^{2}(

*t*)] and $E[X\u02d92(t)]$ consist of the following terms (we assume

*ω*

_{d}≃ 1, i.e., the system is lightly damped):

- Exponentially decaying oscillatory terms with period 2
*π*$e\u2212(\zeta +\alpha )tsin(\omega dt),e\u2212(\zeta +\alpha )tcos(\omega dt)$ - Exponentially decaying oscillatory terms with period
*π*$e\u22122\zeta tsin(2\omega dt),e\u22122\zeta tcos(2\omega dt)$ Monotonically decaying term

$e\u22122\zeta t$Stationary term (constant term)

*t*→ ∞ in the analytical solution, which corresponds to the stationary moment.

### 4.2 Third-Order Moments.

*X*

^{3}(

*t*)] and the third-order velocity moment $E[X\u02d93(t)]$ are given, respectively, as follows:

*U*(

*t*

_{1})

*U*(

*t*

_{2})

*U*(

*t*

_{3})] derived in Sec. 3.1 into Eqs. (37) and (38) and then performing integration. When integrating, we need to take into account that Eq. (18) is valid only for

*t*

_{1}≤

*t*

_{2}≤

*t*

_{3}, and that there are a total of six combinations of the order of

*t*

_{1},

*t*

_{2}, and

*t*

_{3}. As mentioned at the end of Sec. 3.1, E[

*U*(

*t*

_{1})

*U*(

*t*

_{2})

*U*(

*t*

_{3})] with different order of

*t*

_{1},

*t*

_{2}, and

*t*

_{3}can be obtained by replacing the three indices of time

*t*in Eq. (18). With the above in mind, we divide the integral into six parts

*X*

^{3}(

*t*)] and $E[X\u02d93(t)]$ are proportional to the third-order excitation moment E[

*U*

^{3}]. The integrals in Eqs. (40) and (41) were calculated using Matlab to obtain the analytical solutions of the third-order response moments. These solutions are provided in Appendix B. We observe from the analytical solutions that both E[

*X*

^{3}(

*t*)] and $E[X\u02d93(t)]$ consist of the following terms:

- Exponentially decaying oscillatory terms with period 2
*π*$e\u2212(\zeta +\alpha )tsin(\omega dt),e\u2212(\zeta +\alpha )tcos(\omega dt)e\u22123\zeta tsin(\omega dt),e\u22123\zeta tcos(\omega dt)$ - Exponentially decaying oscillatory terms with period
*π*$e\u2212(2\zeta +\alpha )tsin(2\omega dt),e\u2212(2\zeta +\alpha )tcos(2\omega dt)$ - Exponentially decaying oscillatory terms with period 2
*π*/3$e\u22123\zeta tsin(3\omega dt),e\u22123\zeta tcos(3\omega dt)$ - Monotonically decaying term$e\u2212(2\zeta +\alpha )t$
Stationary term (constant term)

*t*→ ∞ in the analytical solution (i.e., the stationary third-order moment).

### 4.3 Fourth-Order Moments.

*X*

^{4}(

*t*)] and the fourth-order velocity moment $E[X\u02d94(t)]$ are written, respectively, as

*t*

_{1}≤

*t*

_{2}≤

*t*

_{3}≤

*t*

_{4}in which Eq. (29) is applicable. Then, multiplying the integral result by 24, which is the total number of permutations of

*t*

_{1},

*t*

_{2},

*t*

_{3}, and

*t*

_{4}, leads to E[

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$, namely,

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$. These solutions are much longer than those of the second- and third-order moments, and thus, the solutions of the fourth-order moments in their complete forms are not shown in this paper. However, the terms contained in E[

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$ are common and indicated below.

- Exponentially decaying oscillatory terms with period 2
*π*$e\u2212(\zeta +\alpha )tsin(\omega dt),e\u2212(\zeta +\alpha )tcos(\omega dt),e\u2212(3\zeta +\alpha )tsin(\omega dt),e\u2212(3\zeta +\alpha )tcos(\omega dt)$ - Exponentially decaying oscillatory terms with period
*π*$e\u22122\zeta tsin(2\omega dt),e\u22122\zeta tcos(2\omega dt),e\u2212(2\zeta +\alpha )tsin(2\omega dt),e\u2212(2\zeta +\alpha )tcos(2\omega dt),e\u22124\zeta tsin(2\omega dt),e\u22124\zeta tcos(2\omega dt),e\u2212(2\zeta +2\alpha \u2212A)tsin(2\omega dt),e\u2212(2\zeta +2\alpha \u2212A)tcos(2\omega dt)$ - Exponentially decaying oscillatory terms with period 2
*π*/3$e\u2212(3\zeta +\alpha )tsin(3\omega dt),e\u2212(3\zeta +\alpha )tcos(3\omega dt)$ - Exponentially decaying oscillatory terms with period
*π*/2$e\u22124\zeta tsin(4\omega dt),e\u22124\zeta tcos(4\omega dt)$ - Monotonically decaying terms$e\u22122\zeta t,e\u2212(2\zeta +\alpha )t,e\u22124\zeta t,e\u2212(2\zeta +2\alpha \u2212A)t$
Stationary term (constant term)

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$ depend upon not only the excitation moment of the same order E[

*U*

^{4}] but also the second-order moment E[

*U*

^{2}] and the third-order moment E[

*U*

^{3}]. This feature is due to the fact that the fourth-order autocorrelation function of the excitation (Eq. (29)) includes E[

*U*

^{2}] and E[

*U*

^{3}] as well as E[

*U*

^{4}].

## 5 Analysis Results

In this section, the magnitude of each component in the analytical solution of the response moment is compared quantitatively to examine its dominant time-varying component. The values of the damping ratio *ζ* of the system and the bandwidth *α* of the excitation power spectrum *S*_{U}(*ω*) are given as follows:

Damping ratio:

*ζ*= 0.05Bandwidth of excitation power spectrum:

*α*= 0.01, 0.05, 1

*ζ*corresponds to the bandwidth of its frequency response function.

*α*is also the non-dimensional bandwidth parameter of

*S*

_{U}(

*ω*). Hence,

*ζ*=

*α*indicates that the bandwidth of the frequency response function of the system and that of

*S*

_{U}(

*ω*) are equal. In the previous papers [29,34], we investigated the stationary and transient responses of the same linear system as in this study, and it was revealed that the response characteristics change significantly according to whether

*α*is larger or smaller than

*ζ*. Based on this knowledge, we consider three types of

*α*, which are smaller, equal, and larger values compared to

*ζ*, respectively.

*a*> 0 and

*b*> 0 are the shape and scale parameters, respectively. The diffusion coefficient

*D*

^{2}(

*u*) corresponding to the shifted gamma distribution (46) is

*a*and

*b*were chosen as follows:

*a*=

*b*= 1 is shown in Fig. 2. In the Monte Carlo simulation procedure, first, the sample functions of the excitation

*U*(

*t*) are generated by solving Eq. (3) with (47) numerically using the Euler–Maruyama scheme [37]. Then, by the use of the generated sample functions of

*U*(

*t*), the linear system (1) is numerically integrated by the fourth-order Runge–Kutta algorithm to obtain the sample functions of the displacement

*X*(

*t*) and the velocity $X\u02d9(t)$, which are employed to evaluate the second-, third-, and fourth-order response moments. The simulation conditions are as follows:

Time-step: Δ

*t*= 0.1 (*α*= 0.01 and 0.05), Δ*t*= 0.01 (*α*= 1)Number of sample functions: 2 × 10

^{7}

*X*

^{2}(

*t*)], $E[X\u02d92(t)]$, E[

*X*

^{3}(

*t*)], and $E[X\u02d93(t)]$, and to excitation with the Pearson distribution family for E[

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$.

### 5.1 Second-Order Moments.

Figure 3 shows the temporal change of each component of the second-order response moments E[*X*^{2}(*t*)] and $E[X\u02d92(t)]$ derived in Sec. 4.1. Since E[*X*^{2}(*t*)] and $E[X\u02d92(t)]$ are proportional to the second-order moment E[*U*^{2}] of the excitation *U*(*t*), in Fig. 3, the results of E[*X*^{2}(*t*)] and $E[X\u02d92(t)]$ are normalized by E[*U*^{2}]. If there are multiple terms in each component, the sum of them is plotted. The black solid line (total) denotes the sum of all components, which is namely E[*X*^{2}(*t*)]/E[*U*^{2}] itself or $E[X\u02d92(t)]/E[U2]$ itself. The circles indicate the corresponding Monte Carlo simulation results for comparison. The analytical results (total) are in good agreement with the simulation results, which prove that the analytical solutions derived in this study are accurate.

The characteristics of the time variation of E[*X*^{2}(*t*)] and $E[X\u02d92(t)]$, which can be seen from Fig. 3, are described below. First, we observe the results in the case that the excitation bandwidth parameter *α* is smaller than the damping ratio *ζ* which corresponds to the bandwidth of the frequency response function. Comparing the components of E[*X*^{2}(*t*)] for *α* = 0.01, the exponentially decaying oscillation with period 2*π* is dominant. Correspondingly, the period of E[*X*^{2}(*t*)] (total) is also 2*π*. On the other hand, for $E[X\u02d92(t)]$, its *π*-period component is superior to other components, so that the period of $E[X\u02d92(t)]$ is also *π*. For *α* = 0.05 (*α* = *ζ*), as in the case of *α* = 0.01, the 2*π*- and *π*-period components dominate in E[*X*^{2}(*t*)] and $E[X\u02d92(t)]$, respectively. It can also be seen that there is no monotonically decaying component in either E[*X*^{2}(*t*)] or $E[X\u02d92(t)]$ (“decay” line is always 0). For *α* = 1 (*α* > *ζ*), in both E[*X*^{2}(*t*)] and $E[X\u02d92(t)]$, the monotonic decay component is dominant. Comparing between the oscillatory components, the *π*-period component is larger, although its magnitude itself is small. Therefore, the second-order response moments for *α* = 1 increase with a small oscillation with period *π*.

### 5.2 Third-Order Moments.

Figure 4 shows each component of the third-order response moments E[*X*^{3}(*t*)] and $E[X\u02d93(t)]$ derived in Sec. 4.2. Similar to the second-order moments, E[*X*^{3}(*t*)] and $E[X\u02d93(t)]$ are proportional to the third-order excitation moment E[*U*^{3}]. Hence, in Fig. 4, E[*X*^{3}(*t*)]/E[*U*^{3}] and $E[X\u02d93(t)]/E[U3]$ are plotted assuming E[*U*^{3}] ≠ 0. (When E[*U*^{3}] = 0, E[*X*^{3}(*t*)] and $E[X\u02d93(t)]$ are equal to 0). The analytical results (total) agree very well with the simulation results.

From Fig. 4, the oscillation component with period 2*π* is dominant for both E[*X*^{3}(*t*)] and $E[X\u02d93(t)]$ irrespective of the excitation bandwidth *α*. Since the other oscillatory and monotonic decay components are all smaller than the 2*π*-period component, these terms do not contribute much to the time variation of E[*X*^{3}(*t*)] and $E[X\u02d93(t)]$.

### 5.3 Fourth-Order Moments.

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$ are not proportional to the fourth-order excitation moment E[

*U*

^{4}]. In addition, E[

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$ are dependent on not only E[

*U*

^{4}] but also the second-order moment E[

*U*

^{2}] and the third-order moment E[

*U*

^{3}]. Therefore, for the fourth-order response moments, it is not possible to show general results that apply to any excitation distribution as those in Figs. 3 and 4. For this reason, we show the results for two types of the excitation probability densities, which are a gamma distribution and a uniform distribution. The gamma distribution considered here is given by Eqs. (46) and (48), and the corresponding diffusion coefficient is expressed by Eq. (47). The uniform distribution, on the other hand, is written as

*a*=

*b*= 1 and $\Delta =3$ are shown below.

Gamma distribution E[

*U*^{2}] = 1 E[*U*^{3}] = 2 E[*U*^{4}] = 9Uniform distribution E[

*U*^{2}] = 1 E[*U*^{3}] = 0 E[*U*^{4}] = 1.8

*D*

^{2}(

*u*) corresponding to the gamma and uniform distributions are linear and quadratic, respectively (see Eqs. (47) and (50)), the fourth-order autocorrelation function of random excitations with these non-Gaussian distributions is expressed by Eq. (29) (

*A*= 0,

*B*= 2

*αb*,

*C*= 2

*αab*

^{2}for gamma distribution and

*A*= −

*α*,

*B*= 0,

*C*=

*α*Δ

^{2}for uniform distribution). Thus, the analytical solutions of E[

*X*

^{4}(

*t*)] and $E[X\u02d94(t)]$ can be obtained by the manner described inSec. 4.3.

Figures 5 and 6 show each component of the fourth-order response moments E[*X*^{4}(*t*)] and $E[X\u02d94(t)]$ for the gamma and uniform distributions, respectively. Note that Figs. 5 and 6 show the fourth-order response moments themselves rather than the normalized ones. The simulation results in Fig. 6 were obtained with the pertinent Monte Carlo simulation using 2 × 10^{7} realizations of the uniformly distributed excitation. These two figures demonstrate that the analytical solutions derived in this study are precise.

From Figs. 5 and 6, it is found that the temporal change of the fourth-order response moments has some common features regardless of the excitation probability distribution. For *α* = 0.01 and 0.05, the 2*π*- and *π*-period components are dominant in E[*X*^{4}(*t*)] and $E[X\u02d94(t)]$, respectively. In the case of *α* = 1, the monotonic decay component is dominant. Among the oscillatory components, the *π*-period component is slightly larger, and the remaining three oscillatory components (period 2*π*, 2*π*/3, and *π*/2) are almost zero. Therefore, the fourth-order response moments for *α* = 1 increase with a small oscillation with period *π*. Finally, the above features are similar to those observed in the second-order response moments. Hence, the characteristics of the temporal change in the response moments are common for the second- and fourth-order moments and different for the third-order moment.

## 6 Conclusions

We have derived the analytical solutions of the response statistics of a single-degree-of-freedom linear system subjected to a class of non-Gaussian random excitation. The non-Gaussian excitation considered in this paper is a zero-mean stationary stochastic process prescribed by an arbitrary probability density and a power spectrum with bandwidth parameter. The excitation is modeled by a one-dimensional Itô-type stochastic differential equation whose drift and diffusion coefficients are determined according to the probability density and the spectral density.

In order to obtain the analytical solutions of the response statistics, first, the third- and fourth-order autocorrelation functions of the non-Gaussian excitation have been derived using its Markov and detailed balance properties. The third-order correlation function is valid for any probability density function of the excitation. The fourth-order correlation function, on the other hand, has been derived under the assumption that the diffusion coefficient of the stochastic differential equation for the excitation is expressed as a quadratic polynomial. However, the quadratic-type diffusion coefficient corresponds to the fact that the probability density belongs to the Pearson distribution family, which is one of the widest classes of probability distributions. Hence, the fourth-order correlation function derived in this study is applicable to excitation with a wide variety of probability density functions. Then, using the autocorrelation functions of the excitation, the second-, third-, and fourth-order response moments have been found based on the convolution integral of the excitation and the impulse response function of the system. It has been shown that the response moments consist of the several exponentially decaying oscillatory components with different periods, the monotonically decaying components and the stationary component.

Finally, we have investigated the dominant time-varying components of the response moments in the case of three different excitation bandwidths. The findings are as follows:

*Second- and fourth-order moments*

When the excitation bandwidth is small or equal to the system damping ratio corresponding to the bandwidth of the system, the displacement moment is dominated by the decaying oscillation component with period 2

*π*. On the other hand, the decaying oscillation with period*π*predominates in the velocity moment.When the excitation bandwidth is large compared to the damping ratio, the monotonically decaying component becomes the most dominant in both the displacement and velocity moments. Among the oscillation components, the

*π*-period component is slightly larger, and other components are nearly zero. Therefore, the moments increase with a small oscillation with period*π*.

*Third-order moment*

Irrespective of the excitation bandwidth, the decaying oscillation component with period 2

*π*dominates for both the displacement and velocity moments.

In this paper, a stationary non-Gaussian process has been considered as a random excitation acting on a system. In some engineering problems, random excitations may have not only non-Gaussianity but also non-stationarity. Therefore, it is of importance to analyze dynamic systems under non-stationary non-Gaussian excitations. The analytical solutions derived in this paper can be applied to the case of non-stationary non-Gaussian excitation if the excitation is given in the form of the stationary non-Gaussian process described by Cai and Lin’s model multiplied by a modulating function introducing amplitude non-stationarity. This topic will be examined concretely in our future work and reported in a later paper.

## Acknowledgment

This work was supported by JSPS KAKENHI Grant Number JP18K13712. The figures were drawn using ZoomPlot (Matlab extension) provided by Qiu, K.^{2}

## Conflict of Interest

The authors have no competing interests to declare that are relevant to the content of this paper.

## Data Availability Statement

The authors attest that all data for this study are included in the paper.

### Appendix A: The Fourth-Order Autocorrelation Function of Excitation When *A* = *α* and 2*α*

*Case of A =*

*α**Case of A = 2*

*α**t*

_{1}≤

*t*

_{2}≤

*t*

_{3}≤

*t*

_{4}.

### Appendix B: Analytical Solutions of the Second- and Third-Order Response Moments

*The second-order displacement moment*

*The second-order velocity moment*

*The third-order displacement moment*

*The third-order velocity moment*