RBF stencils for poisson equation

Đặng Thị Oanh<br /> <br /> Tạp chí KHOA HỌC & CÔNG NGHỆ<br /> <br /> 78(02): 63 - 66<br /> <br /> RBF STENCILS FOR POISSON EQUATION<br /> Dang Thi Oanh*<br /> Faculty of Information Technology - TNU<br /> <br /> ABSTRACT<br /> In the paper we present a method for finding the weight vector called stencil with the help of RBF<br /> interpolation. This stencil is the foundation for constructing meshless finite difference scheme for<br /> boundary value problems. The results of numerical experiments show that the numerical solution<br /> obtained by RBF-FD with the stencils generated by Gauss RBF interpolation is much more<br /> accurate then the solution obtained by FEM.<br /> Keyword: Radial Basis Function (RBF), meshless method, shape parameter, stencil<br /> <br /> INTRODUCTION*<br /> Because of the difficulties to create, maintain<br /> and update complex meshes needed for the<br /> standard finite difference, finite element or<br /> finite volume discretisations of the partial<br /> differential equations, meshless methods<br /> attract growing attention. In particular, strong<br /> form methods such as collocation or<br /> generalised finite differences are attractive<br /> because they avoid costly numerical<br /> integration of the non-polynomial shape<br /> functions on non-standard domains often<br /> encountered in those meshless methods that<br /> are based on the weak formulation of PDE.<br /> Thanks to their excellent local approximation<br /> power, radial basis functions are an ideal tool<br /> to produce numerical differentiation stencils<br /> for the Laplacian and other partial differential<br /> operators on irregular centres, without any<br /> need for a mesh. This leads to exceptionally<br /> promising RBF based generalised finite<br /> difference method.<br /> The essense of the method is to find a weight<br /> vector called stencil corresponding a row of<br /> stiffness matrix in FEM. In this paper we give<br /> formulas for finding the stencil based on<br /> RBF, and then we use them for discretising<br /> the Poisson equation on a nonuniform set of<br /> centres. In [1], [2] for the above purpose we<br /> used RBF interpolation with additional<br /> polynomial term and performed many<br /> numerical examples in complicated geometry<br /> domains with several RBF. In this paper we<br /> *<br /> <br /> Email: dtoanhtn@gmail.com<br /> <br /> concentrate on RBF interpolation without<br /> additional polynomial term and test the<br /> Poisson equation on the disk domain with<br /> RBF-Gaussian. The results of numerical<br /> experiments show that the numerical solution<br /> obtained by RBF-FD with the stencil<br /> generated by Gauss RBF interpolation is<br /> much more accurate then the solution<br /> obtained by FEM.<br /> The paper is organised as follows. In Section<br /> 2 we describe stencils from RBF<br /> interpolation. Section 3 is devoted RBF-FD<br /> discretisation of Poisson equation and finally<br /> Section 4 we provide the results of the<br /> numerical tests with the near-optimal shape<br /> parameter on the disk domain.<br /> STENCILS FROM RBF INTERPOLATION<br /> Definition 1 (stencil) Let D be a linear<br /> n<br /> differential operator, and X  xi i 1 a fixed<br /> irregular set of centres in Rd. A linear<br /> numerical differentiation formula for the<br /> operator D,<br /> n<br /> <br /> Du( x)   wi ( x)u ( xi ) ,<br /> <br /> (1)<br /> <br /> i 1<br /> <br /> is determined by the weights wi = wi (x) .<br /> T<br /> The vector w  w1 ,, wn  is called stencil.<br /> Definition 2 (positive definite function) A<br /> d<br /> continuous function  : R  R is called<br /> positive semi-definite if, for all n  N , all<br /> sets of pairwise<br /> distinct centers<br /> X  x1 ,, xn   R d , and all c = (c1, c2,<br /> …, cn)  Rn, the quadratic form<br /> 63<br /> <br /> Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br /> <br /> http://www.lrc-tnu.edu.vn<br /> <br /> Đặng Thị Oanh<br /> <br /> Tạp chí KHOA HỌC & CÔNG NGHỆ<br /> <br />   c c x<br /> n<br /> <br /> n<br /> <br /> j 1<br /> <br /> k 1<br /> <br /> j k<br /> <br /> j<br /> <br />  xk <br /> <br /> (2)<br /> <br /> is nonnegative. Function  is called positive<br /> definite if the quadratic form is positive for all<br /> c  R n \ 0.<br /> Definition 3 (conditionally positive definite<br /> function)<br /> A<br /> continuous<br /> function<br />  : R d  R is said to be conditionally<br /> positive semi-definite of order m (i.e. to have<br /> conditional positive definiteness of order m)<br /> if, for all n  N , all sets of pairwise distinct<br /> centers X  x1 ,, xn   R d , and all<br /> c  R n satisfying<br /> n<br /> <br /> c<br /> j 1<br /> <br /> j<br /> <br /> p( x j )  0<br /> <br /> n<br /> <br />  c c ( x<br /> <br /> j ,k 1<br /> <br /> j k<br /> <br /> j<br /> <br />  xk )<br /> <br /> is nonnegative,  is said to be conditionally<br /> positive definite of order m if the quadratic<br /> form is positive, unless c is zero.<br /> Definition 4 (Radial Basis Function (RBF))<br /> d<br /> In the Euclidean space R setting, an RBF is<br /> a function of the form<br /> n<br /> <br /> <br /> <br /> x   c j x  x j<br /> j 1<br /> <br /> 2<br /> <br /> ,<br /> <br /> (3)<br /> <br /> where x1 ,, xn are some scattered points,<br /> and x - x 2 denotes the Euclidean<br /> distance between x and x j ; c1,, cn are<br /> some constants, and<br /> function<br /> <br /> <br /> <br /> is a univariate<br /> <br />  : 0,   R.<br /> <br /> (4)<br /> <br /> Proposition 1 Assume D be a linear<br /> differential operator. Let  : R  R is a<br /> possitive definite function, set of pairwise<br /> distinct centers X  x1 ,, xn   R d and a<br /> function u : R d  R , the RBF interpolation<br /> s is sought in the form<br /> <br /> s( x)   a j x  x j , ( x) :   x<br /> n<br /> <br /> j 1<br /> <br /> Then, the stencil w for the numerical<br /> differentiation at x is given by<br /> 1<br /> (6)<br /> w   | X  D( x  ) | X ,<br /> Where<br /> <br /> <br /> <br /> <br /> <br /> (5)<br /> <br /> <br /> <br />  | X  ( xi  x j ) i , j 1 , D( x  ) | X <br /> n, n<br /> <br /> D( x  x1 ),, D( x  nxn )T<br /> <br /> <br /> <br /> <br /> <br /> <br /> <br /> Proof<br /> Since s( x) <br /> a j  x  x j is<br /> j<br /> <br /> 1<br /> interpolation function of function u, we have<br /> <br /> s( xi )  u( xi ), i  1,, n,<br /> that is<br /> n<br /> <br />  a ( x<br /> j 1<br /> <br /> for all real-valued polynomials p( x) of<br /> degree less than m, the quadratic form<br /> <br /> 78(02): 63 - 66<br /> <br /> j<br /> <br /> i<br /> <br />  x j )  u ( xi ), i  1,, n ,<br /> <br /> ( 7)<br /> which can be written in matrix form as the<br /> linear equation<br />  |X a  u |X ,<br /> with a symmetric positive definite matrix,<br /> where<br /> <br /> <br /> <br /> <br /> <br />  | X  ( xi  x j ) i , j 1 ,<br /> T<br /> u | X  u( x1 ),, u( xn ) ,<br /> T<br /> a  a1 ,, an  .<br /> n<br /> <br /> Since the function u is sufficiently smooth<br /> and the set of points x1 ,, xn  R d is<br /> sufficiently dense in a neighbourhood of x<br /> RBF interpolant s(x) provides a good<br /> approximation of u(x) [3]. Moreover, the<br /> derivatives of s are good approximations of<br /> the derivatives of u if  is sufficiently<br /> smooth [3]. Therefore, an approximation of<br /> Du(x), where D is a linear differential<br /> operator annihilating constants, may be<br /> considered in the form<br /> n<br /> <br /> n<br /> <br /> j 1<br /> <br /> i 1<br /> <br /> Du( x)  Ds( x)   a j D( x  x j )   wi u( xi )<br /> (8)<br /> where the weights wi (depending on x) exist<br /> because the coefficients a j of the<br /> interpolation function s defined by (5)-(7)<br /> depend linearly on the data u( xi ), i  1,, n<br /> . These weights can be found by solving the<br /> symmetric positive definite linear system<br />  | X w  D( x  ) | X , (9)<br /> that is<br /> <br /> 64<br /> <br /> Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br /> <br /> http://www.lrc-tnu.edu.vn<br /> <br /> Đặng Thị Oanh<br /> <br /> Tạp chí KHOA HỌC & CÔNG NGHỆ<br /> <br /> w   | X  D( x  ) | X<br /> <br /> ~<br />  DD( x  yi ), i  1, m.<br /> <br /> 1<br /> <br /> Proposition 2. Assume D be a linear<br /> differential operator. Let  : R  R is a<br /> positive definite or conditionally positive<br /> definite function order l  1 , set of pairwise<br /> distinct centers X  x1 ,, xn   R d and a<br /> function u : R d  R , the RBF interpolation<br /> s is sought in the form<br /> <br /> s( x)   a j x  x j   ~p ( x), ( x) :   x  ,<br /> n<br /> <br /> j 1<br /> <br /> where ~<br /> p is a polynomial of degree l . Then,<br /> the stencil w can be found by solving the<br /> symmetric positive definite linear system<br /> <br />  | X<br /> T<br /> 1<br /> where<br /> <br /> <br /> <br /> RBF-FD DISCRETISATION OF POISSON<br /> EQUATION<br /> In the finite difference method stencils are<br /> used for the discretisation of partial<br /> differential equations. Consider the Dirichlet<br /> problem for the Poisson equation in a<br /> bounded domain   R d : given a function f<br /> defined on  , and a function g defined on<br />  find u such that<br /> (10)<br /> u  f on <br /> <br /> u |  g.<br /> <br /> <br /> <br /> that for each    int a set<br /> <br /> n, n<br /> <br /> D( x  x1 ),, D( x  xn ),T , 1 : [11]T<br /> Proposition 3 Assume D be a differential<br /> operator of order k . Let  : R  R is a<br /> positive definite and 2k times continuously<br /> differentiable<br /> function.<br /> Let<br /> X  x1 ,, xn   R d<br /> and<br /> Y  y1 ,, ym   R d are two point sets in<br /> R d (which may coincide or overlap), and<br /> give a function u : R d  R . Consider the<br /> RBF interpolation s is sought in the form<br /> j 1<br /> <br /> ( x) :   x  ,<br /> ~<br /> <br /> the<br /> by<br /> for<br /> be<br /> <br /> n<br /> <br /> m<br /> ~<br /> w<br /> <br /> (<br /> x<br /> <br /> x<br /> )<br /> <br />  j i j  j D( xi  x j )  D( x  xi )<br /> j 1<br /> <br /> j 1<br /> <br /> i  1, n,<br /> n<br /> <br /> m<br /> <br /> j 1<br /> <br /> j 1<br /> <br /> ~<br />  w j D( yi  x j )   j DD( yi  y j )<br /> <br />   <br /> <br /> is<br /> <br /> chosen such that     and<br /> <br /> <br /> <br />  <br /> <br /> (12)<br /> <br />   int<br /> <br /> For each    int , choose a linear numerical<br /> differentiation formula for Laplace operator<br /> ,<br /> (13)<br /> u ( )   w , u ( ),<br />  <br /> w ,   , and replace (10)with stencil<br /> <br /> <br /> <br /> <br /> <br /> <br /> <br /> (11) by the system of linear equations<br /> <br />  w  uˆ( )  f ( ),<br /> <br /> n<br /> m<br /> ~<br /> s H ( x)   a j x  x j    b j D( x  y j ),<br /> <br /> where the operator D is defined by<br /> ~<br /> formula D  RD , where R is defined<br /> Ru ( x)  u( x) . Then, the stencil w, <br /> the numerical differentiation at x can<br /> found from the linear system<br /> <br /> (11)<br /> <br /> This problem can be discretised with the help<br /> of differentiation formulae (1) as follows. Let<br />    be a finite set discretisation centres,<br />  :    and  int   \  . Assume<br /> <br /> 1 w  D( x  ) | X <br />  <br /> ,<br /> 0  v   0<br /> <br /> <br />  | X  ( xi  x j ) i , j 1 , D( x  ) | X <br /> <br /> j 1<br /> <br /> 78(02): 63 - 66<br /> <br />  <br /> <br /> ,<br /> <br />    int (14)<br /> <br /> (15)<br /> uˆ ( )  g ( ),   <br /> If (14)-(15) is nonsingular, then its solution<br /> uˆ :   R can be compared with the vector<br /> u |  u( ) of the discretised exact<br /> solution of (10)-(11).<br /> A standard finite difference method is<br /> obtained from the above if we take   R d<br /> to be a square domain,  a uniformly spaced<br /> grid, and (13) the classical 5-point<br /> differentiation formula for the Laplacian.<br /> NUMERICAL EXAMPLE<br /> In the numerical results of this paper we<br /> restrict our attention to the Gaussian RBF<br /> 65<br /> <br /> Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br /> <br /> http://www.lrc-tnu.edu.vn<br /> <br /> Đặng Thị Oanh<br /> <br /> Tạp chí KHOA HỌC & CÔNG NGHỆ<br /> <br />  (r )  e ( r ) , which is positive definite for<br /> <br /> 78(02): 63 - 66<br /> <br /> 2<br /> <br /> any value of the shape parameter<br />  . For<br /> definiteness we take   1.24 , which is nearoptimal parameter (see [2]). RBF-FD<br /> meshless method for solving Possion equation<br /> as following: For each    int , stencil<br /> w ,   is found by solving the symmetric<br /> <br /> positive definite linear system (6), and replace<br /> (10)-(11) and we have (14)-(15).<br /> Test Problem Consider Dirichlet problem<br /> (10)-(11) with boundary conditions defined<br /> by the restriction of the test function on the<br /> boundary, where the domain  is disk with<br /> centre (0, 0), and radius r  1 , the right hand<br /> side<br /> is<br /> defined<br /> by<br /> 2<br /> f ( x, y)  2 sin(x) sin(y) . The exact<br /> solution is u( x, y)  sin(x) sin(y) .<br /> Numerical experiments To assess the quality<br /> of a discrete solution uˆ of the Dirichlet<br /> problem defined on a set of discretisation<br /> centres  , we consider the root mean square<br /> (rms) error against the values of the exact<br /> solution u on  int ,<br /> <br /> <br /> <br /> <br /> <br />  1<br /> rms : <br />  #  int<br /> <br />  uˆ( )  u( )<br /> <br /> 2<br /> <br />  1/ 2<br /> . (16)<br /> <br /> <br /> <br /> Apart from the RBF-FD solutions, this<br /> formula applies to the standard linear finite<br /> element method with midpoint quadrature<br /> rule on the corresponding triangulation. We<br /> will use rms of the finite element method as<br /> reference.<br /> In Figure 1, FEM corresponds to finite<br /> element method, Gauss-FD corresponds to<br /> RBF-FD method (formulas (14)-(15)) with<br /> <br /> stencil w defined by (6) in Proposition 1 for<br /> RBF Gaussian.<br /> The results of experiments show that the error<br /> in the case of using RBF-FD method based<br /> on the stencil obtained by RBF is much better<br /> in comparison with FEM.<br /> <br /> Figure 1. RBF solution with Gaussian on the<br /> quasi-uniform FEM centres using FEM stencil<br /> support  <br /> REFERENCES<br /> [1]. O. Davydov and D. T. Oanh. Adaptive<br /> meshless centres and RBF stencil for Poission<br /> equation. J. Comput. Phys., (230) 287-304, 2011.<br /> [2]. O. Davydov and D. T. Oanh. On optimal<br /> shape parameter for Gaussian RBF-FD<br /> approximation for Poission equation. Technical<br /> Report 2011- 01. University of Strathclyde,<br /> Department of Mathematics and Statistics. Jan.<br /> 2011,<br /> http://www.mathstat.strath.ac.uk/research/reports/<br /> 2011.<br /> [3]. H. Wendland. Scattered Data Approximation.<br /> Cambridge University Press, 2005.<br /> <br /> TÓM TẮT<br /> VECTƠ TRỌNG SỐ CHO HÀM POISSON<br /> Đặng Thị Oanh*<br /> Khoa Công nghệ thông tin - ĐH Thái Nguyên<br /> <br /> Trong bài báo này chúng tôi giới thiệu phƣơng pháp tìm các véc tơ trọng số stencil với sự trợ giúp<br /> của các hàm nội suy RBF. Các stencil này là cơ sở xây dựng lƣợc đồ không lƣới RBF-FD. Thử<br /> nghiệm cho thấy trên cùng một tập tâm, sai số của phƣơng pháp không lƣới RBF-FD bởi các<br /> stencil RBF với tham số hình dạng gần tối ƣu tốt hơn đáng kể sai số của phƣơng pháp FEM.<br /> Từ khóa: Hàm cơ sở bán kính (RBF); phương pháp không lưới<br /> *<br /> <br /> Email: dtoanhtn@gmail.com<br /> <br /> 66<br /> <br /> Số hóa bởi Trung tâm Học liệu – Đại học Thái Nguyên<br /> <br /> http://www.lrc-tnu.edu.vn<br /> <br />