Numerical solution of time-dependent stochastic partial differential equations using RBF partition of unity collocation method based on finite difference

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Numerical solution of time-dependent stochastic partial differential equations using RBF partition of unity collocation method based on finite difference M. Esmaeilbeigi a,∗, O. Chatrabgoun b, M. Shafa a a b

Department of Mathematics, Faculty of Mathematical Sciences and Statistics, Malayer University, Malayer 65719-95863, Iran Department of Statistics, Faculty of Mathematical Sciences and Statistics, Malayer University, Malayer 65719-95863, Iran

a r t i c l e

i n f o

Keywords: Stochastic partial differential equations Meshless methods RBF collocation Partition of unity approximation

a b s t r a c t Meshfree methods based on radial basis functions (RBFs) are popular tools for the numerical solution of stochastic partial differential equations (SPDEs) due to their nice properties. However, the RBF collocation methods in global view have some disadvantages for the numerical solution of time-dependent SPDEs. Calculation of matrix condition number in the resulting dense linear systems indicates that the meshless method using global RBFs may be unstable at each realization to solve SPDEs. In order to avoid numerical instabilities in global RBF methods, we are interested in the use of RBF methods in local view for the numerical solution of time-dependent SPDEs. In this paper, the RBF partition of unity collocation method based on a finite difference scheme for the Gaussian random field (RBF-PU-FD) as a localized RBF approximation presented to deal with these issues. For this purpose, we simulate the Gaussian field with spatial covariance structure at a finite collection of predetermined collocation points. The matrices formed during the RBF-PU-FD method will be sparse and, hence, will not suffer from illconditioning and high computational cost. We will show that the method is viable through analyzing its numerical accuracy, CPU time, stability and sparsity structure. For the test problems, we perform 1000 realizations and statistical criterions such as mean, standard deviation, lower bound and upper bound of prediction are computed and evaluated using the Monte-Carlo method.

1. Introduction The study of stochastic partial differential equations (SPDEs) attracted a lot of interest due to the existence of different stochastic terms and uncertainty in many phenomena. This is the reason that SPDEs are able to fully capture the behavior of such phenomena. Models based on SPDEs play an important role in a range of application areas, including chemistry, physics, biology, geography, economics, finance, microelectronics, and nowadays also nanotechnology. The theoretical foundation of SPDEs and their applications are described in [11,19]. However, it is difficult to obtain the analytical solutions of SPDEs [1,11]. Thus, the numerical solutions of SPDEs become a fast-growing research area [9,26,40]. Several competitive methods are discussed for approximating SPDEs, including the finite difference (FD) method [27], the finite element (FE) method [29,31]. Some authors have used FDs and FEs for spatial variable discretization and then solved the resulting system of stochastic ordinary differential equations (SODEs) via Euler method or the Crank– Nicolson scheme [27,31].

Some of the existed problems in the numerical solution of SPDEs are solved by FD schemes, FE techniques, BE methods, but all these numerical methods are based on a grid discretization that has to improve the solution through adaptive meshing. It showed by Dehghan and Shirzadi that mesh generation is one of the biggest challenges in meshbased methods [16]. To overcome this difficulty, other tools proposed for the numerical solution of SPDEs [15,17]. In [22], Fasshauer and Ye showed that meshless methods are powerful numerical techniques that have been applied to solving SPDEs. And consequently, the purpose of meshless methods is to eliminate the structure of the mesh. The theoretical framework of the meshless approximation method for the numerical solution of SPDEs is presented in [41]. A class of meshless methods is based on radial basis function (RBF) methods which put radial functions as the basis functions for the collocation [12]. The main advantages of these methods lie in their simplicity and their effectiveness in dealing with high-dimensional problems with complicated geometries since no mesh generation is needed. Some authors like Cialenco et al. [12] and Ye [43] have investigated the numerical solution of time-dependent SPDEs. Furthermore, Fasshauer used an RBF method to approximate the solutions of time-

∗

Corresponding author. E-mail addresses: [email protected] (M. Esmaeilbeigi), [email protected] (O. Chatrabgoun), [email protected] (M. Shafa). https://doi.org/10.1016/j.enganabound.2019.03.013 Received 15 July 2018; Received in revised form 30 January 2019; Accepted 15 March 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Table 1 Examples of some popular RBFs, where (⋅)+ denotes the truncated power function, and 𝜈 ∈ ℕ.

independent SPDEs [23]. Also, RBFs applied for the numerical solution of nonlinear SPDEs [42]. However, the global RBF collocation methods have some disadvantages for the numerical solution of time-dependent SPDEs in [15,16]. Estimates of condition numbers of the matrices in the resulting dense linear systems indicate that the meshless method using RBFs in global view may be unstable at each realization to solve SPDEs. Also, the computational cost is the main barrier when using global RBF methods. In order to avoid numerical instabilities in global RBF methods, we are interested in the use of local RBF methods for the numerical solution of time-dependent SPDEs. For deterministic case, many authors followed this approach and applied local methods to solve different problems. In 2004, diffusion problem has been solved by Chantasiriwan using local RBF [8]. The collocation is made locally on the overlapping domains of influence and hence reduces the collocation matrix size. In the sequel, local RBF methods are used by researchers to solve different PDEs [13,32]. Shu et al. in 2003 gave a local approach to combine the meshfree nature of RBF and the high accuracy and simplicity of Differential Quadrature (DQ) method known as local RBF-DQ method [36]. In local RBF-DQ method, only those nodes are used, which are in the neighbourhood of the node under consideration, and are called supporting nodes. This technique has been used by researchers to solve PDEs [14]. Another promising approach is the RBF-PUM to solve PDEs which combines the partition of unity method with RBFs. The idea of RBFPUM method is to partition the domain into overlapping subdomains. The local approximation is done on the subdomains and combines to get the global approximation. RBF-PUM reduces the computational cost while maintaining high accuracy [7]. Apart from these methods there are other approaches to avoid the dense nature of the interpolation matrix such that these approaches are based on local approximation. In this paper, the RBF partition of unity collocation method based on a FD scheme (RBF-PU-FD) as a localized RBF approximation presented to deal with these issues. In fact, FD approximation leads to time derivative discretization for solving time-dependent SPDE, and finally, we can arrive at a system of equations according to the stepwise structure. For this purpose, we simulate the Gaussian random field with spatial covariance structure at a finite collection of predetermined collocation points. We then apply localized RBF approximation using PU structure on the spatial discretization. The matrices formed during the localized method will be sparse and, hence, will not suffer from ill-conditioning and high computational costs. RBF-PU-FD method is performed by using an implicit FD technique for the time discretization. Then, local RBF approximations on overlapping subdomains or patches that form a cover of the problem domain are weighted together by compactly supported PU weight functions to form the global approximation. When RBFs are used locally instead of globally, the computational cost is reduced because the previously dense linear systems then become sparse at the patch level. Here, a large problem is decomposed into many small problems and therefore, in the approximation procedure we can also work with a large number of nodes. Moreover, the convergence properties of the local approximations can be leveraged, while local couplings between approximations on different patches are enforced through the PU framework. For deterministic case, PU methods were proposed by Babuš ka and Melenk [2,33], combined with compactly supported RBFs (CSRBFs) [38], implemented by Cavoretto et al. for interpolation scattered data [4–6]. In this paper, we investigate the capability of RBF-PU-FD method for solving of time-dependent SPDEs in a two-dimensional space at each realization or simulation. RBF-PU-FD collocation method is the focus of the presented paper to reduce the computational cost of RBF-based methods through localization. The advantage of this approach is that a large problem is decomposed into many small problems. Therefore, in the approximation procedure, we can also deal with a large number of points. This method is proposed to cope with ill-conditioning, overcoming the traditional issues of global RBF-based methods. We will show that the method is viable through analyzing its numerical accuracy, CPU

RBF

𝜙𝜖 (r)

Gaussian C∞ (GA) Inverse MultiQuadric C∞ (IMQ) Maté rn C4 (M4) Maté rn C2 (M2) Wendland C4 (W4) Wendland C2 (W2)

𝑒− 𝜖 𝑟 (1 + 𝜖 2 𝑟2 )−1∕2 𝑒−𝜖𝑟 (𝜖 2 𝑟2 + 3𝜖𝑟 + 3) 𝑒−𝜖𝑟 (𝜖𝑟 + 1) (1 − 𝜖𝑟)6+ (35𝜖 2 𝑟2 + 18𝜖𝑟 + 3) (1 − 𝜖𝑟)4+ (4𝜖𝑟 + 1) 2 2

time, stability and sparsity structure. Finally, in order to show the efficiency of the proposed approach, we compare the results with the global RBF-FD method studied in [15]. The outline of the article is as follows. In the next section, we shortly describe the RBF approximation and also, we briefly review the PU collocation by RBFs. We propose the RBF-PU-FD collocation method for time-dependent SPDEs in Section 3. In Section 4, we investigate the stability of the resulting RBF-PU-FD scheme. In Section 5, the accuracy and efficiency of the proposed method are numerically investigated. Finally, in the last section, we present some conclusions and future work. 2. RBF-PU approximation In this paper, we will apply a spatially local variant of more traditional global RBF methods using PU based on a FD scheme for the numerical solution of time-dependent SPDEs. This RBF-PU-FD method has the advantage of low computational cost due to relatively sparse matrices. In this section, we will shortly describe a global approximation with RBFs and the PU collocation by RBFs. 2.1. RBF approximation Here, we shortly describe the RBF approximation. The technique of RBFs is one of the meshless methods that has been effectively worked for scattered data interpolation and approximation in high dimensions and complicated geometries [3]. Consider a domain Ω ⊆ ℝ𝑑 and a set of distinct points 𝑋 = {𝒙1 , … , 𝒙𝑁 } in Ω. In the interpolation of the scattered data using RBFs the approximation of a function u(x) on this set X may be written as a linear combination of N RBFs. The RBF interpolation problem is to find an interpolant 𝑠𝑢,𝑋 ∶ Ω → ℝ of the form 𝑠𝑢,𝑋 (𝒙) =

𝑁 ∑ 𝑗=1

𝛼𝑗 𝜙𝜖 (||𝒙 − 𝒙𝑗 ||),

𝒙 ∈ ℝ𝑑 ,

(1)

where 𝛼 j are unknown real coefficients, || · || denotes the Euclidean distance, and 𝜙𝜖 ∶ ℝ≥0 → ℝ is a RBF depending on a shape parameter 𝜖 > 0 such that 𝜙𝜖 (||𝒙 − 𝒛||) = 𝜙(𝜖||𝒙 − 𝒛||), for all x, z ∈ Ω. For simplicity, from now on we refer to 𝜙𝜖 as 𝜙. In Table 1, some well-known RBFs with their orders of smoothness which are commonly employed in the literature are listed [20]. We remark that Gaussian, Inverse MultiQuadric, and Maté rn functions are globally supported RBFs (GSRBFs) and strictly positive definite in ℝ𝑑 for any d, whereas Wendland ones are CSRBFs (whose support is [0, 1/𝜖]) and strictly positive definite in ℝ𝑑 for d ≤ 3 [37]. The coefficients 𝛼1 , … , 𝛼𝑁 are determined by enforcing the conditions 𝑠𝑢,𝑋 (𝒙𝑖 ) = 𝑢(𝒙𝑖 ),

𝑖 = 1, … , 𝑁.

Imposing these conditions leads to a linear system of equations 𝐴𝜶 = 𝒖,

(2) ℝ𝑁×𝑁

where the matrix 𝐴 ∈ is given by 𝐴𝑖𝑗 = 𝜙(||𝒙𝑖 − 𝒙𝑗 ||), 𝑖, 𝑗 = 1, … , 𝑁, 𝜶 = (𝛼1 , … , 𝛼𝑁 )𝑇 and 𝒖 = (𝑢(𝒙1 ), … , 𝑢(𝒙𝑁 ))𝑇 . 121

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

For a linear operator , we have from [18] that 𝑢(𝒙) ≃

𝑁 ∑ 𝑗=1

CSRBFs with support on Ωj such as the Wendland C2 function shown in Table 1. Clearly, the partition of unity property

𝛼𝑗 𝜙(||𝒙 − 𝒙𝑗 ||).

𝑀 ∑ 𝑗=1

In the sequel, for simplicity, we refer to the case of RBFs by positive definite RBFs, but the whole presentation could be done for conditionally positive definite RBFs. For the implementation of boundary conditions discussed later, we express the interpolant in Lagrange form, i.e., using cardinal basis functions. The cardinal basis functions, 𝜓 j (x), 𝑗 = 1, … , 𝑁, have the property { 1 if 𝑖 = 𝑗, 𝜓 𝑗 ( 𝒙𝑖 ) = 𝑗 = 1, … , 𝑁, 0 if 𝑖 ≠ 𝑗,

𝝍 𝑇 (𝒙)

where tain

𝑢,𝑋 (𝒙) =

= (𝜓1 (𝒙), … , 𝜓𝑁 (𝒙)). Combining (1), (3), and (2), we ob(4)

𝜓𝑗 (𝒙)𝑢(𝒙𝑗 ).

𝒔 = Ψ𝑇 𝒖 = Φ𝑇 𝐴−1 𝒖, where

In this subsection, we introduce the PU collocation method by RBFs to produce a conforming global approximation in terms of its weight functions and local RBF approximations. In the sequel, for simplicity, we refer to the case of PU collocation by positive definite RBFs, but the whole presentation could be done for conditionally positive definite RBFs. Let Ω ⊆ ℝ𝑑 be an open bounded domain. Let {Ω𝑗 }𝑀 be an open and 𝑗=1 bounded covering of Ω. This means all Ωj are open and bounded and Ω is contained in their union. Further, we define 𝐼(𝒙) = {𝑗 ∶ 𝒙 ∈ Ω𝑗 },

𝑘∈𝐼(𝒙) 𝜑𝑘 (𝒙)

𝜑𝑗 (𝒙) = 𝜙(

𝑅𝑗

),

𝑗∈𝐼(𝒙) 𝑘∈𝐽 (Ω𝑗 )

where C𝝁 is some positive constant, and 𝛿 j = diam(Ωj ). Definition 2.1. Let Ω ⊆ ℝ𝑑 be bounded and let 𝒙1 , … , 𝒙𝑁 ∈ 𝑋 ⊆ Ω be given. An open and bounded covering {Ω𝑗 }𝑀 is called regular for (Ω, 𝑗=1 X) if the following conditions hold:

𝑐𝑎𝑟𝑑 (𝐼(𝒙)) ≤ 𝐾,

,

𝑗 = 1, … , 𝑀,

• for each x ∈ Ω, the number of subdomains Ωj with x ∈ Ωj is bounded by a global constant K; • each subdomain Ωj satisfies an interior cone condition [39]; • the local fill distances ℎ𝑋𝑗 ,Ω𝑗 are uniformly bounded by the global fill distance

(5)

ℎ𝑋,Ω = sup min ||𝒙 − 𝒙𝑖 ||.

where 𝜑j (x) are compactly supported functions on Ωj . To guarantee nonnegativity and compact support in Ωj , we define 𝜑j (x) in (5) as ||𝒙 − 𝝃 𝑗 ||

𝑤𝑗 (𝒙)𝑠𝑢𝑗 ,𝑋𝑗 (𝒙),

‖ 𝝁 ‖ |𝝁| ≤ 𝐶 𝝁 ∕ 𝛿𝑗 , ‖𝐷 𝑤𝑗 ‖ ‖ ‖𝐿∞ (Ω𝑗 )

where the constant K is independent of the number of patches M. For each subdomain the PU weight function 𝑤𝑗 ∶ Ω𝑗 → ℝ are constructed by using Shepard method given by 𝜑 𝑗 ( 𝒙)

𝑗∈𝐼(𝒙)

If we fix 𝒙 = 𝒙𝑖 and k in Eq. (8), we get the ik-element of the global differentiation matrix corresponding to the 𝝁-derivative. For composite linear operators, we take the contributions from each term. To be able to formulate error bounds, we first consider some technical conditions and then define a few assumptions on regularity of Ωj [38]. So we require the partition of unity functions wj to be k-stable. This means that each 𝑤𝑗 ∈ 𝐶 𝑘 (ℝ𝑑 ) satisfies, for every multi-index 𝝁 ∈ ℕ𝑑0 with |𝝁| ≤ k, the inequality

2.2. PU collocation by RBFs

𝑤 𝑗 ( 𝒙) = ∑

∑

𝑗

= (𝜙(||𝒙𝑖 − 𝒙𝑗 ||)), 𝑖, 𝑗 = 1, … , 𝑁.

∀𝑥 ∈ Ω

𝑗=1

𝑤𝑗 (𝒙)𝑠𝑢𝑗 ,𝑋𝑗 (𝒙) =

In a differential problem using the PU approximation (7) needs to compute the effect of applying a spatial differential operator  at the interior data points. Denoting first by 𝝁 and 𝝂 the multi-indices for common rules and using then Leibniz’s rule, we can compute a derivative term of order 𝝁 of the global fit (7), obtaining the derivative rule ∑ ∑ 𝜕 |𝝁| ( ) 𝜕 |𝝁| 𝑤 𝑗 ( 𝒙) 𝜓 𝑘 ( 𝒙) 𝑢 𝑘  ( 𝒙) = 𝝁 𝜕 𝒙𝝁 𝑢,𝑋 𝜕 𝒙 𝑗∈𝐼(𝒙) 𝑘∈𝐽 (Ω𝑗 ) ( ) ( ) ∑ ∑ ∑ 𝝁 𝜕 |𝝁−𝝂| 𝑤𝑗 𝜕 |𝝂| 𝜓𝑘 = ( 𝒙) (𝒙) 𝑢𝑘 . (8) 𝝂 𝜕 𝒙𝝁−𝝂 𝜕 𝒙𝝂 𝑗∈𝐼(𝒙) 𝑘∈𝐽 (Ω ) 𝝂≤𝝁

To evaluate 𝑠𝑢,𝑋 (𝒙) at the nodes, i.e., to evaluate 𝒔 = (𝑠𝑢,𝑋 (𝒙1 ), … , 𝑠𝑢,𝑋 (𝒙𝑁 ))𝑇 , we need the differentiation matrix Ψ𝑇 = (𝜓𝑗 (𝒙𝑖 )), 𝑖, 𝑗 = 1, … , 𝑁. From relation (4), this leads to Φ𝑇

𝑀 ∑

𝑗∈𝐼(𝒙)

This transformation is valid whenever the matrix A is nonsingular, i.e. for distinct data points 𝒙1 , … , 𝒙𝑁 ∈ 𝑋 and positive definite RBFs. In a similar representation as (3), for any linear partial differential operator , 𝑢 may be approximated by

𝑗=1

𝒙 ∈ Ω,

where 𝐽 (Ω𝑗 ) = {𝑘 ∶ 𝒙𝑘 ∈ Ω𝑗 } that defines the set of nodes in Ωj and 𝜓 k are cardinal basis functions. Hence, in the PU, we find for the global approximation ∑ ∑ ∑ 𝑢,𝑋 (𝒙) = 𝑤𝑗 (𝒙)𝑠𝑢𝑗 ,𝑋𝑗 (𝒙) = 𝑤𝑗 (𝒙)𝜓𝑘 (𝒙)𝑢𝑘 . (7)

𝝍 𝑇 (𝒙) = 𝝓𝑇 (𝒙)𝐴−1 .

𝑁 ∑

𝑤 𝑗 ( 𝒙) = 1 ,

𝑘∈𝐽 (Ω𝑗 )

where 𝝓𝑇 (𝒙) = (𝜙(||𝒙 − 𝒙1 ||), … , 𝜙(||𝒙 − 𝒙𝑁 ||)). It is clear that relation (4) leads to the following relation between the cardinal basis and the original radial basis

𝑠𝑢,𝑋 (𝒙) =

𝑗∈𝐼(𝒙)

where 𝑠𝑢𝑗 ,𝑋𝑗 are local interpolations such that 𝑠𝑢𝑗 ,𝑋𝑗 (𝒙𝑘 ) = 𝑢𝑘 for all nodes xk ∈ Ωj . Then, the global PU approximation inherits the interpolation property of the local interpolations, i.e. 𝑢,𝑋 (𝒙𝑘 ) = 𝑢𝑘 . We construct the global approximation built up from local RBF interpolations of type (4). For each patch Ωj , the local RBF approximation 𝑠𝑢𝑗 ,𝑋𝑗 ∶ Ω𝑗 → ℝ is given by ∑ 𝑠𝑢𝑗 ,𝑋𝑗 (𝒙) = 𝜓 𝑘 ( 𝒙) 𝑢 𝑘 , (6)

(3)

𝑠𝑢,𝑋 (𝒙) = 𝝍 𝑇 (𝒙)𝒖 = 𝝓𝑇 (𝒙)𝜶 = 𝝓𝑇 (𝒙)𝐴−1 𝒖,

∑

is satisfied since 𝑤𝑗 (𝒙) = 0, ∀𝑗 ∉ 𝐼(𝒙). The global approximation function 𝑢,𝑋 (𝒙) in Ω to the solution function u(x) using the PU approximation is constructed as

leading to the alternative formulation for the interpolant 𝑠𝑢,𝑋 (𝒙) = 𝝍 𝑇 (𝒙)𝒖,

𝑤 𝑗 ( 𝒙) =

𝒙∈Ω 𝒙𝑖 ∈𝑋

After defining the space 𝐶𝜈𝑘 (ℝ𝑑 ) of all functions u ∈ Ck whose derivatives of order |𝝁| = 𝑘 satisfy 𝐷𝝁 𝑢(𝒙) = (||𝒙||𝜈 ) for ||x||2 → 0, we consider the following convergence result [39].

𝑗 = 1, … , 𝑀,

where {𝝃 𝑗 }𝑀 and {𝑅𝑗 }𝑀 are the centre and radius of the circular, 𝑗=1 𝑗=1 spherical or hyper-spherical patches

{Ω𝑗 }𝑀 𝑗=1

Theorem 2.1. Let Ω ⊆ ℝ𝑑 be open and bounded and 𝒙1 , … , 𝒙𝑁 ∈ 𝑋 ⊆ Ω. Let 𝜙 ∈ 𝐶𝜈𝑘 (ℝ𝑑 ) be a conditionally positive definite function of order m. If

and where 𝜙 is one of the

122

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

{Ω𝑗 }𝑀 is a regular covering for (Ω, X) and {𝑤𝑗 }𝑀 is k-stable for {Ω𝑗 }𝑀 , 𝑗=1 𝑗=1 𝑗=1 then the error between 𝑢 ∈ 𝜙 (Ω) and its PU interpolant 𝑢,𝑋 is bounded by (𝑘+𝜈)∕2−|𝝁|

|𝐷𝝁 𝑢(𝒙) − 𝐷𝝁 𝑢,𝑋 (𝒙)| ≤ 𝐶ℎ𝑋,Ω

𝑢

Rearranging Eq. (12) (or Eq. (13)), we obtain ) ( (𝒙) + 𝜂 𝜅 ▵ 𝑢𝑛+1 (𝒙) + 𝜈 ⋅ ▿𝑢𝑛+1 (𝒙)

𝑛+1

= 𝑢𝑛 (𝒙) + 𝜁 (𝜅 ▵ 𝑢𝑛 (𝒙) + 𝜈 ⋅ ▿𝑢𝑛 (𝒙)) + 𝜉(𝒙),

|𝑢|𝜙 (Ω) ,

(14)

or

for all x ∈ Ω, and all |𝝁| ≤ k/2.

𝑢𝑛+1 (𝒙) + 𝜂𝑢𝑛+1 (𝒙) = 𝑢𝑛 (𝒙) + 𝜁 𝑢𝑛 (𝒙) + 𝜉(𝒙),

If we compare the convergence result reported in Theorem 2.1 with the global error estimates in [39], we can note that the PU interpolant preserves the local approximation order for the global fit. So we can efficiently compute large RBF interpolants problem by solving many small RBF interpolation problems and then glue them together with the global partition of unity. It follows that the PU-based approach is a simple and computationally efficient tool to decompose a large interpolation problem into many small problems, while at the same time ensuring that the accuracy obtained for the local fits is carried over to the global one.

where 𝜂 = −𝜃𝛿𝑡, 𝜁 = (1 − 𝜃)𝛿𝑡. It follows from (14) and the definition of Brownian motion that the noise increment 𝜉(x) at each time instance 𝑡𝑛+1 is independent of the solution un (x) at the previous step, we simulate the Gaussian field with covariance structure 𝜎 2 𝛿tq(x, y) at a finite collection of predetermined collocation points [12]. Referring now to Eq. (7), we get that un (x) can be approximated by 𝑢 𝑛 ( 𝒙) ≃ 𝑝 𝑛 ( 𝒙) =

𝐴 = 𝐴ℑ + 𝐴𝔅 , 𝐴 = [𝑤𝑗 (𝒙𝑖 )𝜓𝑘 (𝒙𝑖 ) for (𝑗 ∈ 𝐼(𝒙𝑖 ), 𝑘 ∈ 𝐽 (Ω𝑗 ), 𝑖 = 1, … , 𝑁) and 0 elsewhere]𝑁×𝑁 , [ ] 𝐴ℑ = 𝑎𝑖𝑗 for (𝑖 ∈ ℑ, 1 ≤ 𝑗 ≤ 𝑁) and 0 elsewhere , [ ] 𝐴𝔅 = 𝑎𝑖𝑗 for (𝑖 ∈ 𝔅, 1 ≤ 𝑗 ≤ 𝑁) and 0 elsewhere .

Let us consider the following stochastic convection-diffusion equation 𝜕𝑊 (𝒙, 𝑡) 𝜕𝑢(𝒙, 𝑡) = 𝜅 ▵ 𝑢(𝒙, 𝑡) + 𝜈 ⋅ ▿𝑢(𝒙, 𝑡) + 𝜎 𝜕𝑡 𝜕𝑡 𝜕𝑊 (𝒙, 𝑡) ≡ 𝑢(𝒙, 𝑡) + 𝜎 , 𝒙 ∈ Ω ⊂ ℝ𝑑 , 𝑡 > 0, (9) 𝜕𝑡

Using the notation 𝐴 to designate the matrix of the same dimension as A, whose elements are 𝑎̂𝑖𝑗 = 𝑎𝑖𝑗 , 1 ≤ 𝑖, 𝑗 ≤ 𝑁, and substituting (16) in Eq. (14) (or Eq. (15)) together with (11), we can write the resulting sparse system, in the matrix form, as 𝐶 𝒖𝑛+1 = 𝐷𝒖𝑛 + 𝒗𝑛+1 + 𝝃,

where  is the convection-diffusion operator, △ and ▽ denote the Laplacian and the gradient operator, respectively, 𝜅 is the diffusion coefficient, 𝜈 is a constant velocity vector, 𝜎 > 0 and u(x, t) may represent concentration or temperature for mass or heat transfer [34]. We assume that W(x, t) is a continuous Wiener process defined on filtered probability space (Ω𝑊 , 𝑊 , { }∞ , ℙ𝑊 ) with mean zero and spatial covariance 𝑡=0 kernel q(x, y) given by

𝐶 = 𝐴 + 𝜂𝜅 ▵ 𝐴ℑ + 𝜂𝜈 ⋅ ▿𝐴ℑ , 𝐷 = 𝐴ℑ + 𝜁 𝜅 ▵ 𝐴ℑ + 𝜁 𝜈 ⋅ ▿𝐴ℑ , [ ]𝑇 𝒗𝑛+1 = 𝑔 𝑛+1 (𝒙𝑖 ) for (𝑖 ∈ 𝔅) and 0 elsewhere , [ ]𝑇 𝝃 = 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) and 0 elsewhere , ( )𝑇 𝒖𝑛 = 𝑢𝑛1 , … , 𝑢𝑛𝑁 . [ ]𝑇 In (18), the random vector 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) is the random right hand [ ]𝑇 side that varies with each realization of 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) [23]. In ad[ ]𝑇 dition, 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) is a random vector which has a multivariate normal distribution, i.e. [ ]𝑇 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) ∼  (𝟎, 𝑸), where 𝑸 = 𝜎 2 𝛿𝑡𝑞 (𝒙𝒊 , 𝒙𝒋 ) for (𝑖, 𝑗 ∈ ℑ).

𝑡, 𝑠 > 0.

(10)

and boundary conditions 𝒙 ∈ 𝜕Ω,

𝑡 > 0,

As a kernel function for the covariance operator q(x, y), we set 1 𝑞(𝒙, 𝒚 ) = . (‖𝒙 − 𝒚 ‖2 + 1)2

(11)

where g(x, t) is a known function, and it is given such that problem (9) has a unique solution [10,11], and  is a boundary operator, which can be of Dirichlet, Neumann or mixed type, and 𝜕Ω denotes the boundary of Ω. In the case of Dirichlet boundary conditions, we may discretize time derivative of the SPDE (9) by usual FD formula, using the following 𝜃-weighted scheme ( ) 𝑢𝑛+1 (𝒙) − 𝑢𝑛 (𝒙) 𝜎𝛿𝑊 𝑛+1 (𝒙) − = 𝜃 𝜅 ▵ 𝑢𝑛+1 (𝒙) + 𝜈 ⋅ ▿𝑢𝑛+1 (𝒙) 𝛿𝑡 𝛿𝑡 + (1 − 𝜃)(𝜅 ▵ 𝑢𝑛 (𝒙) + 𝜈 ⋅ ▿𝑢𝑛 (𝒙)),

The system (18) is obtained by combining Eq. (14) (or Eq. (15)), which applies to the internal points, and Eq. (11) that refers to the boundary points. Therefore, using the initial condition represented by Eq. (10), we can compute 𝒖𝑛+1 by solving Eq. (18). Then, by substituting such values of un in 𝒑𝑛 = 𝐴 𝒖𝑛 , 𝒑𝑛

(19) (𝑝𝑛 (𝒙

), … , 𝑝𝑛 (𝒙

))𝑇 ,

where = we can obtain the approximated solu1 𝑁 tion of the SPDE at time level n. Nevertheless, as well as for the RBF collocation scheme [25], also in the RBF-PU-FD method case we do not have a theoretical proof that shows the matrix C is invertible when 𝜃 > 0. So, here, we deal with the problem numerically, postponing the treatment of more theoretical issues such as well-posedness to future works. For the explicit scheme, i.e. for 𝜃 = 0, we only need to invert the matrix A, whose invertibility is guaranteed provided that the set of collocation points is distinct.

(12)

or equivalently 𝑢𝑛+1 (𝒙) − 𝑢𝑛 (𝒙) 𝜎𝛿𝑊 𝑛+1 (𝒙) − = 𝜃𝑢𝑛+1 (𝒙) + (1 − 𝜃)𝑢𝑛 (𝒙), 𝛿𝑡 𝛿𝑡

(18)

where

The Eq. (9) must be supplemented with an initial condition of the form

𝑢(𝒙, 𝑡) = 𝑔(𝒙, 𝑡),

(17)

where

3.1. Collocation scheme for a stochastic convection-diffusion problem

𝑢(𝒙, 0) = 𝑢0 (𝒙),

(16)

If we assume that ℑ and 𝔅 denote the indexes of internal and boundary points, respectively, and suppose that N is the total number of centers, i.e. 𝑁 = 𝑁ℑ + 𝑁𝔅 , then the N × N matrix A can be split into two matrices 𝐴ℑ and 𝐴𝔅 so that

In this section, we apply the RBF partition of unity collocation method based on FD for two time-dependent SPDEs. Specifically, we present our collocation scheme to firstly solve an stochastic convectiondiffusion equation in Section 3.1, and then a stochastic pseudo-parabolic equation in Section 3.2.

𝒙, 𝒚 ∈ Ω,

∑ ( ) 𝑤𝑗 (𝒙)𝜓𝑘 (𝒙) 𝑢𝑛𝑘 .

𝑗∈𝐼(𝒙) 𝑘∈𝐽 (Ω𝑗 )

3. RBF-PU collocation method based on FD

𝔼(𝑊 (𝒙, 𝑡), 𝑊 (𝒚 , 𝑠)) = 𝑚𝑖𝑛{𝑡, 𝑠}𝑞(𝒙, 𝒚 ),

∑

(15)

(13)

where 0 ≤ 𝜃 ≤ 1, 𝑢𝑛+1 (𝒙) = 𝑢(𝒙, 𝑡𝑛+1 ), 𝛿𝑊 𝑛+1 (𝒙) = 𝑊 𝑛+1 (𝒙) − 𝑊 𝑛 (𝒙), 𝑡𝑛+1 = 𝑡𝑛 + 𝛿𝑡, and 𝛿t is the time step size. Let 𝜉(𝒙) = 𝜎𝛿𝑊 𝑛+1 (𝒙). Here 𝜉(x) is a Gaussian field with 𝔼(𝜉(𝒙)) = 0 and 𝔼(𝜉(𝒙), 𝜉(𝒚 )) = 𝜎 2 𝛿𝑡𝑞 (𝒙, 𝒚 ). 123

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

3.2. Collocation scheme for a stochastic pseudo-parabolic problem

The associated system (18) is obtained by combining Eq. (24), which refers to the internal points, and Eq. (22) that applies to the boundary points. Therefore, using the initial condition represented by Eq. (21), we solve Eq. (18) computing 𝒖𝑛+1 . Finally, replacing such values of un in (19), we get the solution of the SPDE at time level n. In conclusion, however, all the remarks about the invertibility issue of the matrix C done previously for the convection-diffusion equation can be extended to this differential problem as well.

We consider the following third order stochastic equation 𝜕 ▵ 𝑢 ( 𝒙, 𝑡 ) 𝜕𝑊 (𝒙, 𝑡) 𝜕𝑢(𝒙, 𝑡) − 𝛼 ▵ 𝑢(𝒙, 𝑡) − 𝛽 = 𝑓 ( 𝒙, 𝑡 ) 𝜎 , 𝜕𝑡 𝜕𝑡 𝜕𝑡 𝑑 𝒙 ∈ Ω ⊂ ℝ , 𝑡 > 0,

(20)

where △ denotes the Laplacian operator, 𝛼, 𝛽 > 0 are known constants, 𝜎 > 0 and f(x, t) is given continuous function. We express W(x, t) similar to those explained in the previous subsection. The above equation is usually called stochastic Sobolev-type or stochastic pseudo-parabolic equation. To Eq. (20) we add initial condition 𝑢(𝒙, 0) = 𝑢0 (𝒙),

4. Stability analysis In this section, we provide a numerical stability analysis of the RBF-PU-FD collocation scheme considering explicitly the two studied stochastic time-dependent equations. We analyze the stability of the method using mean-square stability analysis [28]. ( ) Definition 4.1. A method is mean-square stable if 𝑙𝑖𝑚𝑛→∞ 𝔼 𝒑𝑛 𝒑𝑛 𝑇 = 0 where 𝔼 denotes the expectation value.

(21)

and boundary conditions 𝑢(𝒙, 𝑡) = 𝑔(𝒙, 𝑡),

𝒙 ∈ 𝜕Ω,

𝑡 > 0,

(22)

where g(x, t) is a known continuous function, and  is a boundary operator, which can be of Dirichlet, Neumann or mixed type, and 𝜕Ω denotes the boundary of Ω. In the case of Dirichlet boundary conditions, we can discretize time derivative of Eq. (20) via 𝜃-weighted FD scheme as ▵ 𝑢𝑛+1 (𝒙)− ▵ 𝑢𝑛 (𝒙) 𝜎𝛿𝑊 𝑛+1 (𝒙) 𝑢𝑛+1 (𝒙) − 𝑢𝑛 (𝒙) −𝛽 − 𝛿𝑡 𝛿𝑡 𝛿𝑡 ) ( = (1 − 𝜃)(𝛼 ▵ 𝑢𝑛 (𝒙) + 𝑓 𝑛 (𝒙)) + 𝜃 𝛼 ▵ 𝑢𝑛+1 (𝒙) + 𝑓 𝑛+1 (𝒙) , 𝑢𝑛+1 (𝒙)

𝑢(𝒙, 𝑡𝑛+1 ),

𝑛+1 (𝒙)

First of all, without loss of generality, suppose that 𝒗𝑛+1 ∈ ℝ𝑁×1 is the zero vector in the previous section. To apply this technique to evaluate the stability region of RBF-PU-FD method aforementioned in Eq. (18), we note that we can rewrite Eq. (18) in the form: 𝐶𝐴−1 𝒑𝑛+1 = 𝐷𝐴−1 𝒑𝑛 + 𝝃.

𝑛+1 (𝒙)

(23)

From (25) it follows that

𝑛 ( 𝒙) ,

𝒑𝑛+1 = 𝐴𝐶 −1 𝐷𝐴−1 𝒑𝑛 + 𝐶𝐴−1 𝝃. ⏟⏞⏞⏞⏞⏟⏞⏞⏞⏞⏟ ⏟⏟⏟

where 0 ≤ 𝜃 ≤ 1, = 𝛿𝑊 =𝑊 −𝑊 𝑡𝑛+1 = 𝑡𝑛 + 𝛿𝑡, and 𝛿t is the time step size. We express 𝜉(x) similar to those explained in the previous subsection. Eq. (23) can then be rewritten in the following form

𝐾

(24) 𝛿𝑡(𝜃𝑓 𝑛+1 (𝒙) +

𝜃)𝑓 𝑛 (𝒙)).

where 𝜂 = 𝜃𝛿𝑡𝛼, 𝜁 = (1 − 𝜃)𝛿𝑡𝛼, and = (1 − From (24) and the definition of Brownian motion, we can deduce considerations similar to those explained in the previous subsection. Similarly to the previous subsection, considering Eq. (7), un (x) can be approximated as in (16). Therefore, always denoting by ℑ and 𝔅 the indexes of internal and boundary points, respectively, and N being the total number of nodes so that 𝑁 = 𝑁ℑ + 𝑁𝔅 , the matrix 𝐴 ∈ ℝ𝑁×𝑁 is given by the sum of matrices 𝐴ℑ and 𝐴𝔅 as in (17), where

𝑄

𝐹𝑛+1

𝐹𝑛

𝐺

where 𝐼 ∈ ℝ𝑁×𝑁 is the identity matrix. On account of the symmetry of the matrices 𝐹𝑛+1 , KFn KT and G, we can write Eq. (27) as a linear system of the from 𝑓̃𝑛+1 = 𝐻 𝑓̃𝑛 + 𝜎 2 𝛿𝑡𝑔̃ 𝑛+1 ,

(28)

where 𝑓̃𝑛+1 , 𝑓̃𝑛 and 𝑔̃ are 𝑁(𝑁 + 1)-dimensional vectors consisting of the free components of 𝐹𝑛+1 , Fn and G, respectively. Also, H is a square matrix where matrix H is constructed by using elements of matrix K. Therefore, using the recessive relation represented by Eq. (28), we have (𝑛−1 ) ∑ 𝑛 𝑛 ̃0 2 ̃ 𝑓 = 𝐻 𝑓 + 𝜎 𝛿𝑡 𝐻 𝑖 𝑔̃ 𝑛+1 . (29)

𝐴 = [𝑤𝑗 (𝒙𝑖 )𝜓𝑘 (𝒙𝑖 ) for (𝑗 ∈ 𝐼(𝒙𝑖 ), 𝑘 ∈ 𝐽 (Ω𝑗 ), 𝑖 = 1, … , 𝑁) and 0 elsewhere]𝑁×𝑁 , [ ] 𝐴ℑ = 𝑎𝑖𝑗 for (𝑖 ∈ ℑ, 1 ≤ 𝑗 ≤ 𝑁) and 0 elsewhere , [ ] 𝐴𝔅 = 𝑎𝑖𝑗 for (𝑖 ∈ 𝔅, 1 ≤ 𝑗 ≤ 𝑁) and 0 elsewhere . Now, substituting (16) in Eq. (24) together with (22), we obtain a sparse matrix system of the form (18), with

𝑖=0

Eq. (29) demonstrates that 𝑙𝑖𝑚𝑛→∞ 𝑓̃𝑛 = 0 if the spectral radius 𝜌 of the matrix H is less than one. Since 𝑓̃𝑛 is the vector consisting of the free ( ) components of 𝔼 𝒑𝑛 𝒑𝑛 𝑇 , so referring to Definition 4.1, we can deduce that mean-square stability is assured if 𝜌(H) < 1. In Section 5, we will investigate numerical stability analysis for the two SPDEs discussed in this paper.

𝐶 = 𝐴 − (𝜂 + 𝛽) ▵ 𝐴ℑ , 𝐷 = 𝐴ℑ + (𝜁 − 𝛽) ▵ 𝐴ℑ , [ ]𝑇 𝒗𝑛+1 = 𝑧𝑛𝑖 +1 for (𝑖 ∈ ℑ) and 𝑔𝑗𝑛+1 for (𝑗 ∈ 𝔅) , [ ]𝑇 𝝃 = 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) and 0 elsewhere , ( 𝑛 )𝑇 𝑛 𝑛 𝒖 = 𝑢1 , … , 𝑢𝑁 .

5. Numerical experiments

[ ]𝑇 The random vector 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) is the random right hand side [ ]𝑇 that varies with each realization of 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) . In addition, [ ]𝑇 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) is a random vector which has a multivariate normal distribution, i.e. [ ]𝑇 𝜉(𝒙𝑖 ) for (𝑖 ∈ ℑ) ∼  (𝟎, 𝑸), where 𝑸 = 𝜎 2 𝛿𝑡𝑞 (𝒙𝒊 , 𝒙𝒋 ) for (𝑖, 𝑗 ∈ ℑ).

In this section, we report the performance of our RBF-PU-FD method which is measured through numerical experiments. All these results illustrated in some tables and figures have been carried out in Matlab on a laptop with a 2.4 GHz Intel Core i5 processor. In the following, we focus on a wide series of experiments, which concern the 2D stochastic convection-diffusion and stochastic pseudoparabolic problems studied in Section 3. In doing that, we discretize the domain Ω taking some sets of two different data point distributions, which consist of a number N of uniformly distributed grid nodes

As a kernel function for the covariance operator q(x, y), we set 𝑞(𝒙, 𝒚 ) =

(26)

Squaring and then taking the expectation of both sides of Eq. (26) coupled with the fact that 𝜉(x) is a Gaussian field with 𝔼(𝜉(𝒙)) = 0 and 𝔼(𝜉(𝒙), 𝜉(𝒚 )) = 𝜎 2 𝛿𝑡𝑞 (𝒙, 𝒚 ), we obtain ( ) ( ) 𝑇 𝔼 𝒑𝑛+1 𝒑𝑛+1 = 𝐾 𝔼 𝒑𝑛 𝒑𝑛 𝑇 𝐾 𝑇 + 𝜎 2 𝛿𝑡 𝑄𝑄𝑇 , (27) ⏟⏟⏟ ⏟⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏟ ⏟⏞⏞⏞⏟⏞⏞⏞⏟

𝑢𝑛+1 (𝒙) − (𝜂 + 𝛽) ▵ 𝑢𝑛+1 (𝒙) = 𝑢𝑛 (𝒙) + (𝜁 − 𝛽) ▵ 𝑢𝑛+1 (𝒙) + 𝑧𝑛+1 (𝒙) + 𝜉(𝒙), 𝑧𝑛+1

(25)

1 . (‖𝒙 − 𝒚 ‖2 + 1)2 124

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 1. Example of partitioning of “regular” and ”irregular” domains with circular patches for uniform and Halton points (left to right).

and quasi-random Halton points. The latter, which are a typical example of scattered node set, are generated by using the Matlab program haltonseq.m [20]. More precisely, in our study, we consider three types of “regular” (convex) and “irregular (non-convex) domains, i.e. the square domain Ω = [0, 1] × [0, 1], and the non-convex and circular domains Ω ⊆ [−1, 1] × [−1, 1], see Fig. 1. The PU covering is composed instead of M circular patches that are centered at a uniform grid of points, where the overlap of the patches is 20% of the distance between the centers [35]. Note that, in the PU scheme, the subdomains might be of any (regular enough) geometric shapes such as circles, squares, rectangles and pentagons. However, many other possible types of patches are allowed: in fact, the only requirement for all shapes of subdomains is that they cover the domain Ω. As said earlier, we here use circular subdomains so that any (possibly also mild) overlap among the different subdomains is guaranteed [35]. An example of uniform and Halton points along with the related PU subdomains for the square, non-convex and circular domains Ω is shown in Fig. 1, top to bottom.

We show the numerical results obtained by applying the PU collocation method based on FD using some of the RBFs contained in Table 1 as local approximants. More precisely, our analysis is based on considering GA, IMQ and M4 as basis functions in (6), and the compactly supported W2 as localizing function of Shepard’s weight in (5). Acting in this way, besides presenting diversified numerical tests, we compute the approximation errors, selecting “optimal” values of the shape parameter 𝜖. Such values of 𝜖 are searched for in the interval [0.05, 13.00], taking as step size 0.05. Notice that the accuracy of RBF-based methods highly depends upon the shape parameter 𝜖 of the basis functions, which is responsible for the flatness of such functions. In particular, for smooth problems, the best accuracy is typically achieved when 𝜖 is small, but then the condition number of the linear system may become very large. In the practice, even if the PU – as highlighted in the next subsections – reveals to be much better conditioned than traditional RBF techniques, the selection of shape parameters may greatly affect the accuracy of the collocation method. However, in some particularly hard situations, in which the 125

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 2. Numerical study of the stability of RBF-PU-FD method regarding the stochastic convection-diffusion equation on the square domain with 8 × 8 uniform (left) and Halton (right) nodes and 2 × 2 patches for IMQ with 𝛿𝑡 = 0.001, and 𝜃 = 1∕2.

ill-conditioning is so big to produce a very negative effect on the results, the ill-conditioning problem due to RBFs can be overtaken by using stable methods; see [21, Chapter 12] for an overview. Such approximation techniques allow stable computations for small values of 𝜖 as well. The numerical algorithm consists of three steps.

where ABS, MAE and RMS denote the absolute error, maximum absolute error, and root mean squared error, respectively. Note that u in (31) is the solution in the absence of Gaussian noise, i.e. the related deterministic solution. In fact, MAE and RMS show the difference between the mean solution and the solution in the absence of Gaussian noise. In the following examples, we assume initial and boundary conditions calculated from the exact solution in the absence of Gaussian noise. The stability is studied by evaluating the condition number (CN) obtained by using the Matlab command condest. Wecalculate the criterion of sparsity as the number of non-zero entries in the coefficient matrix divided by the number of entries in the coefficient matrix. Although the approximation scheme described in Section 3 is generally valid for any value of 𝜃 ∈ [0, 1], in the following, we focus on the case 𝜃 = 1∕2, which identifies the famous Crank-Nicholson scheme. We show that the method is viable through analyzing its numerical accuracy, CPU time, stability and sparsity structure. Finally, in order to show the efficiency of the proposed approach, we compare the results with the global RBF-FD method studied in [15]. Therefore, in the following sections we consider more in detail the two SPDE problems: first, we show results for the stochastic convection-diffusion equation, and then for the stochastic pseudo-parabolic one.

• Step 1: For 𝑖 = 1, … , 𝑠, approximate the Gaussian noise 𝜉(𝒙) ∼  (𝟎, 𝜎 2 𝛿𝑡𝑞 (𝒙, 𝒚 )) at the collocation points, where s is the number of iterations and 𝜎 2 𝛿tq(x, y) is the covariance functions of the Gaussian noise. [ ]𝑇 • Step 2: For 𝑖 = 1, … , 𝑠, solve Eq. (18), with 𝜉(𝒙𝑗 ) for (𝑗 ∈ ℑ) ∼ 2  (𝟎, 𝑸) replaced by 𝜉(x) where 𝑸 = 𝜎 𝛿𝑡𝑞 (𝒙𝒊 , 𝒙𝒋 ) for (𝑖, 𝑗 ∈ ℑ), to obtain un . Then, by substituting such values of un in (19), we can obtain the approximated solution p of the SPDE at time level n. • Step 3: suppose that 𝔼, Dev, Upp and Low denote the statistical criterions such as mean, standard deviation, lower bound ( ) and upper bound for prediction, respectively. Evaluate 𝔼 𝑝𝑛𝑇 (𝒙𝑗 ) , ( 𝑛 ) ( 𝑛 ) ( 𝑛 ) 𝑇 𝑇 𝑇 Dev 𝑝 (𝒙𝑗 ) , Upp 𝑝 (𝒙𝑗 ) and Low 𝑝 (𝒙𝑗 ) using the Monte-Carlo method: 𝑠 ( ) 1∑ 𝑛 𝔼 𝑝𝑛𝑇 (𝒙𝑗 ) ≈ 𝑝 𝑇 ( 𝒙𝑗 ) , 𝑠 𝑖=1 𝑖

(30)

√ √ 𝑠 ( 𝑛 ) √ ( ))2 1 ∑ ( 𝑛𝑇 Dev 𝑝 𝑇 (𝒙𝑗 ) ≈ √ 𝑝 ( 𝒙𝑗 ) − 𝔼 𝑝 𝑛 𝑇 ( 𝒙𝑗 ) , 𝑠 − 1 𝑖=1 𝑖

5.1. Results for stochastic convection-diffusion problem In this subsection, we show numerical results acquired from experiments carried out for the stochastic convection-diffusion equation (9) on the square and non-convex domains, see Fig. 1 (top and center). Taking an appropriate initial condition and Dirichlet boundary conditions, the analytical solution in the absence of Gaussian noise is given by

( ) ( ) ( ) Upp 𝑝𝑛𝑇 (𝒙𝑗 ) = 𝔼 𝑝𝑛𝑇 (𝒙𝑗 ) + 1.96Dev 𝑝𝑛𝑇 (𝒙𝑗 ) , ( ) ( ) ( ) Low 𝑝𝑛𝑇 (𝒙𝑗 ) = 𝔼 𝑝𝑛𝑇 (𝒙𝑗 ) − 1.96Dev 𝑝𝑛𝑇 (𝒙𝑗 ) ,

𝑢(𝑥, 𝑦, 𝑡) = 𝑎 exp(𝑏𝑡)(exp(−𝑐𝑥) + exp(−𝑐𝑦)),

𝑛

where T denotes the maximum time level, i.e. 𝑇 = 𝑡max = 1 and 𝑝𝑖 𝑇 is approximated solution 𝑝𝑛𝑇 of the SPDE at time T and ith realization. We examine the performance of the proposed method for solving stochastic convection-diffusion equations in two space dimension. For test problems, we perform 1000 realizations (or 𝑠 = 1000) and statistical criterions such as mean, variance andstandard deviation are computed. To see the convergence of the proposed schemes we also employ the following errors: ) | ( | ABS(𝒙𝑖 ) = |𝔼 𝑝𝑛𝑇 (𝒙𝑖 ) − 𝑢(𝒙𝑖 , 𝑇 )|, (31) | |

where constants a and b can be chosen freely, and √ 𝜈 ± 𝜈 2 + 4𝑏𝜅 𝑐= > 0. 2𝜅 Here, to exhibit our results after a wide experimentation with different parameters, we assume 𝑎 = 1, 𝑏 = 0.1 and 𝜎 = 1. Also, the convection velocity is chosen to be 𝜈 = 1 in most experiments, and 𝜅 = 1 is used as default diffusion coefficient, but other values are used as well. First of all, focusing on the square domain, we perform a numerical study of the stability of our method by investigating condition 𝜌(H) < 1 that depends on the eigenvalues 𝜆 of H matrix of the RBF-PU-FD method for the stochastic convection-diffusion problem. As an example, in Fig. 2 we highlight that condition |𝜆(H)| < 1 is satisfied for all eigenvalues of matrix H for 𝜃 = 1∕2. Therefore, these results also confirm from a numerical standpoint that our collocation method is unconditionally stable.

MAE = max ABS(𝒙𝑖 ), 1≤𝑖≤𝑁

√ √ 𝑁 √1 ∑ ( )2 RMS = √ ABS(𝒙𝑖 ) , 𝑁 𝑖=1 126

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 3. Sparsity structure of coefficient matrix regarding the stochastic convection-diffusion equation on the square domain with 24 × 24 uniform nodes and 3 × 3 patches (left) and 5 × 5 patches (right) for IMQ with 𝜖 = 1.80.

Table 2 Comparison between RBF-PU-FD and RBF-FD [15] using IMQ for the stochastic convection-diffusion equation on the square domain: MAE, RMS, CPU time (in seconds), CN and sparsity for uniform points with 𝛿𝑡 = 0.001 and 𝑇 = 1. √

𝑁

16 24 32 40

RBF-PU-FD √ 𝑀 𝜖

MAE

RMS

time

CN

sparsity

𝜖

MAE

RMS

time

CN

sparsity

2 3 4 5

7.01 × 10−3 3.81 × 10−3 5.10 × 10−3 7.33 × 10−3

2.56 × 10−3 1.71 × 10−3 2.25 × 10−3 3.28 × 10−3

1.83 5.93 19.41 34.58

4.03 × 10+01 1.00 × 10+02 2.45 × 10+02 3.59 × 10+02

5.16 × 10−1 2.93 × 10−1 1.88 × 10−1 1.29 × 10−1

1.55 2.40 3.25 4.25

1.63 × 10−2 1.87 × 10−2 1.13 × 10−2 2.04 × 10−2

6.77 × 10−3 8.60 × 10−3 4.69 × 10−3 9.05 × 10−3

8.62 30.54 73.46 367.19

4.86 × 10+17 6.68 × 10+17 4.48 × 10+17 3.20 × 10+17

1 1 1 1

1.35 1.80 2.85 3.60

RBF-FD

The RBF-FD method leads to a dense coefficient matrix in the resulting problem. Therefore, we use the RBF-PU-FD method which introduces sparsity in the matrices. In order to stress sparsity of the coefficient matrix, in Fig. 3 we show the sparse structure generated by the RBF-PU-FD collocation scheme for the stochastic convection-diffusion equation considering IMQ and two different domain partitions. The sparsity pattern of the matrix depends on the number of patches. More patches lead to more sparsity for the same number of nodes as can be seen in Fig. 3. Generally, a covering with smaller partitions leads to worse approximation results but, at the same time, it turns out to be less expensive from the computational standpoint. This is essentially due to the greater sparsity of the linear system. For a fixed number of nodal points, we compute the MAE for different numbers of patches. The number of patches corresponding to the the minimal MAE is then used in the numerical results. So when increasing the number of nodal points to refine the grid the number of patches will be changed accordingly. In Tables 2 we present an extensive and detailed analysis of the RBFPU-FD collocation method, which is also compared with the global RBFFD method studied in [15]. More precisely, considering four sets of uniform and Halton points and IMQ as local RBF approximation. Also, we report the MAE and RMS computed with “optimal” 𝜖, the CPU time required to solve the sparse linear system (18), the CN and the criterion of sparsity in the coefficient matrix. Similar results are obtained for the Halton points that have been removed to be shortened. Furthermore, the sparsity structure of cofficient matrix obtained for 40 × 40 uniform points contained in the square domain equals 2560000 number of zero (𝑛𝑧 = 2560000) with 𝜖 = 4.25 for RBF-FD method while for RBF-PU-FD method the number of zero equals 331315 (𝑛𝑧 = 331315) with 𝜖 = 3.60 and 5 × 5 patches, using IMQ with = 0 ∶ 001 and 𝑇 = 1. Also, in Fig. 4 graphically the surface of exact solution in the absence of Gaussian noise, absolute error, statistical criterions such as mean, standard deviation, lower bound and upper bound for prediction in given uniform points is represented. In Fig. 5, the convergence of our RBFPU-FD scheme has been investigated numerically for a fixed number of patches M, uniformly distributed nodes with varying the global fill distance hX,Ω are employed as shown such discussions in [30]. In all cases we fix 𝛿𝑡 = 0.001 and 𝑇 = 1. From this study we can note a quite uniform

behavior: on the one hand, accuracy of the RBF-FD method is almost similar, with usually a slight advantage for the RBF-PU-FD method, but on the other hand we observe a significant reduction of CPU time, CN and criterion of sparsity in coefficient matrix for the RBF-PU-FD method compared to the RBF-FD one. Furthermore, as shown in Fig. 5, increasing the number of local points for a fixed number of patches results high order convergence. Note that, in general, for the RBF-PU-FD scheme we found a good compromise between √ accuracy √ and efficiency assuming that on the square domain the ratio 𝑁 ∕ 𝑀 = 8. Obviously, the flexibility of our method enables also to consider different values as it can occur in some situations. To illustrate the applicability of our method on irregularly shaped domains, we solve the stochastic convection-diffusion problem on a nonconvex domain, whose region is shown in Fig. 1 (center). Specifically, we hereby define the boundary of this irregular domain Ω as follows: 𝜕Ω = {(𝑟, 𝜃) ∣ 𝑟(𝜃) = 1 +

1 (sin(6𝜃) + sin(3𝜃))}. 10

To cover the domain with patches, we use the partition covering scheme proposed in [24]. Hence, in Table 3 we exhibit results obtained with the RBF-PU-FD collocation method on the irregular domain, comparing them with the global RBF-FD method for both uniform points. Similar results are obtained for the Halton points that have been removed to be shortened. Though many experiments have been done by various RBFs, for brevity we are limited to show the results associated with the IMQ function, even so reporting all information provided so far in our previous tests. Furthermore, as an example, using sets of uniform and Halton nodes, in Figs. 6 and 7 we graphically represent surface of the exact solution in the absence of Gaussian noise, absolute error, statistical criterions such as mean, standard deviation for prediction for given uniform (left) and Halton (right) points. As the nonphysical pressure oscillation in the Stokes problem [44], numerical simulations of convection-diffusion problems with small diffusivity present more challenges. To show the efficiency of the proposed method to small diffusivity, 𝜅 = 0.1 is used. Therefore, taking a set of uniform points, in Table 4 we compare our RBF-PU-FD scheme with the global RBF-FD method. Analyzing these experiments, we can thus

127

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 4. Graphs of the exact solution in the absence of Gaussian noise, mean solution, absolute error, standard deviation, lower bound and upper bound of prediction obtained for 40 × 40 uniform points, 𝜖 = 3.60 with 5 × 5 patches contained in the square domain, using IMQ with 𝛿𝑡 = 0.001 and 𝑇 = 1. Table 3 Comparison between RBF-PU-FD and RBF-FD [15] using IMQ for the stochastic convection-diffusion equation on the non-convex domain: MAE, RMS, CPU time (in seconds), CN and sparsity for uniform points with 𝛿𝑡 = 0.001 and 𝑇 = 1. RBF-PU-FD

RBF-FD

N

M

𝜖

MAE

RMS

time

CN

sparsity

𝜖

MAE

RMS

time

CN

sparsity

289 639 1119 1729

6 10 18 27

0.80 1.35 1.75 2.35

7.61 × 10−3 9.84 × 10−3 1.00 × 10−2 1.01 × 10−2

2.77 × 10−3 4.83 × 10−3 4.54 × 10−3 4.04 × 10−3

1.83 8.95 18.35 58.91

1.10 × 10+02 1.87 × 10+03 2.14 × 10+03 5.27 × 10+02

4.32 × 10−1 3.73 × 10−1 1.96 × 10−1 1.41 × 10−1

0.90 1.45 2.05 2.60

1.68 × 10−2 1.98 × 10−2 4.31 × 10−2 1.98 × 10−2

7.45 × 10−3 7.87 × 10−3 1.82 × 10−2 9.19 × 10−3

10.01 36.03 89.51 452.09

4.58 × 10+17 2.86 × 10+18 2.30 × 10+17 1.34 × 10+17

1 1 1 1

Table 4 Comparison between RBF-PU-FD and RBF-FD [15] using IMQ for the stochastic convection-diffusion equation on the square domain: MAE, RMS, CPU time (in seconds), CN and sparsity for uniform points with 𝜅 = 0.1, 𝛿𝑡 = 0.001 and 𝑇 = 1. √ 𝑁 16 24 32 40

RBF-PU-FD √ 𝑀 𝜖

MAE

2 3 4 5

1.47 × 10 1.67 × 10−2 2.53 × 10−2 1.11 × 10−2

1.35 1.80 2.85 3.60

RBF-FD RMS −2

6.41 × 10 5.38 × 10−3 1.04 × 10−2 4.98 × 10−3 −3

time

CN

sparsity

1.66 7.02 10.92 31.37

6.89 × 10 2.98 × 10+01 3.26 × 10+01 1.34 × 10+02 +00

5.16 × 10 2.93 × 10−1 1.88 × 10−1 1.29 × 10−1 −1

128

𝜖

MAE

1.65 2.50 3.55 4.40

3.54 × 10 1.92 × 10−2 1.65 × 10−2 3.74 × 10−2

RMS −2

1.55 × 10 7.93 × 10−3 7.30 × 10−3 1.35 × 10−2 −2

time

CN

sparsity

7.97 31.52 81.35 378.87

1.87 × 10 5.46 × 10+17 3.16 × 10+17 7.34 × 10+17 +17

1 1 1 1

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 6. Surface of the exact solution in the absence of Gaussian noise of (9) given 1762 Halton points contained in the non-convex domain, with 𝑇 = 1.

observe a behavior in terms of MAE, RMS, CPU time, CN and sparsity, similar to that exhibited for 𝜅 = 1 and already beforehand remarked. In addition, taking set of uniform points and using IMQ, Figs. 8 and 9

Fig. 5. The convergence of our RBF-PU-FD scheme as a function of hX,Ω for a fixed number of patches M, uniform points contained in the square domain, using IMQ with 𝛿𝑡 = 0.001 and 𝑇 = 1.

Fig. 7. Graphs of approximate value of mean solution, standard deviation and absolute error obtained and for 1729 uniform points, 𝜖 = 2.35 (left) and 1752 Halton points, 𝜖 = 2.15 (right) with 27 patches contained in the non-convex domain, using IMQ with 𝛿𝑡 = 0.001 and 𝑇 = 1. 129

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 8. Graphs of the approximate value of mean solution, absolute error, lower bound and upper bound of prediction obtained for 40 × 40 uniform points, 𝜖 = 3.60 with 5 × 5 patches contained in the square domain, using IMQ with 𝜅 = 0.1, 𝛿𝑡 = 0.001 and 𝑇 = 1.

while the function f that appears in (20) is given by 𝑓 (𝑥, 𝑦, 𝑡) = (2 + 𝛼 + 2𝛽) exp(2𝑡)(cos(𝑥) + sin(𝑦)). The experiments we report have been performed assuming as parameters the values of 𝛼 = 1, 𝛽 = 0.00025 and 𝜎 = 1. At first, considering again – as at the beginning of the previous subsection – the square domain, we analyze numerically the stability of our RBF-PU-FD scheme by verifying the validity of the inequality 𝜌(H) < 1 for the stochastic pseudo-parabolic equation. In particular, we can show that condition |𝜆(H)| < 1 is satisfied for all eigenvalues of matrix H for 𝜃 = 1∕2 such that our collocation scheme turns out to be unconditionally stable. Now, we compare the sparse structure of the matrix generated by using the RBF-PU-FD collocation method with local IMQ approximations and two different domain partitions for the pseudo-parabolic equation. The sparsity structure of cofficient matrix on the square domain with 24 × 24 Halton nodes and 3 × 3 and 5 × 5 patches with 𝜖 = 1.85 equal to 𝑛𝑧 = 94754 and 𝑛𝑧 = 39192, respectively. As in the stochastic convectiondiffusion case, although the matrix structure is different, we can derive similar conclusions. In Tables 5 and 6 we then compare our RBF-PU-FD scheme with the RBF-FD method studied in [15]. Therefore, for some sets of uniform points, we show the numerical results obtained by using IMQ and GA as local RBF approximants. To be shortened the obtained results for the Halton points removed. Specifically, we report the MAE and RMS with the associated “optimal” value of 𝜖, the CPU time required to solve the sparse linear system (18), the CN and the criterion of sparsity in the coefficient matrix. Furthermore, the sparsity structure of cofficient matrix obtained for 40 × 40 uniform points contained in the square domain equals 2560000 number of zero (𝑛𝑧 = 2560000) with 𝜖 = 4.15 for RBF-FD method while for RBF-PU-FD method the number of zero equals 331315 (𝑛𝑧 = 331315) with 𝜖 = 3.35 with 5 × 5 patches. In addition, taking a set of uniform points and using IMQ and GA, Figs. 10–12 shows the surface of exact solution in the absence of Gaussian noise, absolute error, statistical criterions such as mean, standard deviation, lower

Fig. 9. Surface of the standard deviation of prediction obtained for 40 × 40 uniform points, 𝜖 = 3.60 with 5 × 5 patches contained in the square domain, using IMQ with 𝜅 = 0.1, 𝛿𝑡 = 0.001 and 𝑇 = 1.

shows the surface of the absolute error, statistical criterions such as mean, standard deviation, lower bound and upper bound for prediction for given uniform points. 5.2. Results for stochastic pseudo-parabolic problem In this subsection, we present some of the numerical experiments we made to test our method for the stochastic pseudo-parabolic equation (20) on the square and circular domains see Fig. 1 (top and bottom). Thus, taken a finite domain Ω, we consider Eqs. (20)–(22) with initial condition 𝑢(𝑥, 𝑦, 0) = cos(𝑥) + sin(𝑦),

(𝑥, 𝑦) ∈ Ω.

The exact solution to the stochastic pseudo-parabolic equation in the absence of Gaussian noise is 𝑢(𝑥, 𝑦, 𝑡) = exp(2𝑡)(cos(𝑥) + sin(𝑦)), 130

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Table 5 Comparison between RBF-PU-FD and RBF-FD [15] using IMQ for the stochastic pseudo-parabolic equation on the square domain: MAE, RMS, CPU time (in seconds), CN and sparsity for uniform points with 𝛿𝑡 = 0.001 and 𝑇 = 1. √

𝑁

16 24 32 40

RBF-PU-FD √ 𝑀 𝜖

MAE

2 3 4 5

7.56 × 10 5.24 × 10−3 5.97 × 10−3 9.65 × 10−3

1.35 1.80 2.85 3.35

RBF-FD RMS −3

3.56 × 10 2.33 × 10−3 2.78 × 10−3 4.74 × 10−3 −3

time

CN

sparsity

1.27 7.22 17.12 44.79

5.02 × 10 1.49 × 10+02 3.15 × 10+02 6.85 × 10+02 +01

5.16 × 10 2.93 × 10−1 1.88 × 10−1 1.29 × 10−1 −1

𝜖

MAE

1.60 2.50 3.40 4.15

1.87 × 10 2.10 × 10−2 1.32 × 10−2 1.10 × 10−2

RMS −2

1.10 × 10 9.49 × 10−3 6.26 × 10−3 4.48 × 10−3 −2

time

CN

sparsity

7.91 32.06 81.13 350.54

1.23 × 10 4.49 × 10+17 4.24 × 10+17 6.72 × 10+17 +17

1 1 1 1

Table 6 Comparison between RBF-PU-FD and RBF-FD [15] using GA for the stochastic pseudo-parabolic equation on the square domain: MAE, RMS, CPU time (in seconds), CN and sparsity for uniform points with 𝛿𝑡 = 0.001 and 𝑇 = 1. √ 𝑁 16 24 32 40

RBF-PU-FD √ 𝑀 𝜖

MAE

2 3 4 5

4.92 × 10 2.81 × 10−3 7.64 × 10−3 7.38 × 10−3

3.20 4.95 7.45 9.20

RBF-FD RMS −3

2.02 × 10 1.39 × 10−3 3.52 × 10−3 2.62 × 10−3 −3

time

CN

sparsity

2.51 6.32 16.53 44.58

4.14 × 10 3.67 × 10+03 2.00 × 10+03 2.87 × 10+05 +02

5.16 × 10 2.93 × 10−1 1.88 × 10−1 1.29 × 10−1 −1

𝜖

MAE

4.00 6.90 9.85 12.65

1.10 × 10 1.68 × 10−2 2.06 × 10−2 3.56 × 10−2

RMS −2

5.34 × 10 9.56 × 10−3 7.06 × 10−3 2.56 × 10−2 −3

time

CN

sparsity

6.84 23.31 57.53 197.38

6.20 × 10 1.81 × 10+18 2.66 × 10+17 2.96 × 10+17 +17

1 1 1 1

Fig. 10. Graphs of the exact solution in the absence of Gaussian noise, mean solution, absolute error, standard deviation, lower bound and upper bound of prediction obtained for 40 × 40 uniform points, 𝜖 = 3.35 with 5 × 5 patches contained in the square domain, using IMQ with 𝛿𝑡 = 0.001 and 𝑇 = 1.

131

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 11. Graphs of approximate value of mean solution, absolute error, lower bound and upper bound of prediction obtained for 40 × 40 uniform points and 𝜖 = 9.20 with 5 × 5 patches contained in the square domain, using GA with 𝛿𝑡 = 0.001 and 𝑇 = 1.

Fig. 12. Surface of the standard deviation of prediction obtained for 40 × 40 uniform points, 𝜖 = 9.20 with 5 × 5 patches contained in the square domain, using GA with 𝛿𝑡 = 0.001 and 𝑇 = 1.

Fig. 13. Surface of the exact solution in the absence of Gaussian noise of (9) given 1709 Halton points contained in the circular domain, with 𝑇 = 1.

bound and upper bound for prediction for given uniform points. The convergence of our RBF-PU-FD scheme can be investigated numerically for a fixed number of patches M such that we can employ uniformly distributed nodes with varying the global fill distance hX,Ω . As time parameters we choose again the values of 𝛿𝑡 = 0.001 and 𝑇 = 1. From such tables, we observe that the accuracy of the RBF-PU-FD method is at least

comparable but often better than the one of the RBF-FD scheme. Additionally, the proposed method offers a remarkable decrease of criterion of sparsity in coefficient matrix and CN compared to the global approach and show that we can attain high order convergence. Finally, we conclude this section by solving the stochastic pseudoparabolic problem on another domain Ω, i.e. the circular domain

Table 7 Comparison between RBF-PU-FD and RBF-FD [15] using M4 for the stochastic pseudo-parabolic equation on the circular domain: MAE, RMS, CPU time (in seconds), CN and sparsity for uniform points with 𝛿𝑡 = 0.001 and 𝑇 = 1. RBF-PU-FD N 283 621 1088 1645

M 5 12 16 25

RBF-FD

𝜖

MAE

0.05 0.35 0.10 0.10

1.73 × 10 2.05 × 10−2 2.61 × 10−2 2.48 × 10−2

RMS −2

7.40 × 10 8.91 × 10−3 1.19 × 10−2 1.22 × 10−2 −3

time

CN

sparsity

1.40 7.83 15.90 44.91

5.97 × 10 6.68 × 10+00 1.19 × 10+01 5.19 × 10+01 +00

4.30 × 10 2.94 × 10−1 2.75 × 10−1 1.65 × 10−1 −1

132

𝜖

MAE

0.25 0.35 0.15 0.15

2.87 × 10 5.93 × 10−2 2.99 × 10−2 4.17 × 10−2

RMS −2

1.16 × 10 2.31 × 10−2 1.30 × 10−2 1.66 × 10−2 −2

time

CN

sparsity

7.29 26.69 77.39 324.83

4.43 × 10 1.35 × 10+13 6.39 × 10+15 1.03 × 10+17 +12

1 1 1 1

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

Fig. 14. Graphs of approximate value of mean solution, standard deviation and absolute error obtained and for 1645 uniform points, 𝜖 = 0.10 (left) and 1709 Halton points, 𝜖 = 0.10 (right) with 25 patches contained in the circular domain, using M4 with 𝛿𝑡 = 0.001 and 𝑇 = 1.

depicted in Fig. 1 (bottom). Therefore, taking some sets of uniform points, in Tables 7 we compare our RBF-PU-FD scheme with the global RBF-FD method. Similar results are obtained for the Halton points that have been removed to be shortened. In particular, as an example, we here show the results obtained by using the M4 as local RBF approximant. Analyzing these experiments, we can thus observe a behavior – in terms of MAE, RMS, CPU time, CN and sparsity – similar to that exhibited for a square domain and already beforehand remarked. In addition, taking sets of uniform and Halton points and using M4, Figs. 13 and 14 shows the surface of the exact solution in the absence of Gaussian noise, absolute error, statistical criterions such as mean, standard deviation for given uniform (left) and Halton (right) points.

Gaussian field. The resulting Gaussian field is approximated at the collocation points at each time step. From this study, it derives that the RBF-PU-FD collocation method allows overcoming the high computational cost associated with the ill-conditioned and dense matrices of global RBF method, while maintaining high accuracy, since the matrices formed during the RBF-PU-FD method be sparse. In fact, the RBF-PU-FD technique also enables to reach a given level of accuracy with significantly less computational effort than the global RBF-FD method. The gained results show the capabilities and efficiency of the RBF-PU-FD method for time-dependent SPDEs with high dimensional and complicated geometries. Also, our method gives statistical criterions such as mean, standard deviation, lower bound and upper bound for prediction which are evaluated using the Monte-Carlo method. The fact that RBF methods are meshfree permits an easy implementation of adaptive grids, which can be clustered around critical regions such as the strike area or the free boundary, in order to improve accuracy or reduce overall computational cost. In the case of RBF-PU-FD method, refinements can be made independently within the partitions, increasing the flexibility. However, this topic is out of the scope of the present paper. Therefore, the development of a technique for automatic adaptivity will be part of our future work, where we will also consider higher-dimensional

6. Conclusions We proposed an RBF-PU-FD method used to solve time-dependent SPDEs such as the 2D stochastic convection-diffusion and stochastic pseudo-parabolic equations. The time evolution is discretized by using the 𝜃-weighted FD scheme, while the RBF-PU approximation procedure is employed for the space variables. After time discretization the Wiener process involves in the problem under consideration transforms into

133

M. Esmaeilbeigi, O. Chatrabgoun and M. Shafa

Engineering Analysis with Boundary Elements 104 (2019) 120–134

computational problems studying and analyzing theoretically the invertibility issue of the coefficient matrix of our scheme.

[22] Fasshauer GE, Ye Q. A kernel-based collocation method for elliptic partial differential equations with random coefficients. In: Dick J, Kuo F, Peters G, Sloan I, editors. Monte Carlo and Quasi-Monte Carlo Methods 2012. Springer Proceedings in Mathematics & Statistics, 65. Berlin, Heidelberg: Springer; 2013. [23] Fasshauer GE, Ye Q. Kernel-based collocation methods versus Galerkin finite element methods for approximating elliptic stochastic partial differential equations. In: Meshfree methods for partial differential equations VI, Lecture notes in computational science and engineering, 89; 2013. p. 155–70. [24] Heryudono A, Larsson E, Ramage A, Sydow LV. Preconditioning for radial basis function partition of unity methods. J Sci Comput 2016;67:1089–109. [25] Hon YC, Schaback R. On unsymmetric collocation by radial basis functions. Appl Math Comput 2001;119:177–86. [26] Kaminski M. The stochastic perturbation method for computational mechanics. Chichester: Wiley; 2013. [27] Kamrani M, Hosseini SM. The role of coefficients of a general SPDE on the stability and convergence of a finite difference method. J Comput Appl Math 2010;234:1426–34. [28] Kloeden PE, Platen E. Numerical solutions of stochastic differential equations. Springer; 1992. [29] Kruse R. Optimal error estimates of Galerkin finite element methods for stochastic partial differential equations with multiplicative noise. IMA J Numer Anal 2014;34(1):217–51. [30] Li X. Three-dimensional complex variable element-free Galerkin method. Appl Math Model 2018;63:148–71. [31] Lord GJ, Tambue A. Stochastic exponential integrators for the finite element discretization of SPDEs for multiplicative & additive noise. IMA J Numer Anal 2013;2:515–43. [32] Mavrič B, Šarler B. Local radial basis function collocation method for linear thermoelasticity in two dimensions. Int J Numer Methods Heat Fluid Flow 2015;25(6):1488–511. [33] Melenk JM, Babuška I. The partition of unity finite element method: Basic theory and applications. Comput Methods Appl Mech Eng 1996;139:289–314. [34] Robler A, Seaid M, Zahri M. Method of lines for stochastic boundary-value problems with additive noise. Appl Math Comput 2008;199:301–14. [35] Safdari-Vaighani A, Heryudono A, Larsson E. A radial basis function partition of unity collocation method for convection-diffusion equations. J Sci Comput 2015;64:341–67. [36] Shu C, Ding H, Yeo KS. Local radial basis function-based dierential quadrature method and its application to solve two-dimensional incompressible Navier–Stokes equations. Comput Methods Appl Mech Eng 2003;192:941–54. [37] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4:389–96. [38] Wendland H. Fast evaluation of radial basis functions: Methods based on partition of unity. In: Chui CK, Schumaker LL, J Stöckler e, editors. Approximation theory X: wavelets, splines, and applications. Nashville, TN: Vanderbilt University Press; 2002. p. 473–83. [39] Wendland H. Scattered data approximation, vol. 17. Cambridge: Cambridge Univ. Press; 2005. Cambridge Monographs on Applied and Computational Mathematics. [40] Yang X, Qi R. Nonconforming finite element method for stochastic stokes equations. Appl Math Model 2013;37:6110–18. [41] Ye Q. Analyzing reproducing kernel approximation methods via a green’s function approach. Illinois Institute of Technology; 2012. Ph.d. thesis. [42] Ye Q. Approximation of nonlinear stochastic partial differential equations by a kernel-based collocation method. Int J Appl Nonlinear Sci 2014;1:156–72. [43] Ye Q. Solving high-dimensional linear stochastic partial differential equations via a kernel-based approximation method. arXiv:1303.5381v5[math.NA]. [44] Zhang T, Li X. A generalized element-free Galerkin method for stokes problem. Comput Math Appl. 2018;75:3127–38.

References [1] Allen EJ. Modeling with ito stochastic differential equations. Dordrecht, The Netherlands: Springer; 2007. [2] Babuška I, Melenk JM. The partition of unity method. Int J Numer Methods Eng 1997;40:727–58. [3] Buhmann MD. Radial basis functions: Theory and implementation, vol. 12. Cambridge: Cambridge University Press; 2003. Cambridge Monographs on Applied and Computational Mathematics. [4] Cavoretto R, De Marchi S, De Rossi A, Perracchione E, Santin G. Partition of unity interpolation using stable kernel-based techniques. Appl Numer Math 2017;116:95–107. [5] Cavoretto R, De Rossi A, Perracchione E. Optimal selection of local approximants in RBF-PU interpolation. J Sci Comput 2018;74(1):1–22. [6] Cavoretto R, Schneider T, Zulian P. OpenCL based parallel algorithm for RBF-PUM interpolation. J Sci Comput 2018;74(1):267–89. [7] Garmanjani G, Cavoretto R, Esmaeilbeigi M. A RBF partition of unity collocation method based on finite difference for initial boundary value problems. Comput Math Appl 2018;75:4066–90. [8] Chantasiriwan S. Investigation of the use of radial basis functions in local collocation method for solving diusion problems. Int Commun Heat Mass Trans 2004;31(8):1095–104. [9] Chernov A, Schwab C. First order k-th moment finite element analysis of nonlinear operator equations with stochastic data. Math Comput 2013;82:1859–88. [10] Chow P-L. Explosive solutions of stochastic reaction-diffusion equations in mean lp -norm. J Differ Equ 2011;250:2567–80. [11] Chow P-L. Stochastic partial differential equations. Boca Raton, FL: Chapman and Hall/CRC; 2007. Chapman and Hall/CRC applied mathematics and nonlinear science series. [12] Cialenco I, Fasshauer GE, Ye Q. Approximation of stochastic partial differential equations by a kernel-based collocation method. Int J Comput Math 2012;89:2543–61. [13] Dehghan M, Abbaszadeh M, Mohebbi A. A meshless technique based on the local radial basis functions collocation method for solving parabolic-parabolic patlak-keller-segel chemotaxis model. Eng Anal Bound Elem 2015;56:129–44. [14] Dehghan M, Mohammadi V. The numerical solution of cahnehilliard (CH) equation in one, two and three dimensions via globally radial basis functions (GRBFs) and RBFs-dierential quadrature (RBFs-DQ) methods. Eng Anal Bound Elem 2015;51:74–100. [15] Dehghan M, Shirzadi M. Meshless simulation of stochastic advection-diffusion equations based on radial basis functions. Eng Anal Bound Elem 2015;53:18–26. [16] Dehghan M, Shirzadi M. Numerical solution of stochastic elliptic partial differential equations using the meshless method of radial basis functions. Eng Anal Bound Elem 2015;50:291–303. [17] Dehghan M, Shirzadi M. The modified dual reciprocity boundary elements method and its application for solving stochastic partial differential equations. Eng Anal Bound Elem 2015;58:98–111. [18] Dehghan M, Shokri A. A numerical method for two-dimentional Schrödinger equation using collocation radial basis functions. Comput Math Appl 2007;54:136–46. [19] Etheridge A. Stochastic partial differential equations. Cambridge, UK: Cambridge University Press; 1995. [20] Fasshauer GE. Meshfree approximation methods with matlab. Singapore: World Scientific; 2007. [21] Fasshauer GE, McCourt MJ. Kernel-based approximation methods using matlab. Singapore: World Scientific; 2015.

134