A boundary element and radial basis function approximation method for a second order elliptic partial differential equation with general variable coefficients

A numerical method based on boundary integral equation and radial basis function approximation is presented for solving boundary value problems governed by a second‐order elliptic partial differential equation with variable coefficients. The equation arises in the analysis of steady‐state anisotropic heat or mass diffusion in nonhomogeneous media with properties that vary according to general smoothly varying functions of space. The method requires only the boundary of the solution domain to be discretized into elements. To check the validity and accuracy of the numerical solution, some specific problems with known solutions are solved.


INTRODUCTION
The analysis of two-dimensional steady-state heat or mass diffusion in an anisotropic medium may be formulated in terms of boundary value problems governed by the second-order linear elliptic partial differential equation (see, for example, the works of Bera et al 1 and Costa 2 ) where u is the temperature or concentration which is a function of the Cartesian coordinates x 1 and x 2 and k ij are the diffusion or conduction coefficients satisfying the symmetric relation k ij = k ji and the positive definiteness condition given by the strict inequality k i i > 0 for all real numbers k such that 2 1 + 2 2 ≠ 0.
In nonhomogeneous media such as a functionally graded material, k ij may be taken to be varying smoothly from point to point in space. Researchers such as Al-Jawary and Wrobel, 3 Chakrabortya et al, 4 Fahmi, 5,6 Rangelov et al, 7 and Reutskiy 8 have proposed various numerical methods for the analyses of such nonhomogeneous materials.
Integral equation based methods like the boundary element method are of interest here as they offer certain advantages such as computational efficiency and accuracy in the treatment of the governing differential equations and the boundary conditions. Boundary element solutions for the boundary value problems governed by (1) are well established for homogeneous media, that is, for constant coefficients k ij (see, for example, the works of Clements 9 and Ooi et al 10 ). The development of boundary element methods for the more general case in which the coefficients k ij are continuously varying functions of x 1 and x 2 is, however, a more challenging task.
In many papers on boundary integral methods for solving (1) for nonhomogeneous media, the diffusion or conduction coefficients are taken to be of the form k i = i g(x 1 , x 2 ), where ij are constants, that is, the variation or grading of the coefficients is determined by just a single function g(x 1 , x 2 ). The perturbation boundary element approach proposed by Rangogni 11 may be used to solve approximately (1) for slightly varying coefficients k ij with g(x 1 , where is a constant parameter of an extremely small magnitude and f(x 1 , x 2 ) is a given smoothly varying function. For heat conduction in isotropic solids, Ang et al 12 and Clements 13 derived special fundamental solutions for (1) with k ij = ij X(x 1 )Y(x 2 ) (where ij is the Kronecker-delta and X(x 1 ) and Y(x 2 ) are given smoothly varying functions) and Kassab and Divo 14 introduced the idea of a generalized fundamental solution. With a suitable fundamental solution, a boundary integral equation needed for developing a boundary element method may then be derived for the partial differential equation in (1).
From a mathematical standpoint, suitable fundamental solutions in analytic closed forms are, however, inherently difficult, if not impossible, to derive for general variable coefficients k ij . If the fundamental solution for the special case where k ij are constants (homogeneous media) is used instead to obtain an integral equation for (1) with k ij given by smoothly varying functions, the resulting integral formulation contains not only a boundary integral but also a domain integral containing the unknown function u in the integrand. For k i = i g(x 1 , x 2 ) (with constant ij ), Ang et al 15 and Tanaka et al 16 applied the dual-reciprocity method proposed by Brebbia and Nardini 17 to approximate the resulting domain integral in terms of a boundary integral, in order to avoid the need to discretize the solution domain into elements.
In this paper, we present a boundary element procedure for a more general steady-state anisotropic diffusion equation by including a linear source term and taking the coefficients k ij to be any general smoothly varying functions of space. No restrictive form (such as k i = i g(x 1 , x 2 ) with constant ij , as assumed in, for example, the works of Ang et al 15,18 Dineva et al, 19 and Rangelov et al 7 ) is imposed on k 11 , k 12 , and k 22 here. The coefficients k 11 , k 12 , and k 22 may be individually given by any smoothly varying functions as long as they satisfy the positive definiteness condition (2) (for elliptic partial differential equations) in the solution domain.
Our solution approach here employs radial basis functions to approximate u(x 1 , x 2 ) and related functions in order to rewrite the governing partial differential equation as a constant coefficient linear elliptic partial differential equation which has a standard boundary integral equation. Unlike the dual-reciprocity boundary element approaches for nonhomogeneous media as in, for example, the works of Ang et al, 15 Fahmy, 20 and Tanaka et al, 16 the integral formulation here for the problem under consideration does not involve any domain integral. The radial basis function approximations and the discretization of the boundary integral equation together with the boundary conditions give rise to a system of linear algebraic equations for the approximate solution of the boundary value problem under consideration here. The numerical procedure does not require the entire solution domain to be discretized into elements. Only the boundary of the solution domain is discretized into elements. To check the validity and accuracy of the numerical procedure, specific problems with known solutions are solved numerically.

BOUNDARY VALUE PROBLEM
The boundary value problem of interest here is to solve the steady-state anisotropic diffusion equation with a linear source term, as given by in a two-dimensional bounded region R on the Ox 1 x 2 plane subject to where r and s are suitably prescribed functions, C 1 and C 2 are nonintersecting curves such that C = C 1 ∪ C 2 is the curve bounding the region R and n i is the x i component of the outward unit normal vector to the curve C.
The coefficients q i and k ij in (3) are given by continuous functions of x 1 and x 2 in the solution domain. As mentioned earlier, k ij may be given by any smoothly varying functions as long as the elliptic condition (2) is satisfied.
The special case of the boundary value problem for isotropic media (that is, the case of the Helmholtz equation where k 11 = k 22 and k 12 = k 21 = 0) is considered in the work of Al-Jawary and Wrobel. 3 In the aforementioned work, 3 the partial differential equation is formulated in terms of boundary-domain integro-differential equations and the resulting domain integral, which has a weak singularity, is approximated as boundary integral by using radial basis functions and the radial integration method in the work of Gao. 21 A different approach based on boundary integral equation and radial basis function approximation, which does not involve any domain integral in the solution formulation, is presented here in Section 3 below for solving (3) numerically subject to the boundary conditions in (4).

Reformulation of the partial differential equations
For solving the boundary value problem, we rewrite the elliptic partial differential equation in (3) as where k (0) i are constants which may be obtained by averaging k ij uniformly over the solution domain R and k (1) i are, in general, functions that vary smoothly with x 1 and x 2 in the solution domain. Motivated by the analysis in the works of Ang 18 and Dobroskok and Linkov, 22 we introduce the substitution where v is chosen to be related to u by and w satisfies the constant coefficient linear partial differential equation It is easy to show that (3) is satisfied by (6), (7), and (8).
Our solution approach is as follows. We employ a meshless technique based on radial basis functions to approximate (7) in terms of a system of linear algebraic equations. We discretize the boundary integral equation for (8), as in the work of Clements, 9 into linear algebraic equations. For the numerical solution of the boundary value problem in Section 2, we solve the resulting linear algebraic equations by taking into account the boundary conditions in (4).

Radial basis function approximation
For the meshless technique for approximating (7), we choose P well spaced out collocation points in R ∪ C, denoting the chosen points by ( (1) is the x i coordinate of the jth collocation point.
We make the approximation where a (r) i are constant coefficients and (r) (x 1 , x 2 ) is a radial basis function centered about ( (r) 1 , (r) 2 ). One may refer to the works of Dehghan and Mohammadi 23 and Fasshauer 24 for details on radial basis functions.
Substituting (9) into (7), we obtain If we make the approximations 2 ) for n = 1, 2, … , P, and invert the resulting linear equations for the constants b (r) and c (r) , we obtain where If we substitute (12) into (9) and collocate at (x 1 , where Inversion of (14) gives where If we substitute (16) into (10) and collocate the resulting equation at each of the chosen collocation points, we obtain where The linear algebraic equations in (18) with unknowns u (m) and v (m) (m = 1, 2, … , P) may be regarded as a radial basis function approximation of the partial differential Equation (7).
In the formulation above, the radial basis function (r) (x 1 , x 2 ) is required to be partially differentiable once with respect to x 1 or x 2 . We may take (r) (x 1 , x 2 ) to be the well-known multiquadric radial basis function given by (see the works of Ferreira 25 and Sarra 26 ) Alternatively, one may use the radial basis function proposed in the work of Zhang and Zhu, 27 that is,

Boundary integral approximation
The elliptic partial differential equation in (8) can be recast into the boundary integral equation (see the work of Clements 9 ) where lies in the interior of the solution domain R bounded by the curve C and ( 1 , 2 ) = 1∕2 if ( 1 , 2 ) lies on a smooth part of the curve C, and where Re denotes the real part of a complex number. Note that k (0) (6) into (22), we obtain where To approximate (24), we discretize the boundary C into M straight line elements denoted by C (1) , C (2) , We make the approximations and where p (m) = p( (m) 1 , (m) 2 ) and t (m) = t( (m) 1 , (m) 2 ) for m = 1, 2, … , M. Note that u (n) and v (n) are respectively the values of u (x 1 , x 2 ) and v(x 1 , x 2 ) at the nth collocation point as defined in Section 3.2.

Numerical procedure
For the numerical solution of the boundary value problem in Section 2, we approximate (7) and (8)

by (18) (with P = M+N)
and (28), respectively. The total number of equations in (7) and (8) is given by 2M + 2N. It is less than the number of unknowns involved. The number of unknowns is 4M + 2N. Thus, another 2M equations are required to complete the numerical formulation.
The boundary conditions in (4) can be expressed in terms of M linear algebraic equations given by or with the constant coefficients y (mq) defined by where n (m) i is the x i component of the unit vector that is normal to C (m) and that points out of the solution domain. Note that (30) is derived from the second line in (4) by taking k i = k (0) i + k (1) i and using (12) and the first equation in (25). Another M linear algebraic equations may be derived by using (12) and the second equation in (25). They are given by where For the numerical solution of the boundary value problem in Section 2, we solve the system of 4M + 2N linear algebraic equations given by (18) (with P = M + N), (28), (29) or (30), and (32). If one wishes to solve a smaller system of linear algebraic equations, one may substitute (29) or (30) together with (32) directly into (28). This will reduce the number of unknowns and the number of equations in the formulation to 2M + 2N. The computer coding of the smaller system of equations is, however, more involved. Once the unknowns are determined, the required solution u(x 1 , x 2 ) can be computed approximately in an explicit manner at any general point (x 1 , x 2 ) in the solution domain by using (12).

SPECIFIC PROBLEMS
The numerical procedure described above is applied to solve the following specific problems. The radial basis function given in (21) is used in the numerical computation. Problem 1. Take the coefficients q i and k ij to be given by The problem is to solve the partial differential equation given by (3) together with (34) in the solution domain 0 < x 1 < 1, 0 < x 2 < 1, subject to the boundary conditions It may be verified by direct substitution that the analytical solution of the problem is given by For the numerical solution of the problem, each side of the square solution domain is discretized into M 0 equal length elements and the interior collocation points are chosen to be well spaced out points given by (i∕(N 0 +1), j∕(N 0 + 1)) for i, j = 1, 2, … , N 0 (so that M = 4M 0 and N = N 2 0 ). We take k (0) i to be the average value of k ij over all the interior collocation points, that is, In Table 1 For the calculation of the flux, the approximate formula for u in (12) may be partially differentiated with respect to x i to compute the first-order partial derivatives of u with respect to x 1 and x 2 denoted by u ,1 and u ,2 , respectively. The numerical values of u from the calculation using (M 0 , N 0 ) = (40, 19) are used in (12) to compute the approximate values of u ,1 and u ,2 at selected interior points. In Table 2, the approximate values obtained are compared with the values calculated from the analytical solution in (36). As may be expected, for the same set of boundary elements and interior collocation points, the accuracy of the numerical values of u ,1 and u ,2 is less than that of u since u ,1 and u ,2 are secondary quantities obtained by postprocessing the numerical values of u at the collocation points. Specifically, the average absolute errors of     Table 2 are about 15 times larger than that of the numerical values of u in the third column of Table 1.
The numerical procedure is executed on a desktop personal computer with an entry level Intel Core i3-4150 processor and solve (3) in the solution domain 0 < x 1 < 1, 0 < x 2 < 1 subject to the boundary conditions As in Problem 1, the boundary of the square solution domain is discretized into 4M 0 equal length elements and N 2 0 interior collocation points are selected for the numerical calculation. In addition, the values of k (0) i are computed using the formula in (37).
It is easy to check that the analytical solution of the problem here is given by Once the numerical values of u at all the collocation points, that is, the values of u (m) (m = 1, 2, … , 4M 0 + N 2 0 ) are known, the formula for u in (12) can be used to compute approximately the solution at any point of interest in the solution domain. In Table 3 in the quarter circular domain x 2 1 + x 2 2 < 1, x 1 > 0, x 2 > 0, subject to the boundary conditions Note that c in (41) is a real constant.
For the numerical solution of the boundary value problem, each of the three parts of the boundary, namely, where x 1 = 0, x 2 = 0 and x 2 1 + x 2 2 = 1, is discretized into M 0 boundary elements (so that M = 3M 0 ) and N 2 0 collocation points in the interior solution domain are taken to be given by for m, n = 1, 2, … , N 0 .
For (M 0 , N 0 ) = (40, 6), the numerically obtained values of the u on x 2 1 + x 2 2 = 1, x 1 > 0, x 2 > 0 (where k 11 x 1 u∕ x 1 + k 22 x 2 u∕ x 2 is specified) are plotted against = arctan(x 2 ∕x 1 ) for selected values of c in Figure 1. The plots in the figure  FIGURE 1 Plots of u on x 2 are very close to those given in the work of Ang et al, 15 where the same problem in the context of a functionally graded material undergoing an antiplane deformation was solved numerically using a different boundary element procedure.

SUMMARY AND FINAL REMARKS
A numerical method based on boundary integral equation and radial basis function approximation is derived and successfully implemented in the computer for solving a two-dimensional steady-state anisotropic heat or mass diffusion equation with a linear source term and general variable coefficients. The numerical procedure requires only the boundary of the solution domain to be discretized into elements. However, the formulation does not involve only unknowns at boundary collocation points but also at well-distributed collocation points in the interior of the solution domain. The validity and accuracy of the boundary element procedure is verified by applying it to solve some specific problems with known solutions. The numerical solutions obtained show convergence and agree well with the known solutions.
In this paper, the boundary integral equation in the formulation is discretized using only constant elements. With constant elements, the accuracy of the numerical solution is expected to be O(h) (where h is the length of a typical element). This is reflected in the reduction in the errors of the numerical solutions of the specific test problems when the computation is refined by increasing the number of boundary elements and interior collocation points. In spite of the constant elements, we have managed to obtain reasonably accurate numerical solutions for the test problems even with a relatively low number of boundary elements and interior collocation points. Higher order elements such as the discontinuous linear elements may be employed in the boundary integral equation if a more accurate numerical solution is desired such as in the postprocessing computation of the flux. The implementation of higher order elements is, however, algebraically more tedious.
The proposed solution approach based on boundary integral and radial basis function approximations provides an interesting and viable alternative to the existing numerical methods for analyzing heat and mass transfer in anisotropic media with spatially varying material properties. It may be explored further and extended to solve more complicated problems involving anisotropic media such as those considered in the works of Aksoy andŞenocak, 29 Baron, 30 Lin et al, 31 Marin and Lesnic, 32 and Fahmy. 33,34 How to cite this article: Ang W-T. A boundary element and radial basis function approximation method for a second order elliptic partial differential equation with general variable coefficients. Engineering Reports. 2019;1:e12057. https://doi.org/10.1002/eng2.12057

2
) are the endpoints of the element C (m) , the integrals over C (m) in (28) can be rewritten as (G (m) t + H (mn) )dt A (mn) t 2 + B (mn) t + E (mn) , Note that n (m) i is the x i component of the unit vector that is normal to C (m) and that points out of the solution domain. For m ≠ n, the expression 4A (mn) E (m) − (B (mn) ) 2 is strictly greater than zero, and hence, the integrals in (A1) and (A2) are proper. The proper integrals may be evaluated either approximately by using a numerical integration formula or by using the analytical formulae ∫ ln(at 2 + bt + c)dt = t(ln(a) − 2) + , ∫ (gt + h)dt at 2 + bt + c = g 2a ln(at 2 + bt + c) , which are valid for 4ac − b 2 > 0. For n = m, the integrals on the right-hand side of (A1) and (A2) are improper to be interpreted in the Cauchy principal sense (with an integrable singularity at t = 0) since B (mm) = 0, E (mm) = 0, and H (mm) = 0. The resulting Cauchy principal integrals are, however, straightforward to evaluate analytically since their integrands are in considerably simple forms.