# A Proposed Stochastic Finite Difference Approach Based on Homogenous Chaos Expansion

Authors: O. H. Galal
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

This paper proposes a stochastic finite difference approach, based on homogenous chaos expansion (SFDHC). The said approach can handle time dependent nonlinear as well as linear systems with deterministic or stochastic initial and boundary conditions. In this approach, included stochastic parameters are modeled as second-order stochastic processes and are expanded using Karhunen-Loève expansion, while the response function is approximated using homogenous chaos expansion. Galerkin projection is used in converting the original stochastic partial differential equation (PDE) into a set of coupled deterministic partial differential equations and then solved using finite difference method. Two well-known equations were used for efficiency validation of the method proposed. First one being the linear diffusion equation with stochastic parameter and the second is the nonlinear Burger's equation with stochastic parameter and stochastic initial and boundary conditions. In both of these examples, the probability distribution function of the response manifested close conformity to the results obtained from Monte Carlo simulation with optimized computational cost.


---

## Body

## 1. Introduction

Stochastic processes are prevalent in nature. They affect all physical phenomena both from external and internal sources. Stochastic excitation and stochastic boundary conditions are a couple of examples of these external sources, while internal inherent random material heterogeneity also exists. To resolve these sources of stochastic inputs, stochastic differential equation is used for mathematical modeling. Methods of solution to these models are classified into two, statistical and nonstatistical. Monte Carlo simulation (MCS) is one of the most popular among the statistical methods of solution in use. As for the nonstatistical methodologies, these include stochastic perturbation method (SPM), stochastic finite element method (SFEM), and stochastic difference method (SFDM), the latter of which is the focus of this paper.The nonstatistical methods’ expansion process involves the discretization of input random fields. A projection scheme is further used resulting in deterministic equations solved by merely using one of the deterministic methods of solutions. Among these nonstatistical methods, SFEM is attracting the most attention. It has been applied to solve a wide range of problems especially in the field of structure mechanics (see, e.g., [1, 2]). In contrast, FDM is rarely studied in the field of stochastic analysis. These rare studies include Kaminski [3] which introduced a second-order perturbation, second probabilistic moment analysis in the context of FDM. This technique uses the perturbation method to expand the input random fields and was applied to time dependent linear problems with a small range of variability, where the first two statistical moments of the response were calculated. In Kamiński’s later works, he introduced a generalized version of this technique [4] which calculated higher probabilistic moments for any random dispersion of input variables. Ma et al. [5] presented a numerical method for a class of forward backward stochastic differential equations (FBSDEs). Their method is based on the four step schemeusing a Hermite-spectral method to approximate the solution to the decoupling quasilinear PDE on the whole space. Hovanessian and  Chang  [6] solved a second-order linear differential equation representing a series RLC network with step function, sinusoidal and stochastic inputs. Soheili et al. [7] later introduced two explicit finite difference schemes in order to approximate the solution of stochastic partial differential equations of Ito type.This paper utilizes the homogenous chaos expansion in the context of finite difference method (SFDHC). In this technique, the random inputs are discretized using Karhunen-Loève (KL) expansion, while the response is represented in terms of homogenous chaos (HC) expansion then a Galerkin projection scheme is applied to generate a system of deterministic equations, which can then be solved using classical deterministic FDM. Enriching other popular methods like SFEM, the proposed technique can handle time dependent nonlinear as well as linear problems with stochastic or deterministic initial and boundary conditions using less computational cost. Moreover, results are obtained in the form of a probability distribution function (pdf) of the response at the grid points where the extreme values and statistical moments can be easily evaluated, an obvious improvement over the usual solution representation of merelymean values and standard deviations. Two illustrative problems are used to validate the efficiency of the proposed technique, the first one being linear diffusion equation with stochastic parameter and the second is the nonlinear Burger’s equation with stochastic parameter and stochastic initial and boundary conditions. Both of these illustrative problem results are compared with MCS results. The comparison showed close conformity between SFDHC and MCS with optimized computational cost in favor of SFDHC.


## 2.  Random Field Discretization

When one or more of the physical parameters are modelled as random fields, it is essential, first, to represent them by an enumerable set of random variables. Quite a number of random field discretization techniques were cited in the literature [1, 8]. However, during the past two decades, KL and HC expansions have risen as two of the most practical techniques.


### 2.1. Karhunen-Loève Expansion

KL expansion is dependent on the covariance kernel [1]. Let α(x;θ) be a random process and Cαα(x1,x2) be its correlation function. The discretized version of this field is written as
(1)α(x;θ)=α¯(x)+∑i=1∞λifi(x)ξi(θ),
where α¯(x) denotes the mean value of α(x;θ), {ξi(θ)}i=1∞ is a set of uncorrelated random variables, and λi and fi(x) are the eigen values and the eigen functions, respectively. The said eigen value and eigen function are evaluated by solving Fredholm integral equation:
(2)∫DCαα(x1,x2)fi(x1)dx1=λifi(x2),
where D is the spatial domain over which α(x;θ) is defined, and x1, x2∈D. In the computational implementation of the expansion, (1) is truncated at the Mth term to reach the finite-dimensional approximation:
(3)α(x;θ)=α¯(x)+∑i=1Mλifi(x)ξi(θ).


### 2.2. Homogeneous Chaos Expansion

The solution process for stochastic differential equations is notknown a priori, which means that its covariance function is unknown. KL expansion is, in this case, inapplicable requiring HC expansion instead. HC expansion considers a probability space (Ω,κ,P), where Ω denotes the space of basic outcomes, κ is the σ-algebra associated with Ω, and P is a probability measure on κ. Considering some physical domain D⊂Rd×T, (d=1,2,or3), this can be a combination of spatial and temporal dimensions. A stochastic process (SP) can be seen as a scalar or vector-valued function u(x,t,θ):D×Ω→Rb, where x is an element of the physical space, t denotes the time, and θ is a point in the sample space Ω. Because of the infinite-dimensional nature of the probability space, it has to be discretized by using a finite number of random variables {ξi(θ)}i=1n, n∈N.HC expansion polynomials are functions in the infinite set of random variables{ξi(θ)}i=1i=∞ and therefore of infinite dimensions. A finite subset {ξλi(θ)}i=1i=n was chosen to build these polynomials; hence, the finite-dimensional polynomials of any order p will be a subset of the infinite-dimensional  polynomials of the same order. The convergence of the finite-dimensional polynomials is related to n as well as to the choice of the subset {ξλi(θ)}i=1i=n; the choice of this subset can be taken according to K-L expansion of the random process to ensure the convergence of a polynomial chaos (PC) expansion in the mean-square sense. As shown in [1, 9, 10], the solution process can be expressed as a summation of nonlinear functionals of the set {ξi(θ)}i=1∞ multiplied by deterministic constants. These functionals can be expanded as a set of polynomials of second-order Gaussian random variables, ξi(θ)∈N(0,1), referred to as polynomial chaos. Hence, the random process, which represents the random field, can be expressed as
(4)u(θ)=a0Γ0+∑i1=1∞ai1Γ1(ξi1(θ))+∑i1=1∞∑i2=1∞ai1i2Γ2(ξi1(θ),ξi2(θ))+⋯;
in which Γn(ξi1(θ),…,ξin(θ)) is the polynomial chaos of order n in the set of variables (ξi1(θ),…,ξin(θ)) given in [10]. Truncating at the Pth order, (4) can be then written as
(5)u(x,t,ξM(θ))=∑i=0Pu-i(x,t)ψi[{ξM(θ)}]
or simply
(6)u(x,t,ξ)=∑i=0Pu-iψi.The total number of these polynomials is(P+1), where
(7)P+1=(M+p)!p!(M)!.In (5), u-i(x,t) is a set of deterministic coefficients dependent on both x and t, ψi[{ξM}] is a set polynomials (up to order p) in the set of random variables (ξi1(θ),…,ξiM(θ)), and M refers to the number of terms in KL expansion. For example, if M=p=2, the set of second-order-second-degree PC consists of 6 elements; {1,ξ1,ξ2,ξ12-1,ξ1ξ2,ξ22-1}. Also, if M=4, p=2, the set of second-order-fourth-degree PC consists of 15 elements; {1,ξ1,ξ2,ξ3,ξ4,ξ12-1,ξ1ξ2,ξ1ξ3,ξ1ξ4,ξ22-1,ξ2ξ3,ξ2ξ4,ξ32-1,ξ3ξ4,ξ42-1}. Also, ψi[{ξn}] have the following orthogonality properties:
(8)〈ψi[{ξn}]〉=1fori=0andzerootherwise,〈ψi[{ξn}]ψj[{ξn}]〉=δij〈ψi2[{ξn}]〉.According to the Cameron-Martin theorem [11], for a fixed value of x and t, this expansion converges to any L2(Ω) functional in the L2(Ω) sense. This implies that the application of polynomial chaos is restricted to those stochastic processes with finite second-order moments (finite variance), and this applies to most physical processes yielding
(9)∫θ∈Ω|u(x,t,θ)|2dP(θ)<∞.


## 2.1. Karhunen-Loève Expansion

KL expansion is dependent on the covariance kernel [1]. Let α(x;θ) be a random process and Cαα(x1,x2) be its correlation function. The discretized version of this field is written as
(1)α(x;θ)=α¯(x)+∑i=1∞λifi(x)ξi(θ),
where α¯(x) denotes the mean value of α(x;θ), {ξi(θ)}i=1∞ is a set of uncorrelated random variables, and λi and fi(x) are the eigen values and the eigen functions, respectively. The said eigen value and eigen function are evaluated by solving Fredholm integral equation:
(2)∫DCαα(x1,x2)fi(x1)dx1=λifi(x2),
where D is the spatial domain over which α(x;θ) is defined, and x1, x2∈D. In the computational implementation of the expansion, (1) is truncated at the Mth term to reach the finite-dimensional approximation:
(3)α(x;θ)=α¯(x)+∑i=1Mλifi(x)ξi(θ).


## 2.2. Homogeneous Chaos Expansion

The solution process for stochastic differential equations is notknown a priori, which means that its covariance function is unknown. KL expansion is, in this case, inapplicable requiring HC expansion instead. HC expansion considers a probability space (Ω,κ,P), where Ω denotes the space of basic outcomes, κ is the σ-algebra associated with Ω, and P is a probability measure on κ. Considering some physical domain D⊂Rd×T, (d=1,2,or3), this can be a combination of spatial and temporal dimensions. A stochastic process (SP) can be seen as a scalar or vector-valued function u(x,t,θ):D×Ω→Rb, where x is an element of the physical space, t denotes the time, and θ is a point in the sample space Ω. Because of the infinite-dimensional nature of the probability space, it has to be discretized by using a finite number of random variables {ξi(θ)}i=1n, n∈N.HC expansion polynomials are functions in the infinite set of random variables{ξi(θ)}i=1i=∞ and therefore of infinite dimensions. A finite subset {ξλi(θ)}i=1i=n was chosen to build these polynomials; hence, the finite-dimensional polynomials of any order p will be a subset of the infinite-dimensional  polynomials of the same order. The convergence of the finite-dimensional polynomials is related to n as well as to the choice of the subset {ξλi(θ)}i=1i=n; the choice of this subset can be taken according to K-L expansion of the random process to ensure the convergence of a polynomial chaos (PC) expansion in the mean-square sense. As shown in [1, 9, 10], the solution process can be expressed as a summation of nonlinear functionals of the set {ξi(θ)}i=1∞ multiplied by deterministic constants. These functionals can be expanded as a set of polynomials of second-order Gaussian random variables, ξi(θ)∈N(0,1), referred to as polynomial chaos. Hence, the random process, which represents the random field, can be expressed as
(4)u(θ)=a0Γ0+∑i1=1∞ai1Γ1(ξi1(θ))+∑i1=1∞∑i2=1∞ai1i2Γ2(ξi1(θ),ξi2(θ))+⋯;
in which Γn(ξi1(θ),…,ξin(θ)) is the polynomial chaos of order n in the set of variables (ξi1(θ),…,ξin(θ)) given in [10]. Truncating at the Pth order, (4) can be then written as
(5)u(x,t,ξM(θ))=∑i=0Pu-i(x,t)ψi[{ξM(θ)}]
or simply
(6)u(x,t,ξ)=∑i=0Pu-iψi.The total number of these polynomials is(P+1), where
(7)P+1=(M+p)!p!(M)!.In (5), u-i(x,t) is a set of deterministic coefficients dependent on both x and t, ψi[{ξM}] is a set polynomials (up to order p) in the set of random variables (ξi1(θ),…,ξiM(θ)), and M refers to the number of terms in KL expansion. For example, if M=p=2, the set of second-order-second-degree PC consists of 6 elements; {1,ξ1,ξ2,ξ12-1,ξ1ξ2,ξ22-1}. Also, if M=4, p=2, the set of second-order-fourth-degree PC consists of 15 elements; {1,ξ1,ξ2,ξ3,ξ4,ξ12-1,ξ1ξ2,ξ1ξ3,ξ1ξ4,ξ22-1,ξ2ξ3,ξ2ξ4,ξ32-1,ξ3ξ4,ξ42-1}. Also, ψi[{ξn}] have the following orthogonality properties:
(8)〈ψi[{ξn}]〉=1fori=0andzerootherwise,〈ψi[{ξn}]ψj[{ξn}]〉=δij〈ψi2[{ξn}]〉.According to the Cameron-Martin theorem [11], for a fixed value of x and t, this expansion converges to any L2(Ω) functional in the L2(Ω) sense. This implies that the application of polynomial chaos is restricted to those stochastic processes with finite second-order moments (finite variance), and this applies to most physical processes yielding
(9)∫θ∈Ω|u(x,t,θ)|2dP(θ)<∞.


## 3.  A Proposed Stochastic Finite Difference Approach Based on Homogenous Chaos Expansion (SDFHC)

In this proposed technique, the involved random parameters in SDE are expanded using KL expansion, which is truncated at theMth term. The response function u is, in turn, expanded using HC expansion since its covariance kernel is not known a priori. Applying Galerkin projection scheme through multiplying both sides of the resulting equation by ψj and taking the statistical average yields a deterministic system of equations in u-i(x,t) which can be solved simply using the classical deterministic FDM. Consider
(10)Λ(x,t,ξ(θ);u)=f(x,t,ξ(θ)),
where f(x,t,ξ(θ)) is a random external excitation and Λ is a differential operator involving differentiation in space and/or time which can be nonlinear and defined on the product space D×Ω. Also, u satisfies the initial and boundary conditions, deterministic or stochastic, and they are defined on some section of the boundary ∂D1∈∂D. Using (6) to expand u in (10) yields
(11)Λ(x,t,ξ(θ);∑i=0Pu-iψi)=f(x,t,ξ(θ)).Multiply the different orthogonal polynomials of the finite expansion in (11) and get statistical average results to
(12)〈Λ(x,t,ξ;∑i=0Pu-iψi),ψj〉=〈f(x,t,ξ),ψj〉,j=0,1,2,…,P.From (12), it should be noted that (P+1) set of coupled equations u-i(x,t) is generated. This most recent step is actually Galerkin projection in random space. The projection ensures the orthogonality of the residual to the functional space spanned by the finite-dimensional basis {ψi}i=0i=P. As predicted, a strictly deterministic u-i(x,t) was achieved which would merely require a classical FDM-based solution in space and time. The resulting feature is mainly due to the averaging inherent in the projection, eliminating the randomness of the system model. It could also be noted that no additional errors, except for truncation errors, manifest from this step.For a numerical solution using FDM, the physical domainD is covered by a grid of mesh points (x,t)=(mΔx,nΔt) where Δx and Δt are mesh parameters that are relatively infinitesimal and m and n are positive integers. The approximate solution of u-i(x,t) at these mesh points will be denoted by u-i(m,n); provided that the derivatives of u-i(x,t) have been replaced by finite difference quotients. Moreover, the choice of both Δx and Δt must satisfy all necessary conditions of the deterministic FD scheme like convergence, stability, and consistency. Once u-i(x,t) is obtained at every point in the mesh, (6) can be used to get the solution vector where the mean is given as
(13)〈u(x,t,ξ)〉=u-0(x,t),
and the covariance matrix for two points u and v is(14a)Cov(u,v)=〈(u-〈u〉)(v-〈v〉)T〉,
and the resulting variance vector of each point is
(14b)Var(u(x,t,ξ))=∑i=1P〈ψi2〉u-i2(x,t).


## 4. Illustrative Example I: The Stochastic Diffusion Equation

In order to validate the efficiency of the proposed SFDHC technique, stochastic diffusion equation is evaluated under the effect of random parameterγ(x,θ). The results are then compared with that obtained from MCS method.


### 4.1. Problem Formulation

Consider the diffusion equation [12] with random parameter γ(x,θ) in the form
(15)∂u(x,t,ξ)∂t=γ(x,θ)∂2u(x,t,ξ)∂x2,
where: x∈[0,L] and t≥0. Also, consider the case of deterministic initial and boundary conditions
(16)u(0,t)=q1(t),u(L,t)=q2(t),u(x,0)=h(x).Moreover, the random parameter can be written as:(17)γ(x,θ)=γ(1+α(x;θ)),
where α(x;θ) is SP which expresses the spatial randomness in the parameter γ(x,θ). Expanding both α(x;θ) and u(x,t,ξ) using (3) and (6), respectively, yields
(18)∑i=0P∂u-i∂tψi=γ(1+∑i=1Mλjfj(x)ξj)∑i=0P∂2u-i∂x2ψi.Using Galerkin projection scheme and multiplying both sides of (18) by ψj then taking the statistical average yields
(19)∑i=0P∑j=0Pbij∂u-j∂t=γ∑i=0P∑j=0Pbij∂2u-j∂x2+γ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)∂2u-j∂x2,
where
(20)bij=〈ψiψj〉,dijk=〈ψiψjξk〉.Replacing the derivatives ofu-i by their difference quotients at the mesh points leads to
(21)1Δt∑i=0P∑j=0Pbij(u-j(m,n+1)-u-j(m,n))=γ(Δx)2∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))-γ(Δx)2∑i=0P∑j=0P∑k=1Mdijkλkfk(x)2×(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n)).Introducingμ=γΔt(Δx)-2 and solving the last equation for u-i(m,n+1) yields
(22)∑i=0P∑j=0Pbiju-j(m,n+1)=∑i=0P∑j=0Pbiju-j(m,n)+γμ∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+γμ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)2×(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n)),
which can be solved using an explicit scheme to find u-i(x,t) at the different mesh points.


### 4.2. Numerical Implementation and Results

To implement the results obtained in (22), consider α(x;θ) as a second-order Gaussian SP with exponential covariance kernel given in [1]
(23)Cov(x1,x2)=σ2exp(-|x1-x2|lx),
where σ denotes the coefficient of variation (COV) of the random field (which is constant) and lx is the correlation length. Assume σ2=0.01 and both lx and γ are unity. Also, let the initial and boundary conditions be considered as
(24)q1(t)=q2(t)=0,h(x)=sinπxL,
where L=1, x∈[0,1], and t∈[0,0.04]. Consider also the mesh parameters Δx=0.1 and Δt=0.004 so that μ≤1/2, which is a necessary condition for the convergence of the FD scheme. This, produces a moderately course mesh of 11 joints in both x and t directions, with a total number of unknown joints NP=90 (without the boundary ones). The problem is solved using MCS [13], where α(x;θ) is expanded using KL expansion with M=6. A large sample size of 1×106 is then used to generate the random variables where exact deterministic solution, given in [12], is evaluated for each input. This in turn, becomes a benchmark for the proposed technique. HC of second and fourth degrees, SFD2 and SFD4, respectively, with order two were employed to derive the SFDHC’s solution. The same sample size was then used to generate the HC polynomials. The pdf of the solution process u is plotted using MCS on SFD2 and SFD4 for 2 arbitrary points, point A on the midspan at (0.5,0.024) and B near the boundary at (0.9,0.036), as shown in Figures 1 and 2, respectively. The mean and standard deviation (SD) of u are also plotted at an arbitrary time level (t=0.04) as shown in Figures 3 and 4.Figure 1
The pdf at pointA using SFDHC and MCS.Figure 2
The pdf at pointB using SFDHC and MCS.Figure 3
The mean att=0.04 using SFDHC and MCS.Figure 4
The SD att=0.04 using SFDHC and MCS.


### 4.3. Error Analysis

The resulting response from the SFDHC technique, with two degrees of the polynomial chaos (SFD2 and SFD4), and MCS results (incorporating the exact solution given in [12]) were compared in terms of the maximum and minimum errors, their locations, and average error. From this, the relative accuracy of SFDHC technique is then recorded as shown in Table 1.Table 1
Errors and accuracy for the stochastic diffusion equation.
The mean
The standard deviation
Min. error/location
Max. error/location
Average error/accuracy
Min. error/location
Max. error/location
Average error/accuracy
SFD2
0.059%
0.696%
0.317%
0.1412%
11.09%
2.841%
(
0.1,0.004
),(0.9,0.004)
(
0.5,0.04
)
99.683%
(
0.1,0.04
),(0.9,0.04)
(
0.5,0.04
)
97.159%
SFD4
0.0035%
0.221%
0.128%
0.0024%
6.3645%
1.393%
(
0.1,0.004
),(0.9,0.004)
(
0.5,0.04
)
99.872%
(
0.3,0.024
),(0.7,0.024)
(
0.5,0.04
)
98.607%Note that the average error for the mean value and SD was calculated through the summation of the errors over all points and then dividing the sum by the number of points NP. Table1 shows good convergence of the SFDHC technique with increasing the degree of the used polynomials. Evidently, for SFD2, the mean value and SD errors were maximum at the midspan, while they were minimum near the boundary. For SFD4, however, the mean value error has the same positions as SFD2, while the SD error is minimum around the first and third quarters. This conclusion is made from this particular problem with their corresponding inputs. Generally, from the insignificant errors recorded, the SFDHC technique exhibits very high accuracy for calculating the mean value and SD.


## 4.1. Problem Formulation

Consider the diffusion equation [12] with random parameter γ(x,θ) in the form
(15)∂u(x,t,ξ)∂t=γ(x,θ)∂2u(x,t,ξ)∂x2,
where: x∈[0,L] and t≥0. Also, consider the case of deterministic initial and boundary conditions
(16)u(0,t)=q1(t),u(L,t)=q2(t),u(x,0)=h(x).Moreover, the random parameter can be written as:(17)γ(x,θ)=γ(1+α(x;θ)),
where α(x;θ) is SP which expresses the spatial randomness in the parameter γ(x,θ). Expanding both α(x;θ) and u(x,t,ξ) using (3) and (6), respectively, yields
(18)∑i=0P∂u-i∂tψi=γ(1+∑i=1Mλjfj(x)ξj)∑i=0P∂2u-i∂x2ψi.Using Galerkin projection scheme and multiplying both sides of (18) by ψj then taking the statistical average yields
(19)∑i=0P∑j=0Pbij∂u-j∂t=γ∑i=0P∑j=0Pbij∂2u-j∂x2+γ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)∂2u-j∂x2,
where
(20)bij=〈ψiψj〉,dijk=〈ψiψjξk〉.Replacing the derivatives ofu-i by their difference quotients at the mesh points leads to
(21)1Δt∑i=0P∑j=0Pbij(u-j(m,n+1)-u-j(m,n))=γ(Δx)2∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))-γ(Δx)2∑i=0P∑j=0P∑k=1Mdijkλkfk(x)2×(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n)).Introducingμ=γΔt(Δx)-2 and solving the last equation for u-i(m,n+1) yields
(22)∑i=0P∑j=0Pbiju-j(m,n+1)=∑i=0P∑j=0Pbiju-j(m,n)+γμ∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+γμ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)2×(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n)),
which can be solved using an explicit scheme to find u-i(x,t) at the different mesh points.


## 4.2. Numerical Implementation and Results

To implement the results obtained in (22), consider α(x;θ) as a second-order Gaussian SP with exponential covariance kernel given in [1]
(23)Cov(x1,x2)=σ2exp(-|x1-x2|lx),
where σ denotes the coefficient of variation (COV) of the random field (which is constant) and lx is the correlation length. Assume σ2=0.01 and both lx and γ are unity. Also, let the initial and boundary conditions be considered as
(24)q1(t)=q2(t)=0,h(x)=sinπxL,
where L=1, x∈[0,1], and t∈[0,0.04]. Consider also the mesh parameters Δx=0.1 and Δt=0.004 so that μ≤1/2, which is a necessary condition for the convergence of the FD scheme. This, produces a moderately course mesh of 11 joints in both x and t directions, with a total number of unknown joints NP=90 (without the boundary ones). The problem is solved using MCS [13], where α(x;θ) is expanded using KL expansion with M=6. A large sample size of 1×106 is then used to generate the random variables where exact deterministic solution, given in [12], is evaluated for each input. This in turn, becomes a benchmark for the proposed technique. HC of second and fourth degrees, SFD2 and SFD4, respectively, with order two were employed to derive the SFDHC’s solution. The same sample size was then used to generate the HC polynomials. The pdf of the solution process u is plotted using MCS on SFD2 and SFD4 for 2 arbitrary points, point A on the midspan at (0.5,0.024) and B near the boundary at (0.9,0.036), as shown in Figures 1 and 2, respectively. The mean and standard deviation (SD) of u are also plotted at an arbitrary time level (t=0.04) as shown in Figures 3 and 4.Figure 1
The pdf at pointA using SFDHC and MCS.Figure 2
The pdf at pointB using SFDHC and MCS.Figure 3
The mean att=0.04 using SFDHC and MCS.Figure 4
The SD att=0.04 using SFDHC and MCS.


## 4.3. Error Analysis

The resulting response from the SFDHC technique, with two degrees of the polynomial chaos (SFD2 and SFD4), and MCS results (incorporating the exact solution given in [12]) were compared in terms of the maximum and minimum errors, their locations, and average error. From this, the relative accuracy of SFDHC technique is then recorded as shown in Table 1.Table 1
Errors and accuracy for the stochastic diffusion equation.
The mean
The standard deviation
Min. error/location
Max. error/location
Average error/accuracy
Min. error/location
Max. error/location
Average error/accuracy
SFD2
0.059%
0.696%
0.317%
0.1412%
11.09%
2.841%
(
0.1,0.004
),(0.9,0.004)
(
0.5,0.04
)
99.683%
(
0.1,0.04
),(0.9,0.04)
(
0.5,0.04
)
97.159%
SFD4
0.0035%
0.221%
0.128%
0.0024%
6.3645%
1.393%
(
0.1,0.004
),(0.9,0.004)
(
0.5,0.04
)
99.872%
(
0.3,0.024
),(0.7,0.024)
(
0.5,0.04
)
98.607%Note that the average error for the mean value and SD was calculated through the summation of the errors over all points and then dividing the sum by the number of points NP. Table1 shows good convergence of the SFDHC technique with increasing the degree of the used polynomials. Evidently, for SFD2, the mean value and SD errors were maximum at the midspan, while they were minimum near the boundary. For SFD4, however, the mean value error has the same positions as SFD2, while the SD error is minimum around the first and third quarters. This conclusion is made from this particular problem with their corresponding inputs. Generally, from the insignificant errors recorded, the SFDHC technique exhibits very high accuracy for calculating the mean value and SD.


## 5. Illustrative Example II: Stochastic Nonlinear Burger’s Equation

To illustrate the proposed SFDHC technique on a nonlinear case, consider stochastic nonlinear Burger’s equation (SBE) [14]. In addition, let us involve random parameters and expose them to deterministic or stochastic random excitation with deterministic or stochastic boundary and initial conditions.


### 5.1. Problem Formulation

With a random parameterγ(x,θ) and random excitation p(x,t,θ), the nonlinear stochastic Burger’s equation can then be written as
(25)∂u(x,t,ξ)∂t+u(x,t,ξ)∂u(x,t,ξ)∂x-γ(x,θ)∂2u(x,t,ξ)∂x2=p(x,t,θ).Given that the domain and initial and boundary conditions are consistent with Section4.2 and the excitation functions can be written as
(26)p(x,t,θ)=p(x,t)(1+β(x;θ)),
where β(x;θ) is an SP expressing the spatial randomness in the excitation function. Also, β(x;θ) can be expanded using (3) with γj and gj(x) as the eigen values and the eigen functions, associated with its covariance kernel, respectively. Expanding u(x,t,ξ) as in (6) yields
(27)∑i=0P∂u-i∂tψi+∑i=0Pu-iψi∑j=0P∂u-j∂xψj-γ(1+∑i=1Mλjfj(x)ξj)×∑i=0P∂2u-i∂x2ψi=p(x,t)(1+∑i=1Mγigi(x)ξi).Multiplying both sides of (17) by ψk and taking the statistical average yields
(28)∑i=0P∑j=0Pbij∂u-j∂t+∑i=0P∑j=0P∑k=0Pcijku-j∂u-k∂x-γ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)∂2u-j∂x2-γ∑i=0P∑j=0Pbij∂2u-j∂x2=p(x,t)+p(x,t)∑i=0P∑j=1Meijγjgj(x),
where
(29)cijk=〈ψiψjψk〉,eij=〈ψiξj〉.Replacing the derivatives ofu-i by their difference quotients and solving the last equation for u-i(m,n+1) yields
(30)∑i=0P∑j=0Pbiju-j(m,n+1)=∑i=0P∑j=0Pbiju-j(m,n)-12Δxμ∑i=0P∑j=0P∑k=0Pcijku-j(m,n)(u-k(m+1,n)-u-k(m-1,n))+γμ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)2×(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+γμ∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+Δtp(x,t)+Δtp(x,t)∑i=0P∑i=jMeijγjgj(x).


### 5.2. Numerical Implementation and Results

The proposed technique will be applied to SBE in two cases: Case1 using random parameter and, Case 2 using random initial and boundary conditions. Both Cases 1 and 2 are treated as follows.Case 1.
The same data used in Section4.2 concerning the domain, the initial and boundary conditions, the mesh parameters, and α(x;θ) will be employed here. First, assume the COV as unity for error analysis and for Figures 5 and 6. Then, some practical values (0.05, 0.10, and 0.20) will be considered in Figures 7 and 8. Both of SFD2 and SFD4 are used to solve (30), where the following results are obtained for points C(0.5,0.036) and D(0.9,0.036).Figure 5
The pdf at pointC using SFDHC and MCS.Figure 6
The pdf at pointD using SFDHC and MCS.Figure 7
The pdf at pointC for different values of COV.Figure 8
The pdf at pointD for different values of COV.


### 5.3. Error Analysis

As in Section4.3, errors and accuracy are recorded as shown in Table 2.Table 2
Errors and accuracy for the stochastic Burger's equation.
The mean
The standard deviation
Min. error/location
Max. error/location
Average error/accuracy
Min. error/location
Max. error/location
Average error/accuracy
SFD2
0.059%
0.696%
0.908%
0.1412%
11.09%
4.632%
(
0.6,0.004
)
(
0.1,0.04
)
99.092%
(
0.9,0.04
)
(
0.9,0.004
)
95.368%
SFD4
0.0035%
0.221%
0.863%
0.0024%
6.3645%
2.973%
(
0.8,0.028
)
(
0.1,0.04
)
99.137%
(
0.9,0.04
)
(
0.4,0.04
)
97.027%Although the errors are slightly bigger than the linear case, the SFDHC technique still exhibits high performance for nonlinear problems.Case 2.
The proposed SFDHC technique can handle stochastic initial and boundary conditions expressed as(31)u(x,t,ξ)=∑i=0P1u-i(x,t)χi(ξ),
defined on some section of the boundary ∂D1∈∂D, where P1≤P and both degree and order of χi are less than or equal to those used for ψi. After executing the same procedure as in Case 1, it yields
(32)∑i=0P∑j=0Pbiju-j(m,n+1)=∑i=0P∑j=0Pbiju-j(m,n)-12Δxμ∑i=0P∑j=0P∑k=0Pcijku-j(m,n)(u-k(m+1,n)-u-k(m-1,n))+γμ∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+Δtp(x,t).
As a practical example, assume an initial condition with the form(33)u(x,0,ξ)=x(L-x)(1+α(x;θ))
and boundary conditions
(34)u(0,t,ξ)=α(0;θ),u(L,t,ξ)=α(L;θ),
where α(x;θ) is the same Gaussian SP used in Section 4.2 with COV = 1. Then, u-0(x,0)=x(L-x), u-i(x,0)=σλifi(x) for 1≤i≤M, and u-i(x,0)=0 for i>M. Also, u-0(0,t)=u-0(L,t)=0, u-i(0,t)=σλifi(0), u-i(L,t)=σλifi(L) for 1≤i≤M and u-i(0,t)=u-i(L,t)=0 for i>M. Both of SFD4 and SFD2 are used to solve (32) with initial and boundary conditions given by (33) and (34), respectively. An MCS solution was provided through expanding the included stochastic processes using KL expansion then a sample of size 50,000 random numbers is used to generate the initial and boundary random values. Consequently, the problem is solved using FD forward scheme leading to 50,000 values of the stochastic response at each of the mesh points. A solution pdf for SFDHC is generated and then plotted for points C(0.5,0.036) and D(0.9,0.036) shown in Figures 9 and 10.Figure 9
The pdf at pointC using SFD4 and SFD2.Figure 10
The pdf at pointD using SFDHC and MCS.Comparing with MCS, SFD4 shows better approximation for the response (as expected) since it uses 4 random variables in its expansion, while SFD2 uses only 2.


## 5.1. Problem Formulation

With a random parameterγ(x,θ) and random excitation p(x,t,θ), the nonlinear stochastic Burger’s equation can then be written as
(25)∂u(x,t,ξ)∂t+u(x,t,ξ)∂u(x,t,ξ)∂x-γ(x,θ)∂2u(x,t,ξ)∂x2=p(x,t,θ).Given that the domain and initial and boundary conditions are consistent with Section4.2 and the excitation functions can be written as
(26)p(x,t,θ)=p(x,t)(1+β(x;θ)),
where β(x;θ) is an SP expressing the spatial randomness in the excitation function. Also, β(x;θ) can be expanded using (3) with γj and gj(x) as the eigen values and the eigen functions, associated with its covariance kernel, respectively. Expanding u(x,t,ξ) as in (6) yields
(27)∑i=0P∂u-i∂tψi+∑i=0Pu-iψi∑j=0P∂u-j∂xψj-γ(1+∑i=1Mλjfj(x)ξj)×∑i=0P∂2u-i∂x2ψi=p(x,t)(1+∑i=1Mγigi(x)ξi).Multiplying both sides of (17) by ψk and taking the statistical average yields
(28)∑i=0P∑j=0Pbij∂u-j∂t+∑i=0P∑j=0P∑k=0Pcijku-j∂u-k∂x-γ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)∂2u-j∂x2-γ∑i=0P∑j=0Pbij∂2u-j∂x2=p(x,t)+p(x,t)∑i=0P∑j=1Meijγjgj(x),
where
(29)cijk=〈ψiψjψk〉,eij=〈ψiξj〉.Replacing the derivatives ofu-i by their difference quotients and solving the last equation for u-i(m,n+1) yields
(30)∑i=0P∑j=0Pbiju-j(m,n+1)=∑i=0P∑j=0Pbiju-j(m,n)-12Δxμ∑i=0P∑j=0P∑k=0Pcijku-j(m,n)(u-k(m+1,n)-u-k(m-1,n))+γμ∑i=0P∑j=0P∑k=1Mdijkλkfk(x)2×(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+γμ∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+Δtp(x,t)+Δtp(x,t)∑i=0P∑i=jMeijγjgj(x).


## 5.2. Numerical Implementation and Results

The proposed technique will be applied to SBE in two cases: Case1 using random parameter and, Case 2 using random initial and boundary conditions. Both Cases 1 and 2 are treated as follows.Case 1.
The same data used in Section4.2 concerning the domain, the initial and boundary conditions, the mesh parameters, and α(x;θ) will be employed here. First, assume the COV as unity for error analysis and for Figures 5 and 6. Then, some practical values (0.05, 0.10, and 0.20) will be considered in Figures 7 and 8. Both of SFD2 and SFD4 are used to solve (30), where the following results are obtained for points C(0.5,0.036) and D(0.9,0.036).Figure 5
The pdf at pointC using SFDHC and MCS.Figure 6
The pdf at pointD using SFDHC and MCS.Figure 7
The pdf at pointC for different values of COV.Figure 8
The pdf at pointD for different values of COV.


## 5.3. Error Analysis

As in Section4.3, errors and accuracy are recorded as shown in Table 2.Table 2
Errors and accuracy for the stochastic Burger's equation.
The mean
The standard deviation
Min. error/location
Max. error/location
Average error/accuracy
Min. error/location
Max. error/location
Average error/accuracy
SFD2
0.059%
0.696%
0.908%
0.1412%
11.09%
4.632%
(
0.6,0.004
)
(
0.1,0.04
)
99.092%
(
0.9,0.04
)
(
0.9,0.004
)
95.368%
SFD4
0.0035%
0.221%
0.863%
0.0024%
6.3645%
2.973%
(
0.8,0.028
)
(
0.1,0.04
)
99.137%
(
0.9,0.04
)
(
0.4,0.04
)
97.027%Although the errors are slightly bigger than the linear case, the SFDHC technique still exhibits high performance for nonlinear problems.Case 2.
The proposed SFDHC technique can handle stochastic initial and boundary conditions expressed as(31)u(x,t,ξ)=∑i=0P1u-i(x,t)χi(ξ),
defined on some section of the boundary ∂D1∈∂D, where P1≤P and both degree and order of χi are less than or equal to those used for ψi. After executing the same procedure as in Case 1, it yields
(32)∑i=0P∑j=0Pbiju-j(m,n+1)=∑i=0P∑j=0Pbiju-j(m,n)-12Δxμ∑i=0P∑j=0P∑k=0Pcijku-j(m,n)(u-k(m+1,n)-u-k(m-1,n))+γμ∑i=0P∑j=0Pbij(u-j(m+1,n)-2u-j(m,n)+u-j(m-1,n))+Δtp(x,t).
As a practical example, assume an initial condition with the form(33)u(x,0,ξ)=x(L-x)(1+α(x;θ))
and boundary conditions
(34)u(0,t,ξ)=α(0;θ),u(L,t,ξ)=α(L;θ),
where α(x;θ) is the same Gaussian SP used in Section 4.2 with COV = 1. Then, u-0(x,0)=x(L-x), u-i(x,0)=σλifi(x) for 1≤i≤M, and u-i(x,0)=0 for i>M. Also, u-0(0,t)=u-0(L,t)=0, u-i(0,t)=σλifi(0), u-i(L,t)=σλifi(L) for 1≤i≤M and u-i(0,t)=u-i(L,t)=0 for i>M. Both of SFD4 and SFD2 are used to solve (32) with initial and boundary conditions given by (33) and (34), respectively. An MCS solution was provided through expanding the included stochastic processes using KL expansion then a sample of size 50,000 random numbers is used to generate the initial and boundary random values. Consequently, the problem is solved using FD forward scheme leading to 50,000 values of the stochastic response at each of the mesh points. A solution pdf for SFDHC is generated and then plotted for points C(0.5,0.036) and D(0.9,0.036) shown in Figures 9 and 10.Figure 9
The pdf at pointC using SFD4 and SFD2.Figure 10
The pdf at pointD using SFDHC and MCS.Comparing with MCS, SFD4 shows better approximation for the response (as expected) since it uses 4 random variables in its expansion, while SFD2 uses only 2.


## 6. Discussion and Conclusions

As evidenced in this work, the SFDHC technique was designed to solve a wide range of time dependent nonlinear as well as linear problems with stochastic or deterministic initial and boundary conditions and exposed to deterministic or stochastic excitation. This technique uses HC expansion in compliment with FDM, providing the capacity to present the analysis result as pdf and not just the statistical moments of the response. Using two illustrative examples, as representative of linear and nonlinear problems, the proposed technique has proven to be highly efficient and computationally optimized. Using SFD4, the efficiency of this technique exceeds 98.6% for both mean and SD evaluations in the linear problem and exceeds 97% for the nonlinear one. Moreover, the computation times spent for SFDHC, on a Core-i3 computer, were 3.08 and 3.91 seconds for linear and nonlinear problems, respectively, as against 40.19 and 989.82 seconds for MCS, respectively, which means that 7.66% and 0.395% only of the execution time are in favor of SFDHC. The overwhelming advantages manifested by the proposed technique open the door for the versatility manifested by FDM in the deterministic field into the stochastic field. FDM had shown its versatility in the deterministic field by directly dealing with time dependent and nonlinear systems through explicit or implicit schemes. These capabilities of FDM will be available in the stochastic field by the way of this proposed technique. The SFDHC technique can handle more complicated nonlinear problems in 2 and 3 dimensions with explicit or implicit schemes which are available in the FDM literature. However, the illustrative examples in this paper were one-dimensional problems and were solved explicitly. Further research concerning the range of variability of random inputs, the choice of mesh parameters which satisfy FDM’s essential conditions, and closed forms for the produced error is highly recommended and is deemed worthwhile.



---

*Source: 950469-2013-07-18.xml*