Skip to main content

An accuracy comparison of polynomial chaos type methods for the propagation of uncertainties

Abstract

In (Augustin et al. in European J. Appl. Math. 19:149-190, 2008) we considered the Polynomial Chaos Expansion for the treatment of uncertainties in industrial applications. For many applications the method has been proven to be a computationally superior alternative to Monte Carlo evaluations. In the current overview we compare the accuracy of Polynomial Chaos type methods for the propagation of uncertainties in nonlinear problems and verify them on two examples relevant for industry. For weakly nonlinear time-dependent models, the generalized Kalman filter equations define an efficient method, yielding good approximations if the quantities of interest are restricted to the first two moments of the solution. Secondly, stochastic collocation is discussed. The method is applied to delay differential equations and random ordinary differential equations. Finally, a generalized PC method is discussed which is based on a subdivision of the random space. This approach is even suitable for highly nonlinear models.

1 Introduction

Uncertainty quantification in industrial applications is an interesting and important field of research [14]. Due to uncertainties in system parameters, measurements and in the modeling of physical processes itself, deterministic approaches for the simulation of those processes are not appropriate. Uncertainty quantification became a vast field of research over the past decades and many methods have been proposed to deal with random data [510]. Amongst these methods, approaches based on Polynomial Chaos expansions have proven to be very efficient in dealing with uncertainties, see [3, 10, 11].

Those methods revealed to be efficient alternatives to Monte-Carlo (MC) simulations. They were successfully applied to solve stationary problems e.g. in the field of stochastic finite elements [10] and to solve differential equations with uncertain parameters and initial values [12]. The suitability of PC for industrial applications was considered in [1]. Additionally, an extensive discussion of the Polynomial Chaos approaches, applied to linear differential equations, was presented. In the current paper we consider numerical methods for the propagation of uncertainties through nonlinear differential equations. The main part is concerned with a comparison of Kalman filter techniques for state estimation and Polynomial Chaos Galerkin techniques. The methods themselves are known from literature, the contribution of this work is the comparison of these methods in selected examples relevant for industry. Our main focus is on the accuracies of the methods. Based on these results, we try to answer the question what is the numerical effort in order to achieve a low, medium or high accuracy while, in terms of function evaluations, gaining at least one order compared to a crude Monte Carlo method.

  • The most widespread method for treating differential equations with uncertainties are techniques of Kalman filtering. While the original Kalman filter was developed to identify noisy observations with linear differential equations [13], here we use it only to integrate the underlying system with uncertainties. A direct generalization to nonlinear equations is given by the extended Kalman filter [14]; it is based on the sensitivities of the equation. The unscented Kalman filter, a nonlinear filter which does not use derivative information, was developed in [15]. It directly approximates the covariance matrix and is, in this context, used just for state estimation.

  • Stochastic collocation methods (SCMs) are introduced by [16, 17]. The first practical application of the methods are discussed in [10]. Here, we apply the methods to nonlinear differential equations with delay, which present a more complex behavior than the delay-free counterparts. Stochastic collocation methods have the advantage of using the differential equation as a black box.

  • In [1] we presented a theory of the stochastic Galerkin method (SGM), see [3, 10], for ordinary differential equations with uncertain parameters. Nevertheless, this method suffers from the lack of convergence, especially when time evolves [18]. Thus, we have to consider generalized approaches to scope with this problem. Several methods have been proposed in literature, see for example the multi resolution method [9, 19] and the Multi-Element generalized Polynomial Chaos (MEgPC) method [20]. In [21] a rigorous convergence theory for MEgPC and a fully adaptive scheme in time and random space has been developed. In the current paper we compare this MEgPC approach to the other PC type methods mentioned above.

In Section 2 we introduce two examples from industrial applications, which are the base for the comparison of the discussed methods. The first example concerns crack propagation, while the second considers the quorum sensing in biofilms. Next, we present a summary of the Kalman filter equations for state estimation in Section 3. Especially, nonlinear filters in the context of uncertainty propagation of state are discussed. In Section 4 we consider the SCM, which we apply to the random delay differential equations of quorum sensing. We discuss the multi-element approach of the stochastic Galerkin method in Section 5. Finally we compare the three described methods on the base of the two applications, Sections 6.2 and 6.3. In Section 7 we close this article by drawing our conclusions.

2 Test sets of (delay) differential equations

In this section, we introduce the applications considered for the comparison of the methods of the Kalman filter, the SCM and the SGM.

  • The first application deals with the growth of cracks, which is important in the context of life cycle analysis. Industrial applications are for example crack growth in turbine blades or in train wheelsets. Here the task is to define appropriate inspection schemes in order to prevent severe damages. The crack growth is modelled by a highly nonlinear ordinary differential equation.

  • The second example is a model of the quorum sensing of biofilms, a delay differential equation. The control of biofilms requires a deep understanding of the structured community of microorganisms living on inert surfaces. Industrial applications range from sewer systems, fresh water systems to corrosion prevention in cooling flow networks of (electric) power plants.

We restrict ourselves to a low number of uncertain variables. Nevertheless, due to the nonlinearities and the delay, differences between the methods with respect to their achieved accuracies will become visible. A more complex model of crack growth, for the evaluation of inspection schemes for train wheelsets, has been computed with the unscented approximation in [22]. The biofilm model, in particular the physical background, has been studied in more detail in [23].

2.1 Crack growth example

For the modeling of a growing crack we use a simplified model of fracture mechanics, see Harris [24], where a center cracked plate in tension is studied. The crack growth relation is based on a modified form of the Forman relation [25]. The crack size a is governed by the ordinary differential equation in the abstract form

a N =f(a, K I c , Δ σ , C F ),a(0)= a 0 ,0N N L
(1)

where N denotes the number of load cycles and N L the life cycle. The parameters are the stress intensity factor K I c , the stress range Δ σ , the fatigue growth coefficient C F and the initial crack size a 0 . In more detail the crack growth differential equation reads [24]

a N = C F ( ( 1 f ) Δ K 1 R ) n f ( 1 4 π Δ K 0 Δ K tan 1 ( 1 R ) ) p ( 1 Δ K ( 1 R ) K I c ) q ,
(2)

where

ΔK= Δ σ π a Y ( 2 a h ) ,Y(α)= 1 0.5 α + 0.370 α 2 0.044 α 3 1 α

and

R is the stress ratio between maximum and minimum stress, h denotes the width of the plate and n f , Δ K 0 , p, q, α , s are curve-fitting constants.

2.2 The biofilm model with quorum sensing

The model for antibiotics action on a biofilm follows closely the model in [26]. Since the antibiotics are not produced by the hosting organism, it must be supplied externally. Thus, our equations read

A ˙ = S A μ A A ( t ) , B ˙ = α 20 ( 1 B ( t ) σ ) B ( t ) + α 2 B ( t ) B ( t τ ) α 3 B ( t ) A ( t ) .
(3)

Here, A denotes the concentration of antibiotics and B the concentration of bacteria. Moreover, the parameters are the logistic reproductive rate of the bacteria α 20 , the effective carrying capacity of the environment σ, the reproductive rate due to quorum sensing, α 2 , the death rate of bacteria due to antibiotics S A and death rate for the antibiotics μ A . All these parameters are taken to be positive. The differential equation for B includes a standard logistic term

α 20 ( 1 B ( t ) σ ) B(t),
(4)

and a quorum sensing (QS) term

α 2 B(t)B(tτ).
(5)

It represents the extra increase in population after the QS mechanism has been initiated. The process is not instantaneous. Thus, assuming it takes a certain amount of time τ to produce the signal molecules, the population of bacteria at time t, B(t), will receive these signal molecules from the population of bacteria at time tτ, i.e. B(tτ). The coefficient α 2 accounts for the strength of the mechanism. The larger α 2 , the more efficient the QS mechanism takes place in the biofilm.

3 Solving stochastic differential equations by Kalman filters

In this section we consider stochastic dynamical systems of the form

(6)
(7)

with 0t T L ,

Here, f, h are linear or nonlinear functions, x is the state vector, u is the control vector and y is the observation vector. We assume the process noise w and the measurement noise v to be white, stochastically independent and

wN ( 0 , R + n w ) ,vN ( 0 , R + n v ) .
(8)

For real processes, white noise is a somewhat simplified assumption. There are extensions of Kalman filtering to colored noise (e.g. [27]) but this will not be picked out here. Equations (6) and (7) describe a stochastic system which is continuous in time. In the following, we restrict ourselves to methods for discrete time systems which are obtained by time discretization of (6). Furthermore we concentrate on non-stiff systems. Of course for the numerical integration of stiff systems semi-implicit or implicit methods are required which means a much higher numerical effort. We begin with linear discrete time systems in Section 3.1. Thereafter, we discuss how nonlinear differential equations can be treated by the extended Kalman filter and the unscented Kalman filter in Sections 3.2 and 3.3. After a comparison of both filters in Section 3.4 we give a reinterpretation of the unscented Kalman filter as a sampling method, see Section 3.5.

3.1 Linear discrete time systems

This section adopts the presentation in [14]. A stochastic linear discrete time system can be written as

with w k 1 N(0, Q k 1 ) and v k N(0, R k ) being independent normally distributed random variables. The predictor step of the Kalman filter is given by

x ˆ k = A k 1 x ˆ k 1 + B k 1 u k 1 .

For the a priori and a posteriori estimate errors e k := x k x ˆ k and e k := x k x ˆ k with corresponding covariances P k :=E( e k e k T ) and P k :=E( e k e k T ) holds

e k = A k 1 ( x k 1 x ˆ k 1 )+ W k 1 w k 1 and P k = A k 1 P k 1 A k 1 T + W k 1 Q k 1 .

The corrector of the Kalman filter is given by

x ˆ k = x ˆ k + K k ( y k H k x ˆ k ) ,

thus for the a posteriori error and its covariances holds

e k = x k x ˆ k = x k x ˆ k K k ( y k H k x ˆ k ) = x k x ˆ k K k ( H k x k + V k v k H k x ˆ k ) = ( I K k H k ) ( x k x ˆ k ) K k V k v k

and P k =(I K k H k ) P k ( I K k H k ) T + K k V k R k V k T K k T . Minimization of the trace of P K leads to

K k = P k H k T [ H k P k H k T + V k R k V k T ] 1 .
(9)

With this expression for K k , we derive a simpler formula for P k ,

P k =(I K k H k ) P k .
(10)

The algorithm is summarized in Table 1.

Table 1 Kalman filter

3.2 Extended Kalman Filter (EKF)

The stochastic nonlinear discrete time system can be written as

The predictor of the extended Kalman filter is given by

x ˆ k =f( x ˆ k 1 , u k 1 ,0)

with the a priori estimation error and its covariance being

Here we use the abbreviations

A k 1 := f x | ( x ˆ k 1 , u k 1 , 0 ) and W k 1 := f w | ( x ˆ k 1 , u k 1 , 0 ) .

The corrector of the extended Kalman filter has the form

x ˆ k = x ˆ k + K k ( y k h ( x ˆ k , 0 ) )

and we write the a posteriori estimation error as

e k = x k x ˆ k = x k x ˆ k K k ( y k h ( x ˆ k , 0 ) ) = x k x ˆ k K k ( h ( x k , v k ) h ( x ˆ k , 0 ) ) = ˙ x k x ˆ k K k ( H k ( x k x ˆ k ) + V k v k ) = ( I K k H k ) ( x k x ˆ k ) K k V k v k

with

H k := h x | ( x ˆ k , 0 ) and V k := h v | ( x ˆ k , 0 ) .

For the Kalman gain and the a posteriori error covariance the same formulae (9), (10) as in the linear case apply. The corresponding algorithm of the extended Kalman filter is summarized in Table 2.

Table 2 Extended Kalman filter

3.3 Unscented Kalman Filter (UKF)

For the extended Kalman filter, derivatives of f and h have to be provided which is not possible in many practical applications. So in [15] a nonlinear filter, called unscented filter, without this drawback was introduced. The unscented filter can be used to estimate the mean and covariance of a nonlinear stochastic process f(w) where w R n w is a normally distributed random vector with mean E(w) and covariance P w w R n w × n w . So-called sigma points X ( i ) , together with weights W i mean and W i cov , are constructed and mapped to Z ( i ) =f( X ( i ) ) for i=0,,p. The unscented filter then yields an approximation of the mean μ and covariance P z z of the nonlinear function by

By Taylor expansion one can show second order accuracy of mean and covariance, [15, 28].

Example A set of p+1=2n+1 weights and sigma points with respect to mean E(x) R n and covariance P x x R n × n is given by [29]

for i=1,,n respectively, where ( P x x ) i denotes the i-th column of the matrix root L of P given by P x x =L L T . Thus, the corresponding weights are W i mean = W i cov = W i for i=1,,2n.

Subsequently, we apply the unscented Kalman filter to the recursive estimation of dynamical processes.

3.3.1 Additive process and measurement noise

For the sake of simplicity we begin with the special case of additive process and measurement noise

At first, we construct a set of sigma points X k 1 ( i ) with mean x ˆ k 1 and covariance P k 1 . The dependence of the sigma points on x ˆ k 1 and P k 1 is denoted by X k 1 ( i ) = X k 1 ( i ) ( x ˆ k 1 , P k 1 ). The predictor of the unscented Kalman filter and the error covariance are given by

(11)
(12)

The predicted observation is

(13)
(14)

with innovation covariance and cross covariance

P k y y = i = 0 p W i cov ( Y k ( i ) y ˆ k ) ( Y k ( i ) y ˆ k ) T + R k , P k x y = i = 0 p W i cov ( X k ( i ) x ˆ k ) ( Y k ( i ) y ˆ k ) T .
(15)

The Kalman gain has the form

K k = P k x y ( P k y y ) 1 .
(16)

The corrector of the unscented Kalman filter and the a posteriori covariance are given by

x ˆ k = x ˆ k + K k ( y k y ˆ k ) , P k = P k K k P k y y K k T .
(17)

The relationship between the Kalman gains K k of the extended and unscented filters (9), (16) is given by

P k x y = P k H k T and P k y y = H k P k H k T + R k .

Also P k in (17) corresponds to (10) because substitution of (16) in (17) yields

P k = P k K k P k y y K k T = P k P k x y ( P k y y ) 1 ( P k x y ) T = P k K k H k P k = ( I K k H k ) P k .

The algorithm is summarized in Table 3.

Table 3 Unscented Kalman filter

3.3.2 Nonlinear noise

For the general case of nonlinear noise augmented sets of sigma points for prediction and corrector steps have to be defined [29]. Let

(18)
(19)
(20)
(21)

For the predictor step, define sigma points X ¯ k 1 ( i ) with respect to mean

x ¯ k 1 :=( x ˆ k 1 0 )
(22)

and covariance

P ¯ k 1 =( P k 1 0 0 Q k 1 ).
(23)

Then (11) is substituted by

X k ( i ) = f ¯ ( X ¯ k 1 ( i ) , u k 1 ) ,i=0,,p
(24)

and the Q-term in (12) vanishes. Accordingly define sigma points X k ¯ ¯ ( i ) for the corrector step with mean

x ¯ ¯ k :=( x ˆ k 0 )
(25)

and covariance

P ¯ ¯ k =( P k 0 0 R k ).
(26)

Then (13) is substituted by

Y k ( i ) = h ¯ ( X k ¯ ¯ ( i ) ) ,i=0,,p
(27)

and the R-term in (15) vanishes. The algorithm is summarized in Table 4.

Table 4 General unscented Kalman Filter

3.4 Comparison of EKF and UKF in crack growth example

Figure 1 shows a comparison of the extended Kalman filter, the unscented Kalman filter and a Monte Carlo simulation for the estimation of mean and covariance of the time dependent crack depth. In this test, we only assume the initial crack size a 0 to be uncertain, whereas we set the other parameters in (1) to deterministic values. It can be seen from the figures that both extended and unscented filter are fairly accurate in the linear regime of the crack growth, whereas in the nonlinear regime only the unscented filter shows sufficient accuracy.

Figure 1
figure 1

Comparison of estimated mean values (Upper) and covariances (Lower) of crack depth.

3.5 Reinterpretation of UKF

The extended Kalman filter does not treat the nonlinearities of the above example and is therefore no longer considered. The Kalman filter techniques are developed to deal with uncertainties in the state and the measurements. In order to allow a comparison to the methods based on Polynomial Chaos expansions, cf. Sections 4 and 5, we omit the treatment of measurements and only consider the state estimation. Thus, the resulting equations that we are considering within the following are ordinary differential equations with random parameters, cf. (30) and (33). In this setting, we interpret the construction of sigma points as special stochastic collocation method which is quadratic in each dimension (28). We reconstruct the solution between the sampled points by the polynomial

x( T L , θ 1 ,, θ n w )= i = 1 n w ( q 0 , i + q 1 , i θ i + q 2 , i ( θ i 2 1 ) ) .
(28)

In a post-processing step, we will apply a Monte-Carlo simulation based on this polynomial approximation, to obtain approximations of any quantity of interest, i.e. the moments or density of the solution.

4 Stochastic Collocation method and random delay differential equations

We will consider the delay differential equation (DDE) with one single, time-independent delay, defined as

x ˙ (t)=f ( x ( t ) , x ( t τ ) ) ,t>0,

with f: R n × R n R n and τ:R[0,+]. The parameter τ is called the delay. Given the initial data set (also called history function) x(t)=ϕ(t), t[τ,0] we obtain the initial value problem (IVP)

IVP:={ x ˙ ( t ) = f ( x ( t ) , x ( t τ ) ) , t > 0 , x ( t ) = ϕ ( t ) , t [ τ , 0 ] .
(29)

When one or more parameters in (29) are affected by uncertainty we will be speaking of a random delay differential equation (RDDE)

x ˙ (t,w)=f ( x ( t , w ) , x ( t τ , w ) , w ) ,t>0,
(30)

if w L 2 . Eventually, the initial history function can also exhibit a random behavior, x(t,w)=ϕ(t,w), t[τ,0].

In order to account for the uncertainties the stochastic collocation schema of the Polynomial Chaos Expansion is used. The schema is explained below, along with the numerical method applied to solve for the underlying deterministic equations.

4.1 Wiener expansions

We will denote a probability space by (Ω,F,P) and the space of square integrable random variables on it by L 2 (Ω,F,P). We will abbreviate it by L 2 in the forthcoming. This space is equipped with the inner product

(w,v):= Ω w(θ)v(θ)dP(θ)

and the induced norm w:= ( w , w ) . Let ξ be a standard Gaussian distributed random variable on (Ω,F,P), where is generated by ξ. N. Wiener [16] and later R. Ghanem and P. Spanos [10] showed, that every random variable x L 2 can be expanded into a series

x ( w ( ) ) = i = 0 q i Ψ i ( w ( ) )

with basis polynomials { Ψ i } i = 0 . A detailed discussion on this expansions can be found in [1, 4, 30]. In [11] D. Xiu and G. Karniadakis generalized this expansion to polynomials of the Askey scheme. For convergence results of these generalized expansions we refer to [31].

The Wiener expansions can be used to solve operator equations with square integrable random parameters. Assume the operator T:X R n , mapping from the Banach space X into R n , given by

T(x;w)=0,
(31)

with xX and parameters w L 2 . If the solution x is in L 2 , it can be expanded by a generalized Wiener expansion. A truncation of this expansion after p ˆ +1 terms results in an approximation for x in the finite dimensional subspace linspan{ Ψ i :i=0,, p ˆ } L 2

x(w) i = 0 p ˆ q i Ψ i (w).
(32)

Thus the solution x can be approximated by computing the coefficients { q i } i = 0 p ˆ . Different methods can be used to calculate the coefficients. We opt here for the stochastic collocation.

4.2 Stochastic collocation

Stochastic collocation methods make use of the fact that the Wiener expansions are generalized Fourier series. Therefore, the coefficients in (32) are given by

q i = ( x , Ψ i ) ( Ψ i , Ψ i ) ,i=0,, p ˆ .

For the computation of the coefficients { q i } i = 0 p ˆ , an estimation of the inner products is necessary. This can be accomplished with cubature rules, which approximate the multidimensional integral as a discrete sum

(x, Ψ i ) j = 1 N Q x ( w ( j ) ) Ψ i ( w ( j ) ) W ( j ) ,

where w ( j ) , meant as w( θ ( j ) ), and W ( j ) are the nodes and weights of the N Q -node cubature and are given for each type of cubature, see [4] for more details.

The procedure for the stochastic collocation scheme would be then as follows: For each of the given nodes calculate the solution of (31) and sum them up with the corresponding weights. The products ( Ψ i , Ψ i ) are independent from the operator equations and can be conveniently pre-computed and stored for efficient use in different calculations.

This approach offers some advantages. Contrary to MC techniques, the sampling is not random, but obey the cubature rules to minimize the needed number of nodes for a given precision. Thus, the amount of realizations is significantly reduced. Usually, 10 nodes per dimension are accurate enough. Nevertheless, for increasing number of random variables in the parametrization, the numbers of nodes grows exponentially (curse of dimensionality), and the dimension-independent MC becomes eventually more efficient. Sparse grid cubatures have been proposed to mitigate the curse of dimensionality [4].

Another important feature that makes the stochastic collocation schema attractive within the Wiener expansion methodology is the splitting of the uncertainty quantification and the system governing equations. Since the values of the random variables are chosen in advance, we only need to solve the underlying deterministic equations for the given value of the parameters. Hence, no modification of our existing deterministic software is needed in order to account for uncertainties. This is specially valuable in the case of commercial programs, which cannot be modified, or complex software, whose adaptation for uncertainty quantification can be error-prone. A further advantage is the modularity, making implementation of refinements in the Wiener expansion algorithms and the underlying governing equations independent.

For solving the RDDEs (30) we use the stochastic collocation schema combined with the program RADAR for the computation of the DDEs.

4.3 Numerical solution of the DDEs

The software RADAR is used as the deterministic solver for DDEs. It has been developed by Ernst Hairer and Nicola Guglielmi. For more details a manual with examples is available [32].

In a simplified form for presentation purposes, the DDE

with α(t,x(t)) the (eventually state-dependent) time-lag function, is solved using an implicit Runge-Kutta method called RADAU IIA [33]. If we denote the stepsize by h k = t k + 1 t k and the coefficient matrix by A= ( a i j ) i , j = 1 s , c i = j = 1 s a i j , weights { b i } i = 1 s the implicit schema reads

Y i ( k ) x k = h k j = 1 s a i j f ( Y j ( k ) , Z j ( k ) ) ,

where x k is the approximation of x( t k ) and Y s = x k + 1 , if we use the weights b i = a s , i , i=1,,s. Z j ( k ) is defined by

Z j ( k ) ={ g ( α j ( k ) ) , α j ( k ) < t 0 , u m ( α j ( k ) ) , α j ( k ) [ t m , t m + 1 ] ,

where α j ( k ) =α( t k + c j h k , Y j ( k ) ) and u m is a polynomial approximation of the solution x(t) on the interval [ t m , t m + 1 ]. Here, the Lagrange interpolation formula is used to build the approximation.

RADAR allows steps larger and smaller than the delay and performs stepsize control by estimating the local error at grid points and the error in the continuous numerical approximation to the solution.

5 Stochastic Galerkin method and random ordinary differential equations

Within Section 4 we discussed the treatment of RDDEs by the collocation method based on generalized Wiener expansions. In this section we will present the stochastic Galerkin method (SGM) for the approximation of the solution of differential equations with uncertain parameters, cf. [1]. We again consider the differential equation (30), x ˙ (t,w)=f(x(t,w),x(tτ,w),w), and omit the delay within the forthcoming discussion to keep the notation simple:

x ˙ (t,w)=f ( x ( t , w ) , w )
(33)

with initial condition x(0,w) L 2 . The extension of the SGM including the treatment of delays is straightforward using interpolation for the approximation of the solution at past times. In Section 6 we will apply the SGM to the test example of the crack growth (an RODE) and the biofilm model (an RDDE) as stated in the introductory Section 1. Subsequently, we intend to present the main idea of the SGM and restate the main results of its convergence analysis from [21, 34].

For the numerical approximation of the solution of (33) we search for an approximation of x(t,w) within the finite dimensional subspace

S=linspan{ Ψ i :i=0,, p ˆ } L 2 ,

the approximation space. Without loss of generality we assume that the orthogonal polynomials { Ψ i } i = 0 are normalized, i.e. ( Ψ i , Ψ i )=1.

Besides RODE (33) we introduce the perturbed RODE

t x ˜ (t,w)=f ( π p ˆ x ˜ ( t , w ) , w ) =: f ˜ ( x ˜ ( t , w ) , w ) ,
(34)

where π p ˆ : L 2 S denotes the orthogonal projection from L 2 onto S. We define the restricted stochastic weak form of (34) by

( t x ˜ ( t , w ) , Ψ ( w ) ) = ( f ˜ ( x ˜ ( t , w ) , w ) , Ψ ( w ) )
(35)

for all ΨS. Equation (35) can be written as a system of ordinary differential equations (ODEs) for the coefficients { q ˜ i } i p ˆ :

d d t q ˜ i (t)=( f ˜ , Ψ i )for i p ˆ .
(36)

This system of ODEs is deterministic and can be solved by standard integrators. The projected system (36) can be solved by means of a Runge-Kutta method (RKM) of order p1 and stepsize h>0. We obtain the approximation x ˆ of the solution of (34) in terms of the numerically computed coefficients { q ˆ i } i p ˆ using the SGM and a RKM,

x ˜ (t,w) x ˆ (t,w):= i p ˆ q ˆ i (t) Ψ i (w).

The following scheme depicts the stochastic Galerkin Runge-Kutta method (SGRKM):

t x ˜ = f ˜ ( x ˜ , w ) SGM d d t q ˜ = F ( q ˜ ) RKM x ˆ = i p ˆ q ˆ i Ψ i ( w ) .
(37)

In [34] it is shown that the local discretization error of the SGRKM does not vanish for decreasing stepsize h0. Thus, this global approach of the stochastic Galerkin method fails to converge to the exact solution x. To overcome this problem a generalized approach has to be applied. X. Wan and G.E. Karniadakis proposed the Multi-Element generalized Polynomial Chaos (MEgPC) method to attain convergence of the SGM. We recall that this method is based on a partitioning of the parameter space Ω into subsets, so called random elements. The SGM is then applied locally on the partition of the parameter space, i.e. the problem is conditionally restricted to each random element. We refer to the literature for a detailed description [18].

Using the MEgPC approach it can be shown that the local discretization error of the SGRKM obeys

x ( t + h , w ) x ˆ ( t + h , w ) 2 O ( h Δ p ˆ + 1 ) +O ( h p + 1 ) ,
(38)

where we denote the maximal diameter of the stochastic elements by Δ, see [34].

5.1 The adaptive algorithm for the solution of RODEs

We use (38) to control the stepsize h in time and the maximal size Δ of the stochastic elements { Ω i } i = 1 ν . In a first step, we compute the conditionally restricted solution x ˜ (t+h,w) on every element Ω i and check the local truncation error of the projected restricted system (36). See [35] for details on adaptive stepsize control of RKMs. This procedure results in the term O( h p + 1 ) being negligible up to a prescribed tolerance. In a second step, we check the error ERR, which is introduced by the truncation of x ˜ (t+h,w) to S. Therefore, we split Ω i into two equal subsets and compute the solution of (36) conditionally restricted to those subsets. We take the L 2 -norm of the difference of the refined solution and the solution x ˜ (t+h,w) on Ω i as substitute for the actual truncation error ERR i . We apply this procedure in every random element at every timestep, so we omit the corresponding index i of the subsets in the following. According to (38) it holds

ERR=ch Δ p ˆ + 1 ,
(39)

where we neglect the error term O( h p + 1 ) of the RKM. Furthermore, we require the optimal size Δ opt of the respective element to satisfy

ε=c Δ opt p ˆ + 1
(40)

for prescribed tolerance ε>0. Division of (40) and (39) results in

Δ opt = ε h ERR p ˆ + 1 Δ=:I(ERR)Δ.

We use the estimator I(ERR) to decide about the refinement of the elements { Ω i } i = 1 ν . If I(ERR)< ρ 1 , then the element will be refined. If I(ERR) exceeds a coarsening threshold ρ 2 > ρ 1 , then the element will be coarsened. For a more detailed description of this procedure see [34].

6 Computational aspects and numerical results

In this section, we apply the three methods, the UKF approach, the stochastic collocation and the SGRKM, to the benchmark problems stated in the introduction. From the solution we compute the densities and the first two moments and compare the results. For the numerical time integration in the SGM we use an embedded explicit RKM of order 4(3). After a discussion of the efficiencies of the methods, we first consider the crack growth example (1) with the initial crack size being uncertain. Thereafter, we study the setting of two uncertain parameters in the model of crack growth. We finish this section by the comparison of the three methods applied to the biofilm model (3).

6.1 Computational aspects

For UKF and SCM, the efficiency in terms of function evaluations can directly be given as function of the stochastic dimension n (Table 5). Here the number of sigma points for UKF has been chosen as in the example Section 3.3. For SCM a standard Gaussian cubature rule of order s has been used. SCM and MEgPC show an exponential growth in the number of function evaluations (suitable for 2D problems). Sparse cubature rules which are necessary to scope with high dimensional random spaces are considered in [22]. The MEgPC approach results in ν projected systems of ordinary differential equations (36) of dimension

( n + p ˆ ) ! n ! p ˆ !

if a polynomial basis of total order p ˆ is applied. The number of needed function evaluations depends on the cubature rule applied to compute the projection integrals, for example using a Gaussian cubature tensor product rule we need ( s + 1 ) n function evaluations per timestep of the time integration method. However, the MEgPC approach is most efficient, if the projection integrals can be computed analytically in a preprocessing step before the time-integration starts. If this so called intrusive approach is not feasible, we have to resort to cubature rules to compute the projection. In the case of multiple random dimensions, n2, sparse integration methods have to be applied in order to reduce the curse of dimensionality.

Table 5 Number of function evaluations for UKF, SCM and MEgPC approach based on a Gaussian cubature rule of order s

6.2 Comparison of UKF, Stochastic Collocation and SGM in crack growth example

In this section, we compare the presented methods and Monte Carlo Sampling (106 samples) in the crack example (1).

6.2.1 One dimensional random space

At first, we assume only the initial crack size a 0 to be a random parameter following a uniform distribution with expected value 10−1 and variance 10−6. The other parameters are set to K I c = 10 12 , Δ σ =16 and C F =0. We remark that although the right-hand side function f in (1) is not polynomial in the initial value a 0 , the convergence results from Section 5 apply, see [21, 34] for a detailed discussion.

To compute the numerical approximation a ˆ (t, a 0 ) of the crack size a(t, a 0 ) by means of the MEgPC method, we use the parameters stated in Table 6. The adaptive MEgPC approach uses a truncated generalized Wiener expansion of degree p ˆ =3 for the approximation of the solution. The set of outcomes Ω is refined where the local error exceeds the prescribed tolerance ε= 10 8 . The approximated response surfaces at t{5,000,5,300,5,400,5,420} are shown in Figure 2. The respective elements of the approximation are indicated by dots on the axis of the outcomes Ω.

Figure 2
figure 2

Response surfaces of the solution a ˆ (t, a 0 ) at t{5,000,5,300,5,400,5,420} .

Table 6 Setting of parameters of the adaptive MEgPC approach for the solution of RODE ( 1 )

Although the initial crack size varies only by a small amount, the impact on the outcome is not negligible. A large variation in the crack size occurs at the final time t=5,420. The dependence of the solution a(t, a 0 ) on the uncertain initial crack size a 0 is highly non-linear. This is recognized by the MEgPC algorithm, which refines the parameter space where the response surface can not be represented by a polynomial of degree p ˆ =3 up to the prescribed tolerance.

In Figure 3 we compare the densities of the approximations a ˆ (t, a 0 ), t{5,000,5,420}, computed by the stochastic collocation of order 5 (SCM5), MEgPC and the UKF approach respectively. We also plot the density computed by Monte Carlo simulation as a reference for the approximatively exact density. On the one hand, we observe that the UKF is not able to compute a good approximation of the density of the crack size a. This can be explained by the interpretation of the UKF as a stochastic collocation method of polynomial order 2, see Section 3.5. On the other hand, the MEgPC method, as we expect from our considerations in Section 5, yields a good approximation of the density of the crack size a in dependence on an uncertain initial crack size.

Figure 3
figure 3

Density of the crack size at two final times t{5,000,5,420} , assuming the initial crack size to be an uncertain parameter. Comparison of four different methods: Monte Carlo simulation, UKF, SCM5 and MEgPC.

In addition to the comparison of the densities we also compare the first two moments, see Table 7. We remark that although the UKF approach does not capture the overall density, it is able to approximate the first two moments in good accuracy. Moreover, the first two moments, which we computed by means of the MEgPC method and the Monte Carlo simulation match very good.

Table 7 The expected values and variances of a(t, a 0 ) at times t{5,000,5,420} computed by SCM5, UKF, MEgPC and Monte Carlo simulation

6.2.2 Two dimensional random space

In this section we consider the crack example (1) with two uncertain parameters. In addition to Section 6.2.1 we assume the parameter K I c to be uniformly distributed with expected value 0.25 and variance 0.02.

To compute the numerical approximation a ˆ (t, a 0 , K I c ) of the crack size a(t, a 0 , K I c ) by means of the MEgPC method, we use the parameters as stated in Table 8. The notation p ˆ =(3,3) means that we use a polynomial approximation of degree 3 in both random variables a 0 and K I c .

Table 8 Setting of parameters of the adaptive MEgPC approach for the solution of RODE ( 1 )

In Figure 4 we plot the densities computed by SCM5, UKF and MEgPC at times t{3,000,3,500,3,550}. Additionally, we plot the density of the reference solution computed by Monte Carlo simulation. We see that the difference between the density computed by the UKF and the reference density is negligible at t=3,000 and becomes larger as time evolves. However, the overall structure of the density can be captured by the UKF approach. This is also reflected by the good approximation of the first two moments, see Table 9. Nevertheless, due to the increasing non-linear dependence on the crack size to the uncertain parameters the accuracy of the UKF approach becomes worse with increasing time. Again, this can be explained by the interpretation of the UKF as a polynomial approximation of degree 2 of the solution, see Section 3.5. Similar to the case of only one random parameter the densities computed by the MEgPC method and Monte Carlo simulation show perfect agreement. Moreover, the first two moments of the solutions agree very well, see Table 9.

Figure 4
figure 4

Density of the crack size at 4 final times resulting from two uncertain parameters. Comparison of 3 different methods: Monte-Carlo simulation, UKF, SCM5 and MEgPC. Left: complete density, right: zoom into the density.

Table 9 The expected values and variances of a(t, a 0 , K I c ) at times t{3,000,3,500,3,550} computed by SCM5, UKF, MEgPC and Monte Carlo simulation based on 10 6 samples

6.3 Comparison of UKF, Stochastic Collocation and SGM in biofilm example

The introductory model for the biofilm (3) with delay Quorum Sensing has been studied in [36] in great detail. Here we compare the numerical approximations of the solution by means of its densities and first two moments. Therefore, we compute them by the UKF approach, the stochastic collocation and the MEgPC method at the time t=100. The uncertain parameter is assumed to be the delay τ, a beta(2,2) distributed random variable with expected value 4.8. Moreover, we use the choice of parameters α 20 =0.25, σ= α 2 = α 3 =0.6, S A =0.4 and μ A =0.5 as well as the initial populations B(0)=1 and A(0)=0.

For the MEgPC method we use the choice of parameters as stated in Table 10.

Table 10 Setting of parameters of the adaptive MEgPC approach for the biofilm model

As indicated by the crosses we do not use the adaptive mesh refinement strategy in this example. The reason is that we have to evaluate the solution at past timesteps in order to get an evaluation of the right-hand side function. In order to evaluate the solution at a time in the past the interpolation of the solution between two timesteps is needed. Therefore, we choose to have the solution defined on the same mesh of elements at every timestep. In order to avoid different meshes at different timesteps we switch off the mesh adaptivity and prescribe a sufficiently fine mesh, see the dotted lines in Figure 5, from the beginning of the computations. The response surface of the solution at time t=100, computed by the MEgPC approach, is shown in Figure 5.

Figure 5
figure 5

Response surface of the bacteria concentration B(100,τ) .

From the solution of the 3 methods we compute the density at time t=100, see Figure 6 (left). The right plot shows a magnification of the density to the tail of the density, i.e. to the interval [3,7]. On the one hand, we see that the density computed by the stochastic collocation of order 5 simulation and the MEgPC method match very good (therefore we here omit the Monte-Carlo solution) . On the other hand, we again observe a difference to the density computed by the UKF approach. The UKF density shows a much smaller weight in [0,0.5] than the two other methods. This can be explained by the polynomial approximation of degree 2, where more mass of the response surface lies within this range [0,0.5] of outcomes. Nevertheless, the overall shape of the density can be captured well by the UKF. In Table 11 we state the first two moments, respectively computed by the three methods. We see good agreement of the first two moments.

Figure 6
figure 6

Density (left) of the biofilm model at the final time t=100 resulting from one uncertain parameter. Comparison of 3 different methods: SCM, UKF and MEgPC approach. Left: complete density, right: zoom into the density.

Table 11 The expected values and variances of the biofilm model at time t=100 computed by the MEgPC approach, the UKF and the stochastic collocation method

7 Conclusions

In the present article, we discussed the numerical treatment of problems with uncertain parameters. We gave a detailed introduction to the unscented Kalman filter and its application to random ordinary differential and random delay differential equations, see Section 3. We dropped the measurements in the equations and model the uncertainty by random variables. In this setting, the unscented Kalman filter can be interpreted as a stochastic collocation method of order 2, see Section 3.5. Sampling of the resulting surrogate polynomial model yields an approximation of the moments of the solution. Due to the fact that the sampling is restricted to the directions given by the columns of the root of the covariance matrix P x x , only 2n+1 realizations of the underlying random differential equation (RDE), resp. an RDDE or an RODE, are needed in the UKF approach. Here, n denotes the number of uncertain parameters. Thus, it is a cheap, in terms of computational costs, method to approximate the random behavior of the solution.

We compare this simplified UKF approach with the SCM and the MEgPC method. In the SCM more realizations, as compared to the UKF, of the solution of the RDE have to be computed. The number of realizations depends on the order of the approximation and the underlying cubature method. For most cubature rules, even in sparse grid approximations, the number of function evaluations is considerably greater than 2n+1 and therefore the SCM is computationally more expensive. The SGM, as opposed to the SCM and the UKF, is not based on sampling, but on a spectral expansion of the solution of the RDE into orthogonal polynomials. An approximate solution is found by a Galerkin approach which results in the stochastic weak form (34). If we are able to compute it in a pre-processing step, only a system of deterministic differential equations for the coefficients of the truncated generalized Wiener expansion has to be solved. The dimension of the projected RDE is small, because we only need a low order approximation within the respective stochastic elements. The overall computational costs scale linearly with the number of stochastic elements. In the present article we did not take advantage of the pre-processing of the projection in the stochastic weak form. On the one hand, the projection in the examples from Section 6 can be efficiently performed by an adaptive Simpson cubature rule. On the other hand, pre-processing the stochastic weak form would result in an projection error, which can not be easily estimated a priori.

In Section 6 we applied the three methods to the benchmark problems of the crack growth model (1) and the biofilm model (3). In the case of the crack growth, which is an RODE, we discussed the prescription of up to two random parameters. It revealed that the UKF approach results in a non-sufficient representation of the overall density, but gives a good approximation of the first two moments. Opposed to that, the MEgPC method computes a good approximation of the density and of the first two moments. This is in agreement of the convergence results in Section 5. In the second benchmark problem, the example of the biofilm (3), we again compared the densities of the solutions computed by the three methods. The SCM, as well as the SGM method showed perfect agreement of the results - in the density and in the first two moments. Again, opposed to that the UKF approach yielded only a rough approximation of the overall density. Nevertheless, it is able to compute a good approximation of the first two moments.

From the numerical results in Section 6 we see that the UKF approach is a cheap alternative to get a rough impression of the random behavior of the solution of the RDE. If detailed information about the solution process is needed, more sophisticated methods like SCM and multi-element SGM have to be applied. This is for example in failure detection the case, where the tail of the probability distribution is of great interest.

References

  1. Augustin F, Gilg A, Paffrath M, Rentrop P, Wever U: Polynomial chaos for the approximation of uncertainties: chances and limits. Eur. J. Appl. Math. 2008, 19: 149–190.

    Article  MATH  MathSciNet  Google Scholar 

  2. Xiu D: Numerical Methods for Stochastic Computations: A Spectral Method Approach. Princeton University Press, Princeton; 2010.

    Google Scholar 

  3. Karniadakis GE, Su CH, Xiu D, Lucor D, Schwab C, Todor RA: Generalized polynomial chaos solution for differential equations with random inputs. Report Nr. 2005–01, ETH Zürich. Seminar für Angewandte Mathematik; 2005.

    Google Scholar 

  4. Le Maître OP, Knio OM Scientific Computation. In Spectral Methods for Uncertainty Quantification. Springer, New York; 2010.

    Chapter  Google Scholar 

  5. Fox BL: Strategies for Quasi-Monte Carlo. Springer, Berlin; 1999.

    Book  Google Scholar 

  6. Babuska I, Tempone R, Zouraris G: Galerkin finite element approximation of stochastic elliptic partial differential equations. SIAM J. Numer. Anal. 2004,42(2):800–825. 10.1137/S0036142902418680

    Article  MATH  MathSciNet  Google Scholar 

  7. Pulch R: Polynomial chaos for multirate partial differential algebraic equations with random parameters. Appl. Numer. Math. 2009,59(10):2610–2624. 10.1016/j.apnum.2009.05.015

    Article  MATH  MathSciNet  Google Scholar 

  8. Xiu D: Fast numerical methods for stochastic computations. Commun. Comput. Phys. 2008, 5: 242–272.

    MathSciNet  Google Scholar 

  9. Maître OPL, Knio OM, Najm HN, Ghanem RG: Uncertainty propagation using Wiener-Haar expansions. J. Comput. Phys. 2004, 197: 28–57. 10.1016/j.jcp.2003.11.033

    Article  MATH  MathSciNet  Google Scholar 

  10. Ghanem R, Spanos P: Stochastic Finite Elements: A Spectral Approach. Springer, New York; 1991.

    Book  MATH  Google Scholar 

  11. Xiu D, Karniadakis G: The Wiener-Askey polynomial chaos for stochastic differential equations. SIAM J. Sci. Comput. 2002,24(2):619–644. 10.1137/S1064827501387826

    Article  MATH  MathSciNet  Google Scholar 

  12. Schevenels M, Lombaerts G, Degrande G: Application of the stochastic finite element method for Gaussian and non-Gaussian systems. ISMA2004 International Conference on Noise and Vibration Engineering Leuven, Belgien 2004.

    Google Scholar 

  13. Kalman R: A new approach to linear filtering and prediction problems. Trans. ASME, J. Basic Eng. 1960, 82: 35–45.

    Article  Google Scholar 

  14. Welch G, Bishop G: An introduction to the Kalman filter. Tr 95–041, UNC-Chapel Hill; 2006.

    Google Scholar 

  15. Julier SJ, Uhlmann JK, Durrant-Whyte HF: A new method for the nonlinear transformation of means and covariances in filters and estimators. IEEE Trans. Autom. Control 2000,45(3):477–482. 10.1109/9.847726

    Article  MATH  MathSciNet  Google Scholar 

  16. Wiener N: The homogeneous chaos. Am. J. Math. 1938, 60: 897–936. 10.2307/2371268

    Article  MathSciNet  Google Scholar 

  17. Cameron R, Martin W: The orthogonal development of nonlinear functionals in series of Fourier-Hermite functionals. Ann. Math. 1947,48(2):385–392. 10.2307/1969178

    Article  MATH  MathSciNet  Google Scholar 

  18. Wan X, Karniadakis G: An adaptive multi-element generalized polynomial chaos method for stochastic differential equations. J. Comput. Phys. 2005, 209: 617–642. 10.1016/j.jcp.2005.03.023

    Article  MATH  MathSciNet  Google Scholar 

  19. Maître OPL, Njam HN, Ghanem RG, Knio OM: Multi-resolution analysis of Wiener-type uncertainty propagation schemes. J. Comput. Phys. 2004,197(2):502–531. 10.1016/j.jcp.2003.12.020

    Article  MATH  MathSciNet  Google Scholar 

  20. Wan X, Karniadakis G: Multi-element generalized polynomial chaos for arbitrary probability measures. SIAM J. Sci. Comput. 2006,26(3):901–928.

    Article  MathSciNet  Google Scholar 

  21. Augustin F, Rentrop P: Stochastic Galerkin techniques for random ordinary differential equations. Numer. Math. 2012, 122: 399–419. 10.1007/s00211-012-0466-8

    Article  MATH  MathSciNet  Google Scholar 

  22. Paffrath M, Wever U: Stochastic integration methods and their application to reliability analysis. ASME Turbo Expo Copenhagen, Denmark 2012.

    Google Scholar 

  23. Villegas M: Random delay differential equations: application to biofilm modeling. Phd thesis. Techn. Universität München (Fachbereich Mathematik, Prof. Rentrop); 2011.

    Google Scholar 

  24. Harris DO: Probabilistic crack growth and modeling. In Reliability-Based Mechanical Design. Edited by: Cruse T. Marcel Dekker, New York; 1997.

    Google Scholar 

  25. Forman RG: Fatigue crack growth computer program, NASA/FLAGRO version 2.0. Johnson Space Center Report JSC-22267A, NASA; 1993.

    Google Scholar 

  26. Juan Z, Zhongua Z, Yaohong S: Hopf bifurcation of a bacteria-immunity system with delayed quorum sensing. Appl. Math. Comput. 2010, 215: 3936–3949. 10.1016/j.amc.2009.11.042

    Article  MATH  MathSciNet  Google Scholar 

  27. Kuhlmann H: Kalman-filtering with coloured measurement noise for deformation analysis. Proceedings, 11th FIG Symposium on Deformation Measurements Santorini, Greece 2003.

    Google Scholar 

  28. Weber B: Numerik nichtlinearer Systeme mit Unsicherheiten. Master thesis. Technische Universität München, Prof. Rentrop; 2008.

    Google Scholar 

  29. Julier S, Uhlmann J: Unscented filtering and nonlinear estimation. Proc. IEEE 2004,92(3):401–422. 10.1109/JPROC.2003.823141

    Article  Google Scholar 

  30. Janson S: Gaussian Hilbert Spaces. Cambridge University Press, Cambridge; 1997.

    Book  MATH  Google Scholar 

  31. Ernst OG, Mugler A, Starkloff HJ, Ullmann E: On the convergence of generalized polynomial chaos expansions. ESAIM Math. Model. Numer. Anal. 2012, 46: 317–339. 10.1051/m2an/2011045

    Article  MATH  MathSciNet  Google Scholar 

  32. Guglielmi N, Hairer E: Users’ guide for the code RADAR5 - version 2.1. Technical Report - available at http://www.unige.ch/~hairer/software.html 2005.

    Google Scholar 

  33. Hairer E, Wanner G Stiff and Differential-Algebraic Problems. In Solving Ordinary Differential Equations II. Springer, Berlin; 1996.

    Chapter  Google Scholar 

  34. Augustin F: Generalized Wiener Expansions for the Numerical Solution of Random Ordinary Differential Equations. Dr. Hut, München; 2012.

    Google Scholar 

  35. Hairer E, Nørsett SP, Wanner G: Solving Ordinary Differential Equations I. Springer, Berlin; 1993.

    MATH  Google Scholar 

  36. Villegas M: Biofilms with quorum sensing: A numerical approach for delay differential equations with uncertainties. Berichte aus den numerischen Arbeitsgruppen. Technische Universität München; 2011, 8/2011.

    Google Scholar 

Download references

Author information

Authors and Affiliations

Authors

Corresponding authors

Correspondence to Florian Augustin, Meinhard Paffrath, Manuel Villegas or Utz Wever.

Authors’ original submitted files for images

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Augustin, F., Gilg, A., Paffrath, M. et al. An accuracy comparison of polynomial chaos type methods for the propagation of uncertainties. J.Math.Industry 3, 2 (2013). https://doi.org/10.1186/2190-5983-3-2

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/2190-5983-3-2

Keywords