1 Introduction

Uncertainty quantification is important in many scientific and engineering applications, e.g., fluid dynamics^[1-2], engineering mechanics^[3-5], and material science^[6-7]. During the past decades, the rapid developments in computational and numerical methods have enabled system characteristics to be quantified under uncertain conditions. Many methods have been developed to quantify the uncertainties, e.g., the generalized polynomial chaos (GPC) expansion approach^[8-11], the stochastic collocation method^[12-14], the fuzzy set theory^[15-16], the interval analysis method^[17-19], and the random matrix theory^[20-21].

Among the above methods, the GPC method^[22] has attracted increasing attention^[23-26], and is being used widely in many fields^[27-32]. In the GPC method, the stochastic response is represented by the GPC bases, which are time invariant and determined by the probability density functions of uncertain variables. Owing to the time-invariant bases, the GPC method can be simply implemented, and it is efficient in representing uncertainty in some complex models. However, according to Refs. [25] and [33], the GPC method will cause large errors in long-time simulations for some differential equations, e.g., the Kraichnan-Orszag three-mode ordinary-differential-equation (ODE) system^[34] and the Navier-Stokes equations for unsteady noisy flows^[28]. Therefore, many methods are proposed to overcome this computational bottleneck, e.g., the time-dependent GPC method (TD-GPC)^[25]. In the TD-GPC method, when the nonlinear terms of the GPC are relatively significant compared with the linear coefficients, the time integration of the GPC procedure will be terminated, and a new random variable will be introduced to represent the current stochastic response. Through this process, new orthogonal polynomials are constructed with respect to the new random variable. Subsequently, new ODEs can be constructed to propagate the uncertainty through the system. The TD-GPC method is effective for problems with a small number of degrees of freedom. However, it will encounter a convergence breakdown when uncertain parameters are absent in the initial conditions^[35]. In addition, the method is not suitable for problems with a continuous spatial field, e.g., heat transfer or computational fluid problems.

In this study, a new method is developed to modify the stochastic bases adaptively to yield a good approximation for the stochastic response surface as time progress. In the proposed method, the stochastic bases are assumed to be time dependent and orthonormal. The orthonormal property renders it convenient to evaluate the expansion coefficients. To maintain orthornormality, the derivatives of these bases are assumed to be orthogonal with the current bases. The orthonormality is not necessary, but it simplifies the evaluation of the dynamical stochastic bases. Based on the orthonormal property, the equations of time-dependent stochastic bases can be derived with the variational principle. Since the proposed method uses time-dependent bases (TDBs) instead of time-invariant bases, we will call it TDB for short.

The following sections are organized as follows. In Section 2, the problem statement presenting a brief introduction of the GPC method is described, the problem of long-term simulation error regarding the GPC method is discussed, and the time-dependent stochastic bases are introduced. In Section 3, the time-dependent stochastic basis-based method is provided, where the stochastic bases are assumed to be orthonormal, and the differential equations of TDBs are derived with the variational principle. In Section 4, some numerical examples are presented, e.g., a linear autonomous dynamic system and a clamped-clamped Euler-Bernoulli beam, to illustrate the efficiency of the proposed method. Finally, in Section 5, some conclusions are presented.

2 Problem statement

In general, a stochastic dynamical system can be expressed as follows:

(1)

where is the differential operator on the physical domain D, and ξ=(ξ₁, ξ₂, …, ξ_m)^T∈Ω represents the vector of m-dimensional uncertain parameters. The uncertainty may arise from the differential operator or stochastic forces.

In the GPC method, the stochastic response is approximated by the GPC bases as follows:

(2)

where H_ip^(p)(ξ_p) is the GPC basis with respect to the random variable ξ_p. These bases are time invariant and determined by the probability density functions of uncertain parameters. For example, we should use the Hermite polynomial chaos for Gaussian variables and the Legendre polynomial chaos for uniformly distributed random variables. To simplify the description of the proposed method, we introduce the multi-index I=(i₁, i₂, …, i_m) so that

(3)

Then, u(x, t, ξ) can be written as follows:

(4)

The GPC method is an important method in uncertainty quantification, and has been successfully used to solve a wide range of problems. However, there are some problems for the effective applications. When the problems are highly nonlinear or discontinuous, the GPC method will lead to slow convergence or convergence failure. Another well-known problem is that, the errors in the long-time simulation of the dynamical system will increase significantly over time. In Ref. [25], the long-term behavior of the GPC method was analyzed from the aspect of statistic properties. It is found that the GPC method cannot approximate the second-order statistics of the dynamic response. In this study, we will reconsider this problem from the aspect of stochastic response surface.

To obtain a better understanding of why the GPC method fails for some stochastic dynamical system, we consider a simple problem as follows:

(5)

where ξ₁ and ξ₂ are independent random variables distributing uniformly between -1 and 1, i.e., ξ₁, ξ₂∈ U[-1, 1].

The exact solution to the above initial value problem is

(6)

Figure 1 shows the solution of the GPC method compared with the exact solution. According to Fig. 1(b), the variance of the solution diverges from the exact solution after 7 s, and errors develop quickly with time.

To obtain a better understanding of how the stochastic response surfaces develop over time, some snapshots of the responses are plotted in Fig. 2. In the initial time, the stochastic responses are smooth functions of uncertain parameters. They can be represented well by low-order orthogonal polynomials. However, a peak develops around (ξ₁, ξ₁)=(1, 1) over time, thereby complicating the approximation of the stochastic response of the surface. In fact, polynomials cannot effectively represent uncertain responses, e.g.,

with a large t. The error of the response surface approximated by the GPC method is plotted in Fig. 3. It is obvious that the errors increase quickly over time.

In the GPC method, increasing the expansion order is a commonly used strategy to improve the accuracy of the results. However, a high-order expansion will lead to computational difficulty since the computational complexity increases exponentially with the order. In this problem, the GPC method fails because the bases are only determined by the uncertain parameters. These bases contain no information regarding the solution. To overcome this problem, we propose to use TDBs as an alternative. These bases are modified adaptively to yield a good approximation of the stochastic response surface over time.

3 Method based on time-dependent orthogonal bases

In this section, time-dependent stochastic bases are used to overcome some drawbacks of the GPC method. This method achieves a better approximation for the stochastic response surface by modifying the stochastic bases adaptively. These stochastic bases are modified according to the uncertain differential equation (see Eq. (1)) , where contains the derivatives of TDBs.

Using time-dependent orthogonal bases, the stochastic response can be expressed by

(7)

where c_i1i2…im (x, t) is the coefficient field, and ψ_ip^(p) (ξ_p, t) is a TDB of the random variable ξ_p. In the GPC method, the orthogonal polynomials are used to approximate the stochastic solution. The orthogonal property renders it simple to be implemented using the stochastic Galerkin method. In addition, the statistical properties of the solution, e.g., the mean and variance, can be evaluated directly from the result with little computational effort. Then, it is reasonable to maintain the orthornormality of the TDB ψ_ip^(p) (ξ_p, t) during numerical simulations as follows:

(8)

where E(·) is the expectation operator. It is noteworthy that the orthonormal property is not necessary. However, it is beneficial to simplify the computational procedure of the method.

The basis vector of the κth uncertain parameter is denoted by

(9)

Then, the orthonormal property in Eq. (8) is equal to

(10)

Differentiating this equation with respect to time yields

(11)

This implies that is an anti-symmetric matrix. To simplify the method, we strengthen this constraint to be

(12)

implying that the derivatives of bases are orthogonal to the current bases. In fact, the off-diagonal part of describes the transformation between the bases. If we consider the subspace S_κ spanned by the basis vector Ψ_κ (t), S_κ will remain invariant under the transformation between the bases. It is reasonable to ignore the off-diagonal part of ^[36].

Using TDBs, (x, t, ξ) can be expressed as

(13)

where

(14)

and T^(r)(t) is the vector of τ_ir^r(t) (1≤ i_r≤ N_r), i.e.,

(15)

Substituting Eq. (13) into the original differential equation (1) yields

(16)

Based on the variational principle, the evolution equations of the dynamical bases in Eq. (16) can be evaluated by (see Appendix A for the details regarding derivation)

(17)

In this equation, some notations are introduced to simplify the expression. The partial expectation operator E_i (·) with respect to a single random variable ξ_i is denoted by

(18)

where P_ξk (·) is an operator to project the stochastic function to a certain stochastic domain ξ_i. It is defined by

(19)

Hence, it is obvious that

(20)

In addition, the spatial integration operator, denoted by 〈·〉, can be expressed as follows:

(21)

The composition operators of P_ξk (·) and 〈·〉, denoted by 〈·〉_ξk, can be expressed by

(22)

Remark 1   Several methods can be used to determine the initial bases and propagate them over time. Typically, the initial conditions are simple functions of the uncertain parameters (such as the examples herein). In these cases, it is recommended to start with the GPC bases. In the GPC method, the bases are chosen to be a set of orthogonal polynomials.

Remark 2   When the initial condition of the problem cannot be represented by the chosen order of the polynomial chaos, we can approximate the stochastic field by using tensor decomposition methods such as higher-order singular value decomposition (HOSVD)^[37-38], which is an extension of the SVD method of high-order tensors and provides near-optimal bases. The numerical implementation of the time-dependent stochastic basis method requires a proper strategy to represent the TDB ψ_i^(κ) (ξ_k, t). Herein, the GPC bases are used. The TDB ψ_i^(κ) (ξ_k, t) can be approximated by the truncated polynomial chaos as follows:

(23)

where {H_i^(k) (ξ_k)}_i=1^∞ are the orthogonal polynomials with respect to the distribution function ρ_k (ξ_k) of the random variable ξ_k, and

Remark 3   In Eq. (17), 〈T^(k)(t)T^(k)T(t)〉_ξk may be ill-conditioned, and some regularized methods are required to calculate its inverse. We denote

It is obvious that A_ξk is a symmetrical matrix. The following regularized method is recommended:

(24)

where is a small number.

4 Numerical simulations 4.1 A linear autonomous dynamical system

A one-dimensional example is used to illustrate the procedure of the proposed method and its efficiency. The model is given by a linear autonomous ODE as follows:

(25)

where ξ₁ and ξ₂ are independent random variables distributing uniformly between -1 and 1, i.e., ξ₁∈ U[-1, 1] and ξ₂∈ U[-1, 1]. The exact solution to the above initial value problem is

(26)

which in separated form is

(27)

Here, the separated representation of the stochastic response surface implies that x(t, ξ₁, ξ₂) can be approximated by only one term with the appropriate basis. We demonstrate that the high-precision approximation can be achieved by using the proposed method with only one term expansion.

The mean and variance of the stochastic solution can be, respectively, calculated exactly by

(28)

(29)

With the proposed method, the stochastic response is approximated by the dynamic bases as follows:

(30)

where P_m (t, ξ₁) and Q_n (t, ξ₂) are the dynamical bases, C=(c_ij)_M×N is the coefficient matrix, and P and Q are vectors of the dynamical bases expressed by

(31)

(32)

If the dynamical bases are represented by Legendre polynomials and truncated by the R term, we obtain

(33)

(34)

where A_m=(a_{m, i})_{0≤ i≤R} and B_n=(b_{n, i})_{0≤ i≤R} are row vectors of the coefficients. Then, P(t) and Q(t) can be rewritten as follows:

(35)

(36)

Further, the stochastic response can be expressed by

(37)

where A(t) and B(t) are matrices of the coefficients expressed by

The numerical implementation of this problem is (see Appendix B for its derivation)

(38)

These equations can be evaluated well without the regularized method.

Figure 4 compares the time-dependent stochastic basis approximation of the results with those of the GPC method. It is obvious that the proposed method can achieve a high-precision approximation of the mean and standard deviations with only a first-order expansion. For the GPC method, the error of the variance increases significantly as time evolves. This implies that a high-order GPC expansion is required to yield accurate results.

In Fig. 5, the convergence rates of the mean and variance at t=30 s are compared with different orders. The errors of the proposed method converge much faster than those of the GPC method. Using the time-dependent stochastic bases, we can calculate the accurate results with a low-order expansion.

To further examine the accuracy of the stochastic response surfaces obtained, some snapshots of the approximation errors of the response surfaces are shown in Fig. 6. From Fig. 6, we can see that, the time-dependent stochastic bases yield nearly the exact approximation of the stochastic response surface. Even when the time increases to 30 s, the approximation error is acceptable. This indicates that the proposed method can modify the bases adaptively to approximate the stochastic response surface.

In Fig. 7, the dynamical basis P₁ (t, ξ₁) is plotted and compared with the optimal basis. It can be observed that P₁ (t, ξ₁) is nearly optimal. This means that the time-dependent stochastic bases can track the optimal bases automatically in this example.

4.2 Clamped-clamped Euler-Bernoulli beam

The dynamical response of a clamped-clamped Euler-Bernoulli beam (see Fig. 8) with uncertain damping and Young's modulus is considered in this example. The material of the beam is assumed to be steel. The system is discretized by using the finite element method with 10 elements. The uncertain damping and Young's modulus are assumed to be independent random variables of the form

(39)

where ξ₁ and ξ₂ are independent random variables distributed uniformly in the interval. The excitation force is assumed to be uniform along the axial direction.

(40)

With the proposed method, the stochastic response surface can be evaluated at different time. To examine the efficiency of the proposed method, the response surfaces are compared with those evaluated by the GPC method. The response surfaces computed by different methods are plotted in Fig. 9, in which how time-dependent stochastic bases adaptively approximate the stochastic response as time progresses is clearly demonstrated. The initial response is a plane that can be well approximated with both the GPC method and the proposed method. As time progresses, a peak develops along ξ₁ at approximately t = 1.5 s. The response surface of the GPC begins to diverge from the exact one. From t = 2.5 s to t = 3.5 s, another peak appears along ξ₁, and the GPC fails to approximate the exact response surface. However, the response surface of the proposed method yields a good approximation at different time. The GPC method is typically used to approximate the moments of the stochastic responses, e.g., the mean and standard deviation. From this example, we observe that the GPC method cannot approximate the response surface of the stochastic dynamical system well. However, the proposed method overcomes this drawback by introducing time-dependent stochastic bases.

In Figs. 10 and 11, we plot the cumulative distribution functions (CDFs) of the stochastic response at the middle point of the beam at t = 2 s and 3 s. The reference results are evaluated by using the Monte-Carlo (MC) method with 4 000 samples. As shown in these figures, the proposed method provides better CDFs than the GPC method, especially when a low-order expansion is used. The CDF curves of the GPC method with low order are far away from the real ones, which indicates that the time-dependent stochastic bases can improve the GPC method effectively to yield more accurate CDFs. With the time-dependent stochastic bases, more robust results can be achieved as they are less sensitive to the expansion order.

Figure 12 shows the approximation error of the proposed method and the GPC method, where errors are filtered with a low-pass filter when they oscillate with time. It is obvious that the approximation error by the proposed method is much smaller than that by the GPC method. The error of the GPC method increases significantly with time, whereas our method yields a good control of approximation error. It is obvious that the proposed method can yield a better approximation of the stochastic response surface with a smaller expansion order. As shown in the figure, the error of the proposed method of the 2nd-order is comparable with the GPC method of the 4th-order.

4.3 Beam frame network

In this example, a five-story frame structure with an uncertain Young's modulus is studied to demonstrate the efficiency of the proposed method with high-dimensional uncertainties (see Fig. 13). The structure is analyzed by using the finite element method, and each beam is modeled with three elements. The material properties are specially chosen to illustrate how the proposed method can track the response surface better. The uncertain density and Young's modulus are assumed to be independent random variables as follows:

(41)

The 3rd-order time-dependent stochastic bases are chosen to approximate the stochastic response surface. Figure 14 depicts the mean response of the vertical displacement at Node A. The reference results are evaluated by using the MC method with 4 000 samples.

In Fig. 15, the response surfaces of the vertical displacement of Node A at different selected time are plotted. Initially, the response surface is a plane. As time progresses, it becomes increasingly complicated. This figure illustrates that the proposed method approximates the stochastic response surface adaptively.

4.4 A nonlinear case: stochastic Burgers' equation

This example considers Burgers' equation driven by stochastic forces. Burgers' equation is a well-known nonlinear system. The stochastic Burgers' equation has been studied well by many researchers.

We consider a one-dimensional Burgers' equation as follows:

(42)

where the initial condition u(x, 0, ξ) is deterministic, and the boundary condition is periodic. The stochastic force f is assumed to be a Brownian motion that can be approximated with a suitable orthonormal basis. Following Ref. [33], the stochastic force is expressed by

(43)

where is the spatial basis, and

(44)

in which M_i (t) (i= 1, 2, …) are the orthonormal bases of L² ([0, T]), and ξ_i (i=1, 2, …) are independent Gaussian random variables.

(45)

To simplify the implementation, the order r in Eq. (44) is set to be 3. The viscosity is chosen to be ν=0.005.

We first consider the error of the results. The reference result is evaluated by the collocation method of the 10th-order. In Figs. 16 and 17, the errors of the proposed method of different orders are compared with those of the GPC method. These two figures reveal that the proposed method yields a significant reduction in the computed error. The error decreases with the expansion order. With the time-dependent stochastic bases, smaller errors can be achieved. This indicates that the proposed method can reduce the expansion order, achieve acceptable results, and yield significant computational efficiency for the stochastic problems with many uncertain parameters.

5 Summary

In this study, a method to approximate the stochastic dynamical systems based on time-dependent stochastic bases is developed. The TDBs are modified under the dynamically orthogonal constraint. Under the assumption of orthogonal bases, differential equations of time-dependent stochastic bases are derived by using the variational principle. The method is applied to both linear and nonlinear examples. In the first example, the proposed method could effectively approximate the response surface with only a 1st-order expansion, while the GPC method needs the 4th-order expansion. In the second and third cases, the vibrations of a beam and a beam network are better approximated by the proposed method than by the GPC method with the same order expansion. The time-dependent stochastic bases provide an effective strategy to improve the computational efficiency of the GPC method, especially when more accurate response surfaces are demanded.

Appendix A   Derivation of evolution equations

In this section, the variational principle is used to derive the evolution equations. We must derive the differential equations of the coefficient field and TDBs .

From the orthonormal condition in Eq. (10), we can see that the multi-index bases should satisfy

(A1)

Additionally, from the orthogonal increment condition in Eq. (12), we obtain

(A2)

This implies that the multi-basis Ψ_I and the basis incremental part of (x, t, ξ) are orthogonal, and

(A3)

Equations (A1) and (A3) are important for the derivation of differential equations of time-dependent orthogonal bases. They will be used in the following.

Using the variational principle, the equations for the dynamical bases can be written as follows:

(A4)

where δu is the variation of u, i.e.,

(A5)

and

(A6)

In addition, the variational principle can be written in the optimization form as follows:

(A7)

subject to arbitrary and with . The equations of motion of TDBs are derived from the variational principle in the following part.

First, δu is chosen to be

with arbitrary multi-index I. From Eq. (A4), we have

(A8)

It equals

(A9)

Because it holds for arbitrary δc_I (x, t) and index I, we can obtain

(A10)

Using the orthonormal condition of the bases and orthogonal incremental constraint, we can calculate the derivatives of the coefficients by

(A11)

After is evaluated, the approximation error without basis increments can be written as follows:

(A12)

In the GPC method, Σ_u is the local approximation error. With TDBs, we aim at reducing the approximation error further to improve the precision. It is obvious that Σ_u is orthogonal to the space spanned by the current basis Ψ_κ(t) (κ=1, 2, 3, …, m), i.e.,

(A13)

To evaluate , we choose

From Eq. (A4), we have

(A14)

where

(A15)

Substituting Eq. (A15) into Eq. (A14) yields

(A16)

The orthogonal incremental constraint implies

(A17)

Subsequently, Eq. (A16) can be simplified to

(A18)

Rearranging the above equation yields

(A19)

(A20)

The expectation operator in the above equation is equal to the combination of E_i (·) and P_ξi (·). Therefore, we have

(A21)

(A22)

As it holds for arbitrary δΨ_i^T(t), we obtain

(A23)

To evaluate , we integrate the equation above in the spatial domain. Subsequently, we obtain

(A24)

Subsequently, Eq. (A24) can be written as

(A25)

Finally, the differential equations of TDBs can be evaluated by

(A26)

This completes the derivation of the evolution equation of the dynamical bases.

From the above formulas, we can rewrite the evolution equations of the dynamical bases as follows:

(A27)

Appendix B   Derivation of TDBs in Subsection 3.1

The orthonormal constraints of P and Q imply

(B1)

(B2)

Therefore, A(t) and B(t) are orthogonal matrices. From the orthogonal increment constraint in Eq. (12), we have

(B3)

Using the matrix notation, we can write the original differential equation as follows:

(B4)

Multiplying the equation above by L(ξ₁) from the left and L^T (ξ₂) from the right and subsequently evaluating the expectation yield

(B5)

where

(B6)

(B7)

To evaluate matrices , , and , we first rewrite Eq. (B5)

(B8)

where the dependence on time t is ignored in this equation to simplify the expression. can be evaluated by simply multiplying A(t) on the left and B^T (t) on the right, i.e.,

(B9)

Substituting the equation above into Eq. (B8), we have

(B10)

Multiplying B^T on the right and simplifying the equation, we have

(B11)

Similarly, can be evaluated by

(B12)

Finally, the differential equations to evaluate A(t), B(t), and C(t) can be written as follows:

(B13)