Poly-Sinc Solution of Stochastic Elliptic Differential Equations

TL;DR

This paper presents a novel numerical method combining polynomial chaos and Sinc point approximation to efficiently solve stochastic elliptic PDEs with high accuracy using few collocation points.

Contribution

It introduces a new approach integrating polynomial chaos with Sinc point interpolation for solving high-dimensional stochastic PDEs efficiently.

Findings

Accurate solutions achieved with few collocation points.

Method effectively handles high-dimensional stochastic systems.

Numerical examples demonstrate high accuracy with Legendre polynomials.

Abstract

In this paper, we introduce a numerical solution of a stochastic partial differential equation (SPDE) of elliptic type using polynomial chaos along side with polynomial approximation at Sinc points. These Sinc points are defined by a conformal map and when mixed with the polynomial interpolation, it yields an accurate approximation. The first step to solve SPDE is to use stochastic Galerkin method in conjunction with polynomial chaos, which implies a system of deterministic partial differential equations to be solved. The main difficulty is the higher dimensionality of the resulting system of partial differential equations. The idea here is to solve this system using a small number of collocation points. Two examples are presented, mainly using Legendre polynomials for stochastic variables. These examples illustrate that we require to sample at few points to get a representation of a model that is sufficiently accurate.

Equations

$$\begin{array}{rl}\mathcal{L}(u)\equiv -\nabla \cdot (a(x,y,\Theta )\nabla u(x,y,\Theta ))& =f(x,y)\text{ in }Q\times \Omega \text{ and }\\ u& =0\text{ on }\partial Q\times \Omega ,\end{array}$$

$$a(x,y,\Theta )=a_0(x,y)+b_0\sum _{k=1}^K{\xi }_k{a}_k(x,y),$$

$$u(x,y,\Theta )=\sum _{i=0}^{\infty }{u}_i(x,y){\Phi }_i(\Theta ).$$

$$u(x,y,\Theta )\simeq \tilde{u}(x,y,\Theta )=\sum _{i=0}^m{u}_i(x,y){\Phi }_i(\Theta )$$

$$u_i(x,y)=⟨u(x,y,\Theta ),{\Phi }_i(\Theta )⟩,i=0,1,\dots ,m.$$

$$⟨{\Phi }_i(\Theta ),{\Phi }_j(\Theta )⟩={\int }_{I^K}{\Phi }_i(\Theta ){\Phi }_j(\Theta )W(\Theta )d\Theta =c_i{\delta }_{ij},\text{for all }i,j,$$

$${\Phi }_i(\Theta )={\Phi }_{\mathbf{i}}(\Theta ):=\prod _{r=1}^K{\Phi }_{i_r}^{(r)}({\xi }_r),\mathbf{i}=(i_1,\dots ,i_K),$$

$${\mathcal{R}}_P=\left\{{\Phi }_i(\Theta ):\sum _{r=1}^K{i}_r\le P\right\},$$

$$m+1=\frac{(K+P)!}{K!P!}.$$

$$⟨{\Phi }_i(\Theta ),{\Phi }_i(\Theta )⟩=1.$$

$$E(u(x,y,\Theta ))$$

$$E\left(u(x,y,\Theta )\right)=u_0(x,y).$$

$$V(u(x,y,\Theta ))$$

$$\displaystyle=\leavevmode\resizebox{375.80542pt}{}{$\int_{I^{K}}[u^{2}(x,y,\Theta)+E^{2}(u(x,y,\Theta))-2u(x,y,\Theta)E(u(x,y,\Theta))]\,W(\Theta)\,d\Theta$}$$

$$={\int }_{I^K}[u^2(x,y,\Theta )+u_0^2(x,y)-2u_0(x,y)u(x,y,\Theta )]W(\Theta )d\Theta$$

$$=⟨u(x,y,\Theta ),u(x,y,\Theta )⟩-u_0^2(x,y).$$

$$V(u(x,y,\Theta ))\approx \sum _{i=1}^m{u}_i^2(x,y).$$

$$⟨\mathcal{L}(\tilde{u})-f(x,y),{\Phi }_j(\Theta )⟩=0\text{for }j=0,1,\dots ,m.$$

$$\mathcal{L}(\tilde{u})=-\nabla \cdot \left(\nabla \sum _{i=0}^m{u}_i(x,y){\Phi }_i(\Theta )\right)-\frac{1}{2}\sum _{k=1}^K{\xi }_k\nabla \cdot \left(a_k\nabla \sum _{i=0}^m{u}_i(x,y){\Phi }_i(\Theta )\right).$$

$$-{\nabla }^2{u}_j-\frac{1}{2}\sum _{k=1}^K\sum _{i=0}^m⟨{\xi }_k{\Phi }_i(\Theta ),{\Phi }_j(\Theta )⟩\nabla \cdot \left(a_k\nabla u_i\right)=F_j$$

$$P_n(x_k)=u_k,k=-M,...,N.$$

$$P_n(x)=\sum _{k=-M}^N{b}_k(x)u(x_k),$$

$$b_k(x)=\frac{g(x)}{(x-x_k)g^{\prime }(x_k)},g(x)=\prod _{j=-M}^N\left(x-x_j\right).$$

$$x_k=\frac{a+b{\mathrm{e}}^{kh}}{1+{\mathrm{e}}^{kh}}.$$

$$\mathbf{B}(x)$$

$$V_{n}u$$

$$u(x)\simeq \mathbf{B}(x)V_{n}u.$$

$$\underset{x\in (a,b)}{\sup }|u(x)-\mathbf{B}(x)V_{n}u|\le A\frac{\sqrt{N}}{B^{2N}}\exp \left(\frac{-{\pi }^2{N}^{\frac{1}{2}}}{2}\right),$$

$${\Lambda }_n\approx \frac{1}{\pi }\log (n+1)+\mathrm{1.07618.}$$

$$(P_{n}u)(X)=\sum _{k_1=-M_1}^{N_1}\sum _{k_2=-M_2}^{N_2}\dots \sum _{k_l=-M_l}^{N_l}u(X_{\mathbf{k}})b_{k_1}(x_1)b_{k_2}(x_2)\dots b_{k_l}(x_l),$$

$$u(X)\simeq \underset{i=1}{\overset{l}{⨀}}{\mathbf{B}}_i(X)U,i=1,2,\dots ,l.$$

$$E_n=\underset{X\in Q}{\sup }|u(X)-\underset{i=1}{\overset{l}{⨀}}{\mathbf{B}}_i(X)U|\le \sum _{i=0}^{l-1}\left(C_i{\log }^{i}N\right)\frac{\sqrt{N}}{{\gamma }_i^{2N}}\exp \left(\frac{-{\pi }^2{N}^{\frac{1}{2}}}{2}\right),$$

Full text

Keywords: Poly-Sinc methods, Collocation method, Galerkin method, Stochastic Differential Equations, Polynomial Chaos, Legendre Polynomials.

MSC Classification: 65N35, 65N12, 65N30, 65C20, 35R60.

1 Introduction

In many applications the values of the parameters of the problem are not exactly known. These uncertainties inherent in the model yield uncertainties in the results of numerical simulations. Stochastic methods are one way to model these uncertainties and shall model this by random fields [1]. If the physical system is described by a partial differential equation (PDE), then the combination with the stochastic model results in a stochastic partial differential equation (SPDE). The solution of the SPDE is again a random field, describing both the expected response and quantifying its uncertainty. SPDEs can be interpreted mathematically in several ways.

In the numerical framework, the stochastic regularity of the solution determines the convergence rate of numerical approximations, and a variational theory for this was earlier devised in [2]. The ultimate goal in the solution of SPDEs is usually the computation of response statistics, i.e. a functional of the solution. Monte Carlo (MC) methods can be used directly for this purpose, but they require a high computational effort [3, 5]. Quasi Monte Carlo (QMC) and variance reduction techniques [3] may reduce the computational effort considerably without requiring much regularity. However, often we have high regularity in the stochastic variables, and this is not exploited by QMC methods.

Alternatives to Monte Carlo methods have been developed. For example, perturbation methods [4], methods based on Neumann-series [6], or the spectral stochastic finite element method (SSFEM) [7, 9]. Stochastic Galerkin methods have been applied to various linear problems, see [7, 8, 11]. Nonlinear problems with stochastic loads have been tackled in [10]. These Galerkin methods yield an explicit functional relationship between the independent random variables and the solution. In contrast with common MC methods, subsequent evaluations of functional statistics like the mean and covariance are very cheap.

We consider an elliptic PDE in space including a random field as material parameters. The polynomial chaos approach and the stochastic Galerkin method yield a deterministic system of PDEs in space [14]. In this paper, we introduce a spatial collocation technique based on polynomial approximation by Lagrange interpolation. For the interpolation points we use a specific set of non-uniform points created by conformal maps, called Sinc points. Later, we use a small number of Sinc points as collocation points to compute a very accurate solution of the PDEs, see [23].

The paper is organized as follows: In Section 2, we introduce a model problem, the structure of its polynomial chaos model and the stochastic Galerkin solution. In Section 3, we illustrate the main theorem of Poly-Sinc approximation. In Section 4, we review a Poly-Sinc collocation technique with the main collocation theorem. Finally, in Section 5, we investigate numerical examples. We start with a simple example in one stochastic variable and then we discuss the general model from Section 2.

2 Stochastic Model Problem

In this paper, we are interested to solve the following stochastic partial differential equations:

[TABLE]

where $\Theta=\left(\xi_{1},\xi_{2},\text{...},\xi_{K}\right)$ is a vector of stochastic parameters. These parameters are independent and uniformly distributed in $I=[-1,1]$ and thus $\Theta:\Omega\longrightarrow[-1,1]^{K}$ with an event space $\Omega$. Moreover the domain of the spatial variables $x$ and $y$ is $Q=(0,1)^{2}$. The function $a(x,y,\Theta)$ is defined as

[TABLE]

where $a_{k}$’s are functions in $x$ and $y$ only, $b_{0}$ is a constant and, $\xi_{k}$’s are the random variables. Without loss of generality, we consider $a_{0}=1$ and $b_{0}=1/2$. We assume that $a(x,y,\Theta)\geq\alpha>0$ for all $(x,y)\in Q$ and all $\Theta\in[-1,1]^{K}$. Thus the differential operator in (1) is always uniformly elliptic.

In the rest of the section, we introduce the main concepts used in the solution of (1) with (2). Basically, we discuss the polynomial chaos in one- and multidimensional cases and the stochastic Galerkin method.

2.1 Polynomial Chaos Expansion

Generalized Polynomial Chaos (gPC) is a particular set of polynomials in a given random variable, with which an approximation of a finite second-moment random variable is computed. This procedure is named Polynomial Chaos Expansion (PCE). This technique exploits orthogonal properties of polynomials involved, to detect a representation of random variables as series of functionals. Now, the function $u$ can be expressed as an infinite series of orthogonal basis functions $\Phi_{i}$ with suitable coefficient functions $u_{i}$ as

[TABLE]

The expansion in (3) converges in the mean square of the probability space. The truncation form including $m+1$ basis functions leads to

[TABLE]

with coefficients functions

[TABLE]

A fundamental property of the basis functions is the orthogonality,

[TABLE]

where $c_{i}$ are real positive numbers and $\delta_{ij}$ is the Kronecker-delta. In general, the inner product in (5) can be defined for different types of weighting function $W$; however, it is possible to prove that the optimal convergence rate of a gPC model can be achieved when the weighting function $W$ agrees to the joint probability density function (PDF) of the random variables considered in a standard form [8, 12]. In this framework, an optimal convergence rate means that a small number of basis functions is sufficient to obtain an accurate PC model (4). Hence, the choice of the basis functions depends only on the probability distribution of the random variables $\Theta$, and it is not influenced by the type of system under study. In particular, if the random variables $\Theta$ are independent, their joint PDF corresponds to the product of the PDFs of each random variable: in this case, the corresponding basis functions $\Phi_{i}$ can be calculated as product combinations (tensor product) of the orthogonal polynomials corresponding to each individual random variable [13, 14, 15]:

[TABLE]

where $\Phi_{i_{r}}^{(r)}$ represents the univariate basis polynomial of degree $i_{r}$ associated to the $r$th random parameter and with one-to-one correspondence between the integers $i$ and the multi-indices $\mathbf{i}$. We assume ${\rm{degree}}(\Phi_{i})\leq{\rm{degree}}(\Phi_{i+1})$ for each $i$. Now let

[TABLE]

be the set of all multivariate polynomials up to total degree $P$ as used in a Taylor expansion. Furthermore, for random variables with specific PDFs, the optimal basis functions are known and are formed by the polynomials of the Wiener-Askey scheme [8]. For example, in the uniform probability distribution, the basis functions are the Legendre polynomials.

Using (6) and (7), it is possible to show that the total number of basis functions $m+1$ in (4) is expressed as

[TABLE]

The total degree of the PC (the maximum degree) $P$ can be chosen relatively small to achieve the desired accuracy in the solution.

In the case of the orthogonal polynomials, we can see that $\Phi_{0}(\Theta)=1$ and for orthonormal polynomials

[TABLE]

Once a PC model in the form of (4) is obtained, stochastic moments like the mean $E(u)$ and the variance $V(u)$ can be analytically calculated by the PC expansion coefficients as, see [15],

[TABLE]

Using the PC expansion of $u(x,y,\Theta)$ given in (3) and the orthogonality of the basis functions $\Phi_{i}(\Theta)$ to get

[TABLE]

The variance can be derived by

[TABLE]

Again, use the PC expansion in (4) and orthonormal polynomials basis satisfying (9), to get

[TABLE]

It is clear now that, in order to obtain a PC model in (4) and the stochastic moments, the coefficients functions $u_{i}(x,y)$ must be computed. The PC coefficient estimation depends on the type of the resulting system from the chaos expansion, not only the PC truncation.

2.2 Stochastic Galerkin Method

To solve the problem in (1) and (2), a Galerkin method is used along side the PC. The main idea is to assume that the solution of (1) and (2) is written as expansion in (4) and then use the PC theory introduced in the previous section. This process transform the SPDE (1) and (2) into a deterministic system of PDEs.

To recover the coefficient functions $u_{i}(x,y)$ we apply the inner product of (1) with the basis polynomial $\Phi_{j}(\Theta)$

[TABLE]

Substituting (4) in (1) we obtain

[TABLE]

Now applying the inner residual product in (10) and use the orthogonality property of the multivariate basis $\Phi_{i}$’s to get

[TABLE]

where $F_{j}(x,y)=\left\langle f(x,y),\Phi_{j}(\Theta)\right\rangle$ forms an $(m+1)$ vector and the array $\left\langle\xi_{k}\Phi_{i},\Phi_{j}\right\rangle$ is a triple tensor of dimension $K\times(m+1)\times(m+1)$. (11) is a system of elliptic PDEs with unknown variables $u_{i}(x,y)$, $i=0,1,2,\ldots,m$. With large number of random variables $K$ (say $K>4$) the size of the system in (11) becomes huge due to (8). One of our targets in the solution of the system in (11) is to use a collocation method to achieve a high accuracy with small numbers of collocation points. The proposed method in this report is to use Sinc points in a Lagrange interpolation.

2.3 Quadrature

The inner product $\left\langle.,.\right\rangle$ is defined by an integral. For the integration of polynomials analytic methods are used. Alternatively, we can use highly accurate quadrature techniques to evaluate the integrals exactly except for round-off errors. We omit the details of these techniques, since they can be easily found in several textbooks. For example, descriptions of Gaussian quadrature can be found in most texts on numerical analysis [18], while [16] contains descriptions of Sinc quadratures over finite, semi-infinite, infinite intervals and contours.

3 Poly-Sinc Approximation

In this section we introduce the Lagrange approximation at Sinc points as interpolation points. This approximation is called Poly-Sinc approximation [23]. It was first introduced to provide a uniform approximation for a function and its derivatives as well [17]. In [23], the main results of this approximation have been extended and have been used to solve differential equations.

Given a set of data $\left\{x_{k},u(x_{k})\right\}_{k=-M}^{N}$ where the $x_{k}$ are interpolation points on $(a,b)$. Then there is a unique polynomial $P_{n}(x),\,n=M+N+1$ of degree at most $n-1$ satisfying the interpolation condition,

[TABLE]

In this case $P_{n}(x)$ can be expressed by the Lagrange polynomials as

[TABLE]

with,

[TABLE]

Now $x_{k}$’s are Sinc points on $(a,b)$ defined as [16]

[TABLE]

Corresponding to such a scheme, we define a row vector $\bm{B}$ of basis functions and an operator $V_{m}u$ that maps a function $u(x)$ into a column vector of dimension $n=M+N+1$ by

[TABLE]

This notation enables us to write the above interpolation scheme in simple operator form, as

[TABLE]

This approximation, like regular Sinc approximation, yields an exceptional accuracy in approximating the function that is known at Sinc points, [16]. Unlike Sinc approximation, it gives a uniform exponential convergence rate when differentiating the interpolation formula given in (13), see [17]. Next, we assume that $M=N$ and that $n=2N+1$ is the total number of Sinc points. Then the upper bound of error for Poly-Sinc approximation is given as

[TABLE]

where $A>0$ and $B>1$ are two constants, independent of $N$. For the proof of (14), see [17].

Another criterion to discuss the convergence and stability of the Poly-Sinc approximation is the Lebesgue constant. In [21] an estimate for the Lebesgue constant for Lagrange approximation at Sinc points (12) has been derived as

[TABLE]

Next we extend these results from the one-dimensional case to the multi-dimensional case.

Let $X=(x_{1},.....,x_{l})$ be a point in $Q=[a,b]^{l}$, then Lagrange approximation of a function $u(X)$ can be defined by a nested operator as

[TABLE]

where $u(X_{\boldsymbol{k}})=U=u(x_{1,k_{1}},\ldots,x_{l,k_{l}})$ with $k_{i}=-M_{i},\ldots,N_{i}$. We can write the approximation (15) in the operator form

[TABLE]

Next, we assume $M_{i}=N_{j}=N,\,i,j=1,\ldots,l$ and $n=2N+1$ is the number of Sinc points in each dimension $i=1,2,\ldots,l$. The convergence and stability of the approximation (16) are discussed in [19] and [21]. For the upper bound of the error $E_{n}$, we have

[TABLE]

where $C_{i}>0$, $\gamma_{i}>1$, $i=1,\ldots,l$ are two sets of constants, independent of $N$.

The notation $\Lambda_{n,l}$ is used to denote the Lebesgue constant using $n$ interpolation points in each dimension $i=1,2,\ldots,l$, i.e. $n^{l}$ Sinc points in total. If $P_{n}(X)$ is defined as in (15), then:

[TABLE]

4 Poly-Sinc Collocation Method

In [20], a collocation method based on the use of bivariate Poly-Sinc interpolation defined in (16) is introduced to solve elliptic equations defined on rectangular domains. In [22], Poly-Sinc collocation domain decomposition method for elliptic boundary value problems is investigated on complicated domains. The idea of the collocation method is to reduce the boundary value problem to a system of algebraic equations which have to be solved subsequently. To start let us introduce the following collocation theorem.

Theorem 1**.**

Let $u:\overline{Q}\rightarrow\mathbb{R}$ be an analytic bounded function on the compact domain $\overline{Q}$. Let $U=\left\{u(x_{j},y_{k})\right\}^{N}_{j,k=-N}$ be a vector, where $x_{j}$ and $y_{k}$ are the Sinc points. If ${\widetilde{U}}=\left\{{\widetilde{u}_{jk}}\right\}^{N}_{j,k=-N}$ is a vector satisfying

[TABLE]

then

[TABLE]

where $n=2N+1$, $E_{n}$ from (17), and $\Lambda_{n,2}$ from (18).

Proof**.**

We apply triangle inequality

[TABLE]

which is the statement of the theorem. $\hfill\square$

This theorem guarantees an accurate final approximation of $u$ on its domain of definition provided that we know a good approximation to $u$ at the Sinc points.

To set up the collocation scheme, let us consider the following partial differential operator,

[TABLE]

where $Q=\left\{a<x<b,\,c<y<d\right\}$ and $u_{x\,x}=\frac{\partial^{2}u}{\partial{x^{2}}}$, $u_{y\,y}=\frac{\partial^{2}u}{\partial{y^{2}}}$.

The first step in the collocation algorithm is to replace $u(x,y)$ in Eq. (20) by the Poly-Sinc approximation defined in (16). Next, we collocate the equation by replacing $x$ and $y$ by Sinc points

[TABLE]

and

[TABLE]

In this case, we have,

[TABLE]

where,

[TABLE]

and $B^{{}^{\prime\prime}}(j,h)(x_{i})$ defines an $n\times n$ matrix, with $n=M+N+1$

[TABLE]

So,

[TABLE]

where ${\cal{M}}_{1}$ is a $n^{2}\times n^{2}$ matrix defined as,

[TABLE]

and where ${\cal{U}}_{x\,x}$ is collected in a vector of of length $n^{2}$. Likewise, it holds that

[TABLE]

where ${\cal{M}}_{2}$ is defined in the same way as ${\cal{M}}_{1}$.

The differential equation has been transformed to a system of $n^{2}$ algebraic equations,

[TABLE]

where $\mathcal{U}$ is the vector of length $n^{2}$ including the unknowns $u_{i\,q}$ and

[TABLE]

The right hand side $\mathcal{F}$ is a vector of Length $n^{2}$ and defined as

[TABLE]

Now the PDE (20) is transformed to a system of $n^{2}$ algebraic equations in $n^{2}$ unknowns. The boundary conditions are collocated separately to yield $4n$ algebraic equations. More precisely,

[TABLE]

where $x_{i}$ and $y_{j}$ are the Sinc points defined on $(a,b)$ and $(c,d)$, respectively. Adding these $4n$ equations to the $n^{2}\times n^{2}$ algebraic system, produced from the collocation of the PDE, yields a rectangular system of linear equations. Finally, solving this least squares problem yields the desired numerical solution.

Note:

• •

In our calculations we used a multiplier factor $\tau=10^{3}$ in the collocation steps of the homogenous boundary conditions. This factor emphasizes the boundary values and improve the error behavior at the boundaries.

• •

The Poly-Sinc collocation technique is based on the collocation of the spatial variables using Sinc points. This means that it is valid also for PDEs with space-dependent coefficients. Moreover, it can be generalized to solve a system of PDEs.

5 Numerical Results

In this section, we present the computational results. Mainly, we discuss two examples. The first simple example includes one stochastic parameter. In the second example we solve the model problem introduced in Section 2.

5.1 One Stochastic Variable

Consider the Poisson equation in two spatial dimensions with one random parameter. This problem is described by the following SPDE

[TABLE]

where $Q=(-1,1)^{2}$ is the spatial domain and $\Omega$ is an event space and $\xi:\Omega\rightarrow[-1,1]$ is a random variable. The function $a(\xi)=\xi+2$ is a linear function of a uniformly distributed random variable $\xi$ and $f(x,y)=1$ for all $(x,y)\in Q$.

Now, we use the PC representation in (4) with $m=3$ to have

[TABLE]

where $\Phi_{i}$’s are the univariate orthonormal Legendre polynomials defined on $[-1,1]$. Substitution of (22) in the SPDE (21) yields the residual

[TABLE]

We then perform a Galerkin projection and use the orthogonality of Legendre polynomials, which yields the system of elliptic PDEs

[TABLE]

where $\mathcal{L}{u_{i}}=(u_{i})_{xx}+(u_{i})_{yy}$. It holds that $\left\langle 1,\Phi_{k}\right\rangle=\delta_{1k}$.

The computational results of this example are given in the following experiments.

Experiment 1**.**

$\mathbf{E(u)}$* and $\mathbf{V(u)}$

In this experiment, we use Poly-Sinc collocation from Section 4 to solve the system of PDEs in (23). In our computation, we use $N=5$, i.e. $11\times 11$ of 2D grid of Sinc points defined as in (12) on the domain $Q$. As a result of the Poly-Sinc solution, the coefficient functions $u_{i}(x,y)$ are obtained. In Fig. 1, the expectation $E(u)=u_{0}(x,y)$ and its contour plot are represented while in Fig. 2, the variance calculations are presented.*

Experiment 2**.**

Coefficients functions*

As we mentioned above, to get an accurate result, just a small number of orthogonal polynomials, $\Phi_{i}$, is needed. In our computations, we used $m=3$, i.e. four orthonormal Legendre polynomials. The $4$ coefficients functions, $u_{i}(x,y),\,i=0,\ldots,3$, are given in Fig.3. In addition, we verify that this number is sufficient by showing that the coefficient functions $u_{i}$ tend to zero as $m$ increases. The results are given in Fig.4. In Fig.4, the dots represent the maximum of the coefficient functions $u_{i}(x,y)$ on the spatial domain. We then use these maximum values in a least square estimation to find the coefficients of the decaying rate function $\alpha\,\exp(-\beta s)$, where $\alpha$ and $\beta$ are constants. In Fig.4, the solid line represents the best fitting function with $\alpha=0.14$ and $\beta=1.2$. This means that the coefficient functions $u_{i}(x,y)$ follow an exponentially decay relation.*

Experiment 3**.**

Error*

To discuss the convergence of Poly-Sinc solution, we need a reference (nearly exact) solution. For that, we create a discrete list of PDEs of the equation (21) at a finite set of instances of $\xi\in[-1,1]$. We choose $100$ points of Gauss-Legendre nodes as values of $\xi\in[-1,1]$ and create corresponding $100$ PDEs. To solve each one of these $100$ equations we use Mathematica Package NDSolve. NDSolve uses a combination of highly accurate numeric schemes to solve initial and boundary value PDEs 111For more information about NDSolve, see Wolfram documentation center at https://reference.wolfram.com/language/ref/NDSolve.html. We then calculate the expectation and variance of the solutions of our set of boundary value problems of PDEs. In Fig.5, the errors in the calculations of $E(u)$ and $V(u)$ using $m=3$ and Poly-Sinc (with orthonormal Legendre) and the references from the $100$ PDEs are presented. Using the spatial $L_{2}$-norm error, calculating the error in both $E(u)$ and $V(u)$ delivers error of order $\mathcal{O}(10^{-4})$ and $\mathcal{O}(10^{-6})$, respectively. In Fig.6, the error between the solution of the SPDE in (21), using the method in this paper, and the reference solution is presented. We choose four instances of $\xi$.*

Experiment 4**.**

Comparison*

In this experiment we compare the Poly-Sinc solution with the classical finite difference (FD) solution. In 5-point-star FD method [24], we use an $11\times 11$ meshing with constant step size for the spatial variables $x$ and $y$, which is the same number of Sinc points used in the Poly-Sinc solution. The error between finite difference solution and the reference exact solution is given in Fig.7. Using the spatial $L_{2}$-norm error, calculating the error in both $E(u)$ and $V(u)$ delivering error of order $\mathcal{O}(10^{-2})$. These calculations shows that for the same number of points, Poly-Sinc delivers better approximation for the solution of the SPDE. In Fig.8 we run the calculations for different numbers of Sinc points $n=2N+1$ and use the same number of points in the FD method. We then calculate the $L_{2}$-norm error. Fig.8 shows that the decaying rate of the error, in both mean and variance, is better in Poly-Sinc than the FD method. Moreover, the Poly-Sinc decaying rates of errors are following qualitatively the upper bound in formula (19).*

5.2 Multiple Stochastic Variables

We solve the model problem defined in Section 2 for five stochastic variables, cf. [25]. Consider the SPDE defined in (1) with $K=5$ in (2) and where,

[TABLE]

$\Theta=\left\{\xi_{k}\right\}^{5}_{k=1}$ is a set of independent random variables uniformly distributed in $[-1,1]$. For this SPDE we run four experiments.

Experiment 5**.**

$\mathbf{E(u)}$* and $\mathbf{V(u)}$

In this experiment, we perform the Galerkin method along side the multivariate PC. For the PC parameters, we choose $K=5$ and $P=3$. Due to (8), the number of multivariate Legendre polynomials is $m+1=56$. As a result the three-dimensional array $\left\langle\xi_{k}\Phi_{i}(\Theta),\Phi_{j}(\Theta)\right\rangle$ is of dimension $5\times 56\times 56$. For the Poly-Sinc solution of the resulting system of PDEs, we use $N=5$, i.e. $n=11$ Sinc points. In Fig. 9 and Fig.10 the expectation and variance plots are presented.*

Experiment 6**.**

Coefficients Functions*

Similar to the second experiment in Example 1, we would like to study the accuracy of the polynomial expansion. In other words, study the decaying rate, to zero, of these functions. In Fig.12, the first six coefficients functions of the Poly-Sinc solution are given. These six coefficient functions are associated to the basis polynomials of degree zero and one. In Fig. 12, the logarithmic plot of the maximum of the absolute value of the coefficient functions $u_{i-1}(x,y),\,\,i=1,\ldots,56$ on the spatial domain is presented. We can see the fast decaying rate to zero.*

Experiment 7**.**

Error*

The idea of creating a set of (exact) instance solutions we used in the previous example is not applicable here as we have a set of $5$ random variables. For that we need to find a different reference to check the accuracy of our solution. We use the Finite Element (FE) solution with cell meshing $10^{-3}$ to solve the resulting system of PDEs. The FE element method is a part of the package NDSolve”FEM” in Mathematica 11 that uses the rectangular meshing of the domain and Dirichlet boundary conditions 222For more information about NDSolve ”FEM”, see Wolfram documentation center at https://reference.wolfram.com/language/FEMDocumentation/guide/FiniteElementMethodGuide.html. In Fig.13, the error for the expectation and variance is presented. Using the $L_{2}$-norm error, calculating the error in both $E(u)$ and $V(u)$ deliver error of order $\mathcal{O}(10^{-4})$ and $\mathcal{O}(10^{-8})$, respectively.*

Experiment 8**.**

Comparison*

In this experiment we compare the Poly-Sinc solution with the 5-point-star FD method. The reference solution is the Finite Element (FE) solution with cell meshing $10^{-3}$. In Fig.14 we run the calculations for different numbers of Sinc points $n=2N+1$ and use the same number of points in FD. We then calculate the $L_{2}$-norm error. These calculations show that the decaying rate of the error, in both mean and variance, is better in Poly-Sinc than the FD method. Moreover, the Poly-Sinc decaying rates of errors are following qualitatively the exponential decaying rate in (19).*

6 Conclusion

In this work we have formulated an efficient and accurate collocation scheme for solving a system of elliptic PDEs resulting from an SPDE. The idea of the scheme is to use a small number of collocation points to solve a large system of PDEs. We introduced the collocation theorem based on the error rate and the Lebesgue constant of the 2D Poly-Sinc approximation. As applications, we discussed two examples, the first example with one random variable while the other with five random variables. For each case the expectation, variance, and error are discussed. The experiments show that using Poly-Sinc approximation to solve the system of PDEs is an efficient method. The number of Sinc points needed to get this accuracy is small and the error decays faster than in the classical techniques, as the finite difference method.