PolySinc Solution of Stochastic Elliptic Differential Equations
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.
Keywords: PolySinc 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 Neumannseries [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 nonuniform 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 PolySinc approximation. In Section 4, we review a PolySinc 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:
(1) 
where is a vector of stochastic parameters. These parameters are independent and uniformly distributed in and thus with an event space . Moreover the domain of the spatial variables and is . The function is defined as
(2) 
where ’s are functions in and only, is a constant and, ’s are the random variables. Without loss of generality, we consider and . We assume that for all and all . 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 secondmoment 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 can be expressed as an infinite series of orthogonal basis functions with suitable coefficient functions as
(3) 
The expansion in (3) converges in the mean square of the probability space. The truncation form including basis functions leads to
(4) 
with coefficients functions
A fundamental property of the basis functions is the orthogonality,
(5) 
where are real positive numbers and is the Kroneckerdelta. In general, the inner product in (5) can be defined for different types of weighting function ; however, it is possible to prove that the optimal convergence rate of a gPC model can be achieved when the weighting function 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 , and it is not influenced by the type of system under study. In particular, if the random variables are independent, their joint PDF corresponds to the product of the PDFs of each random variable: in this case, the corresponding basis functions can be calculated as product combinations (tensor product) of the orthogonal polynomials corresponding to each individual random variable [13, 14, 15]:
(6) 
where represents the univariate basis polynomial of degree associated to the th random parameter and with onetoone correspondence between the integers and the multiindices . We assume for each . Now let
(7) 
be the set of all multivariate polynomials up to total degree 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 WienerAskey 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 in (4) is expressed as
(8) 
The total degree of the PC (the maximum degree) can be chosen relatively small to achieve the desired accuracy in the solution.
In the case of the orthogonal polynomials, we can see that and for orthonormal polynomials
(9) 
Once a PC model in the form of (4) is obtained, stochastic moments like the mean and the variance can be analytically calculated by the PC expansion coefficients as, see [15],
Using the PC expansion of given in (3) and the orthogonality of the basis functions to get
The variance can be derived by
It is clear now that, in order to obtain a PC model in (4) and the stochastic moments, the coefficients functions 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 we apply the inner product of (1) with the basis polynomial
(10) 
Now applying the inner residual product in (10) and use the orthogonality property of the multivariate basis ’s to get
(11) 
where forms an vector and the array is a triple tensor of dimension . (11) is a system of elliptic PDEs with unknown variables , . With large number of random variables (say ) 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 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 roundoff 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, semiinfinite, infinite intervals and contours.
3 PolySinc Approximation
In this section we introduce the Lagrange approximation at Sinc points as interpolation points. This approximation is called PolySinc 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 where the are interpolation points on . Then there is a unique polynomial of degree at most satisfying the interpolation condition,
In this case can be expressed by the Lagrange polynomials as
with,
Now ’s are Sinc points on defined as [16]
(12) 
Corresponding to such a scheme, we define a row vector of basis functions and an operator that maps a function into a column vector of dimension by
This notation enables us to write the above interpolation scheme in simple operator form, as
(13) 
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 and that is the total number of Sinc points. Then the upper bound of error for PolySinc approximation is given as
(14) 
Another criterion to discuss the convergence and stability of the PolySinc approximation is the Lebesgue constant. In [21] an estimate for the Lebesgue constant for Lagrange approximation at Sinc points (12) has been derived as
Next we extend these results from the onedimensional case to the multidimensional case.
Let be a point in , then Lagrange approximation of a function can be defined by a nested operator as
(15) 
where with . We can write the approximation (15) in the operator form
(16) 
Next, we assume and is the number of Sinc points in each dimension . The convergence and stability of the approximation (16) are discussed in [19] and [21]. For the upper bound of the error , we have
(17) 
where , , are two sets of constants, independent of .
The notation is used to denote the Lebesgue constant using interpolation points in each dimension , i.e. Sinc points in total. If is defined as in (15), then:
(18) 
4 PolySinc Collocation Method
In [20], a collocation method based on the use of bivariate PolySinc interpolation defined in (16) is introduced to solve elliptic equations defined on rectangular domains. In [22], PolySinc 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 be an analytic bounded function on the compact domain . Let be a vector, where and are the Sinc points. If is a vector satisfying
Proof.
We apply triangle inequality
which is the statement of the theorem.
This theorem guarantees an accurate final approximation of on its domain of definition provided that we know a good approximation to at the Sinc points.
To set up the collocation scheme, let us consider the following partial differential operator,
(20)  
where and , .
The first step in the collocation algorithm is to replace in Eq. (20) by the PolySinc approximation defined in (16). Next, we collocate the equation by replacing and by Sinc points
and
In this case, we have,
where,
and defines an matrix, with
So,
where is a matrix defined as,
and where is collected in a vector of of length . Likewise, it holds that
where is defined in the same way as .
The differential equation has been transformed to a system of algebraic equations,
where is the vector of length including the unknowns and
The right hand side is a vector of Length and defined as
Now the PDE (20) is transformed to a system of algebraic equations in unknowns. The boundary conditions are collocated separately to yield algebraic equations. More precisely,
where and are the Sinc points defined on and , respectively. Adding these equations to the 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 in the collocation steps of the homogenous boundary conditions. This factor emphasizes the boundary values and improve the error behavior at the boundaries.

The PolySinc collocation technique is based on the collocation of the spatial variables using Sinc points. This means that it is valid also for PDEs with spacedependent 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
(21)  
where is the spatial domain and is an event space and is a random variable. The function is a linear function of a uniformly distributed random variable and for all .
Now, we use the PC representation in (4) with to have
(22) 
where ’s are the univariate orthonormal Legendre polynomials defined on . Substitution of (22) in the SPDE (21) yields the residual
We then perform a Galerkin projection and use the orthogonality of Legendre polynomials, which yields the system of elliptic PDEs
(23)  
where . It holds that .
The computational results of this example are given in the following experiments.
Experiment 1.
and
In this experiment, we use PolySinc collocation from Section 4 to solve the system of PDEs in (23). In our computation, we use , i.e. of 2D grid of Sinc points defined as in (12) on the domain . As a result of the PolySinc solution, the coefficient functions are obtained. In Fig. 1, the expectation 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, , is needed. In our computations, we used , i.e. four orthonormal Legendre polynomials. The coefficients functions, , are given in Fig.3. In addition, we verify that this number is sufficient by showing that the coefficient functions tend to zero as increases. The results are given in Fig.4. In Fig.4, the dots represent the maximum of the coefficient functions on the spatial domain. We then use these maximum values in a least square
estimation to find the coefficients of the decaying rate function , where and are constants. In Fig.4, the solid line represents the best fitting function with and . This means that the coefficient functions follow an exponentially decay relation.
Experiment 3.
Error
To discuss the convergence of PolySinc 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 . We choose points of GaussLegendre nodes as values of and create corresponding PDEs. To solve each one of these equations we use Mathematica Package NDSolve. NDSolve uses a combination of highly accurate numeric schemes to solve initial and boundary value PDEs ^{1}^{1}1For 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 and using and PolySinc (with orthonormal Legendre) and the references from the PDEs are presented. Using the spatial norm error, calculating the error in both and delivers error of order and , 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 .
Experiment 4.
Comparison
In this experiment we compare the PolySinc solution with the classical finite difference (FD) solution. In 5pointstar FD method [24], we use an meshing with constant step size for the spatial variables and , which is the same number of Sinc points used in the PolySinc solution. The error between finite difference solution and the reference exact solution is given in Fig.7. Using the spatial norm error, calculating the error in both and delivering error of order . These calculations shows that for the same number of points, PolySinc delivers better approximation for the solution of the SPDE. In Fig.8 we run the calculations for different numbers of Sinc points and use the same number of points in the FD method. We then calculate the norm error. Fig.8 shows that the decaying rate of the error, in both mean and variance, is better in PolySinc than the FD method. Moreover, the PolySinc 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 in (2) and where,
is a set of independent random variables uniformly distributed in . For this SPDE we run four experiments.
Experiment 5.
and
In this experiment, we perform the Galerkin method along side the multivariate PC. For the PC parameters, we choose and . Due to (8), the number of multivariate Legendre polynomials is . As a result the threedimensional array is of dimension . For the PolySinc solution of the resulting system of PDEs, we use , i.e. 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 PolySinc 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 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 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 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 ^{2}^{2}2For 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 norm error, calculating the error in both and deliver error of order and , respectively.
Experiment 8.
Comparison
In this experiment we compare the PolySinc solution with the 5pointstar FD method. The reference solution is the Finite Element (FE) solution with cell meshing . In Fig.14 we run the calculations for different numbers of Sinc points and use the same number of points in FD. We then calculate the norm error. These calculations show that the decaying rate of the error, in both mean and variance, is better in PolySinc than the FD method. Moreover, the PolySinc 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 PolySinc 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 PolySinc 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.
References
 [1] R. J. Adler, The Geometry of Random Fields, John Wiley and Sons, Chichester, (1981).
 [2] T. G. Theting, Solving Wickstochastic boundary value problems using a finite element method, Stochastics and Stochastics Reports 70(3â4), 241270, (2000).
 [3] R. E. Caflisch, Monte Carlo and quasiMonte Carlo methods, Acta Numerica, 7, 149, (1998).
 [4] M. Kleiber, T. D. Hien, The Stochastic Finite Element Method, Basic Perturbation Technique and Computer Implementation, JohnWiley and Sons, Chichester, (1992).
 [5] H. Niederreiter, Random Number Generation and QuasiMonte Carlo Methods, Philadelphia, PA, SIAM, (1992).
 [6] I. Babuka, P. Chatzipantelidis, On Solving Elliptic Stochastic Partial Differential Equations, Comp. Meth. Appl. Mech. Engrg. 191, 40934122, (2002).
 [7] R. Ghanem, P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer, Berlin, (1991).
 [8] D. Xiu, G. E. Karniadakis, The WienerAskey Polynomial Chaos for Stochastic Differential Equations, SIAM J. Sci. Comput., 24, 619644, (2002).
 [9] R. Ghanem, Ingredients for a general purpose stochastic finite elements implementation, Comp. Meth. Appl. Mech. Engrg. 168, 1934, (1999).
 [10] D. Xiu, D. Lucor, C.H. Su, G. E. Karniadakis, Stochastic Modeling of FlowStructure Interactions Using Generalized Polynomial Chaos, ASME J. Fluid Engrg. 124, 5169, (2002).
 [11] R. Pulch, Stochastic Collocation and Stochastic Galerkin Methods for Linear Differential Algebraic Equations, J. Comput. Appl. Math., 262, 281291, (2014).
 [12] J. A. S. Witteveen, H. Bijl, Modeling Arbitrary Uncertainties Using GramSchmidt Polynomial Chaos, In Proceedings of the 44th AIAA Aerospace Sciences Meeting and Exhibit, Number AIAA20060896, Reno, NV, USA, 9â12, (2006).
 [13] M.S. Eldred, Recent Advances in NonIntrusive Polynomial Chaos and Stochastic Collocation Methods for Uncertainty Analysis and Design, In Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, CA, USA, 47 May 2009.
 [14] R. Pulch, Polynomial chaos for boundary value problems of dynamical systems, Appl. Numer. Math. 62(10), 14771490, (2012).
 [15] D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press: Princeton, NJ, USA, (2010).
 [16] F. Stenger, Handbook of Sinc Methods, CRC Press, (2011).
 [17] F. Stenger, M. Youssef, J. Niebsch, Improved Approximation via Use of Transformations, In: Multiscale Signal Analysis and Modeling, Eds. X. Shen and A.I. Zayed, NewYork: Springer, pp. 2549, (2013).
 [18] J. Stoer, R., Bulirsch, Introduction to Numerical Analysis, Springer, New York, 3rd ed., (2002).
 [19] M. Youssef, H. A. ElSharkawy, G. Baumann, Multivariate PolySinc Approximation, Error Estimation and Lebesgue Constant, Journal of Mathematics Research, Canadian Center of Sc. and Ed., 8(4), (2016). http://dx.doi.org/10.5539/jmr.v8n4p118.
 [20] M. Youssef, G. Baumann, Collocation Method to Solve Elliptic Equations, Bivariate PolySinc Approximation, Journal of Progressive Research in Mathematics (JPRM), ISSN: 23950218, 7(3), pp. 10791091 (2016).
 [21] M. Youssef, H. A. ElSharkawy, G. Baumann, Lebesgue Constant Using Sinc Points. Advances in Numerical Analysis, 2016, Article ID 6758283, 10 pages, (2016). http://dx.doi.org/10.1155/2016/6758283
 [22] M. Youssef, G. Baumann, On Bivariate PolySinc Collocation Applied to Patching Domain Decomposition, Applied Mathematical Sciences, 11(5), pp. 209226, (2017).
 [23] M. Youssef, PolySinc Approximation Methods, PhD thesis, Math. Dept. German University in Cairo, (2017).
 [24] Ch. Grossmann, H.G. Roos, M. Stynes, Numerical Treatment of Partial Differential Equations, Springer, (2007).
 [25] C. J. Gittelson, An Adaptive Stochastic Galerkin Method For Random Elliptic Operators, Math. of Comp. 82(283), pp. 15151541, (2013).