# Research on Error Estimations of the Interpolating Boundary Element Free-Method for Two-Dimensional Potential Problems

Authors: Jufeng Wang; Fengxin Sun; Ying Xu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The interpolating boundary element-free method (IBEFM) is a direct solution method of the meshless boundary integral equation method, which has high efficiency and accuracy. The IBEFM is developed based on the interpolating moving least-squares (IMLS) method and the boundary integral equation method. Since the shape function of the IMLS method satisfies the interpolation characteristics, the IBEFM can directly and accurately impose the essential boundary conditions, which overcomes the shortcomings of the original boundary element-free method in enforcing the essential boundary approximately. This paper will study the error estimations of the IBEFM for two-dimensional potential problems and the relationship between the errors and the influence radius and the condition number of the coefficient matrix. Two numerical examples are presented to verify the correctness of the theoretical results in this paper.


---

## Body

## 1. Introduction

The meshless method is an important numerical method after the finite element method, which has been widely used in scientific engineering calculation [1–5]. Based on different methods of constructing the shape function and problem discretization, scholars have proposed a series of meshless methods, such as the smoothed particle hydrodynamics (SPH) method [6, 7], element-free Galerkin (EFG) method [8], the singular boundary method (SBM) [9], radial basis functions [10, 11], boundary knot method (BKM) [12], meshless local Petrov–Galerkin (MLPG) method [13], complex variable element-free Galerkin (ICVFG) method [14, 15], point interpolation method (PIM) [16], method of fundamental solutions [17], and dimension splitting meshless method [18].The traditional least-squares method (LSM) is a commonly used method in numerical computation, which obtains the optimal solution by seeking the minimum deviation of all data [19, 20]. When studying surface fitting, based on the LSM, Lancaster and Salkauskas proposed the moving least-squares (MLS) method [21], which is one of the most important methods to construct the trial function in the meshless method. The meshless method based on the MLS method has high numerical accuracy [22, 23]. However, the shape function of the MLS method does not satisfy the property of the Kronecker δ function [24], which increases the difficulty of dealing with the boundary conditions, and also raises a lot of extra computational burden. The meshless method of the MLS method class cannot apply essential boundary conditions directly and simply as the finite element method [25–27]. To solve this problem, by using the singular weight function, Lancaster and Salkauskas proposed the interpolation moving least square (IMLS) method [21], in which the shape function meets the interpolation characteristics at the nodes. The meshless method based on the IMLS method can directly and simply impose the essential boundary conditions [28], which improve calculation accuracy and efficiency. Based on the IMLS method, different meshless methods are proposed to solve various problems [29–33].Based on the MLS method and boundary integral equation method, the boundary element-free method (BEFM) was developed by Cheng and Peng [34]. The BEFM is a meshless boundary integral equation method. Compared with the boundary point method, the boundary integral equation adopted in BEFM is a common and more regular form, and there is need neither to set Gauss points in the neighbourhood of the source point nor to choose the so-called “evaluation point,” so the solution process of the BEFM is more straightforward. The BEFM uses the real solutions of the node variables as the unknown quantities, while the boundary point method and the local boundary integral equation method both use the approximate solutions of the node variables as the unknown, so the BEFM is a direct numerical mesh-free method of the boundary integral equation and has greater computational precision [35–38].In the BEFM, the essential boundary conditions are directly imposed. However, because the shape function of the MLS method does not satisfy the interpolation characteristics at the nodes, the essential boundary conditions are only approximately satisfied, which may lead to reduction in calculation accuracy. In order to improve this deficiency, by taking into account the interpolation characteristics of the IMLS method, Ren et al. also proposed the interpolating boundary element-free method (IBEFM) [39], which has higher calculation accuracy than the BEFM [40, 41]. By using the improved IMLS method with a nonsingular weight function, Wang et al. also developed an enhanced IBEFM for the potential problems [42]. As a direct numerical mesh-free method, the IBEFM has both the advantages of high accuracy and the ability to directly and accurately apply boundary conditions.The BEFM has been widely used in various fields [37, 39, 43]. However, there is still little research on the error estimation of the BEFM. Based on the error estimation of the IMLS method and related matrix theory, this paper will study the error estimation of the IBEFM for the two-dimensional potential problems. The relationship between the errors and the influence radius and the condition number of the coefficient matrix is studied, and some numerical examples are presented to verify the correctness of the theoretical results.


## 2. The IBEFM for the Potential Problems

Consider the following two-dimensional Poisson equation,(1)∇2ux+bx=0,x=x1,x2∈Ω,with the Dirichlet boundary condition,(2)ux=u¯x,x∈Γu,and the Neumann boundary condition,(3)qx=∂ux∂n=q¯x,x∈Γq,where Ω is the problem field, ux represents the potential of the field function, bx is the given source function, Γ is the boundary of Ω satisfying Γ=Γu∪Γq, u¯ denotes the known potential function on the Dirichlet boundary Γu, q¯ is the known potential gradient on the Neumann boundary Γq, and n is the cosine of the outer normal direction on the boundary Γ.The integral equation can be obtained using the weighted residual method:(4)∫Ω∇2ux−bxu∗ξ,xdΩ=∫Γqqx−q¯xu∗ξ,xdΓ−∫Γux−u¯xq∗ξ,xdΓ,where x and ξ represent, respectively, the field point and source point and u∗ is the basic solution of the Poisson equation (1) satisfying(5)∇2u∗=−δξ,x,where δξ,x is the Dirac δ function. Then, it follows that(6)u∗ξ,x=12πln1rξ,x,(7)q∗ξ,x=∂u∗ξ,x∂nx,(8)rξ,x=ξ1−x12+ξ2−x22.Performing the partial integration on (4) and using the properties of the basic solution of the Poisson equation, when the source point ξ is within Ω, the following boundary integral equation can be obtained:(9)uξ=∫Γqxu∗ξ,xdΓ−∫Γuxq∗ξ,xdΓ−∫Ωbxu∗ξ,xdΩ,when the source point ξ is on the boundary Γ, and the boundary integral equation is(10)Cξuξ=∫Γqxu∗ξ,xdΓ−∫Γuxq∗ξ,xdΓ−∫Ωbxu∗ξ,xdΩ,where Cξ is a free term determined by the geometry of the point ξ, and(11)Cξ=12,whenξis a smooth boundary point,θ2π,whenξis a boundary corner point,where θ=θ2−θ1 and according to the boundary integration direction, θ1 and θ2 are the angles between the boundary tangent and the coordinate axis away from the corner point and near the corner point, respectively.For the two-dimensional potential problem, the points on the boundary can be seen as points in one-dimensional space.∀x∈Γ can be expressed by curvilinear coordinates s, that is x=xs. In this article, the curvilinear coordinates are taken as the arc length of the curve. Let X=x1,x2,…,xN denote the coordinate set of all nodes on the boundary, where N represents the total number of nodes. Suppose the arc length coordinate of xI to be sI, and this node is also called a node sI. The union of the influence domains of all nodes must cover the entire boundary, and for xI∈Γ, it is assumed that the influence domain RxI is(12)RxI:=xs∈Γ:s−sI≤ρI,where ρI is the influence radius of node xI.Supposes1,s2,…,snI are the nodes with the influence domain covering the xs. Then, from the IMLS method, the approximation functions can be obtained as follows:(13)ux≈uhxs=∑I=1nIΦIsuI,(14)qx≈qhxs=∑I=1nIΦIsqI,where uI=uxI and qI=qxI.Substituting (13) and (14) into (10), it follows that(15)CxIuxI=∫Γu∗xI,x∑J=1nIΦJsqJdΓ−∫Γq∗xI,x∑J=1nIΦJsuJdΓ−∫Ωbxu∗xI,xdΩ.To calculate the definite integral, the boundaryΓ is divided into some integral subdomains Γn, n=1,2,…,Ne. Ne is the total number of boundary integral subdomains, and only one point of any adjacent subdomains is connected. These integral subdomains are independent of the nodes distribution on the boundary, and they are only used to calculate the numerical integral.To calculate∫Ωbxu∗xI,xdΩ, the field Ω is divided into Nc integral subfields. Ac is the area of the subfield c, and K is the number of Gauss integral points in the integral subfield c. wk is the integral weight coefficient corresponding to the k-th Gauss integral point, and bu∗k is the value of bu∗. Therefore, equation (15) can be written as(16)CxIuxI=∑n=1Ne∫Γnu∗xI,x∑J=1nIΦJsqJdΓ−∑n=1Ne∫Γnq∗xI,x∑J=1nIΦJsuJdΓ−∑c=1Nc∑k=1Kwkbu∗kAc.By substituting each node on the boundary into (16) and combining the equations, the following equation can be obtained:(17)B+HU=GQ,where H=HIJN×N, G=GIJN×N, B=BIN×1, and(18)HIJ=∑n=1Ne∫Γnq∗xI,xΦJsdΓ+δIJCxI,GIJ=∑n=1Ne∫Γnu∗xI,xΦJsdΓ,BI=∑c=1Nc∑k=1Kwkbxu∗xI,xkAc,U=u1,u2,…,uNT,Q=q1,q2,…,qNT.In the calculation of the boundary integral equation (16), in most cases, the boundary source point xI is not on the integral subdomain Γn, and the integration has no singularity. The Gauss integration method can be used directly. If the source point coincides with the integration point, the basic solution will lead to the integration singularity. At this time, the Gauss integration point can be rearranged to eliminate the singularity.Because the shape function of the IMLS method satisfies the property of the Kroneckerδ function, the boundary condition can be easily and directly applied. Substituting the potential and potential gradient at the nodes on the boundary into (17), the potential and potential gradient at all the nodes on the boundary can be solved. Then, from (9), the potential at ∀ξ∈Ω can be obtained as(19)uξ=∑n=1Ne∫Γnu∗ξ,x∑I=1nIΦIsqIdΓ−∑n=1Ne∫Γnq∗ξ,x∑I=1nIΦIsuIdΓ−∑c=1Nc∑k=1Kwkbu∗kAc.Calculating differential for (19), the potential gradient of the interior point can be obtained:(20)u,iξ=∑n=1Ne∫Γnu,i∗ξ,x∑I=1nIΦIsqIdΓ−∑c=1Nc∫Γnq,i∗ξ,x∑I=1nIΦIsqIdΓ−∑c=1Nc∑k=1Kwkbu,i∗kAc,i=1,2.The abovementioned equation is the IBEFM for two-dimensional potential problems.


## 3. Error Estimations of the IBEFM for Two-Dimensional Potential Problems

For simplicity, letE1=x1,x2,…,xN1 be the set of all nodes on the Dirichlet boundary, and E2=xN1+1,xN1+2,…,xN be the set of all nodes on the Neumann boundary. For ∀xI∈X, assume that there is a constant cI such that ρ≤cIρI, and the basis functions used in the IMLS method are complete polynomials of order m.Since there is no intersection betweenΓu and Γq, then in (17), when the i-th element in U is unknown, the i-th element in Q must be known; otherwise, when the i-th element in U is known, the i-th element in Q must be unknown. By shifting the term, all the unknown quantities are moved to the left end, and the known quantities are moved to the right end. Suppose the unknown quantities to be solved at the right end be X, and its coefficient matrix be A, and at the right end, the known column vector is denoted by f after the combining operation. Then, (17) can be transformed into the following linear equation:(21)AX=f,where A is an N×N matrix and X and f are N×1 column vectors.SupposeU^=u^1,u^2,…,u^NT and Q^=q^1,q^2,…,q^NT are the numerical solutions of the potential and potential gradient obtained by the IBEFM. Since the IMLS method satisfies interpolation characteristics, it follows on the Dirichlet and Neumann boundaries that(22)u^I=u¯xI,I=1,2,…,N1,q^I=q¯xI,I=N1+1,N1+2,…,N.LetX^=q^1,…,q^N1,u^N1+1,…,u^NT be the solution of the linear system (21), that is,(23)AX^=f.Letuxs and qxs be the exact solutions to the potential and potential gradient, respectively. Also, let(24)Xe=qx1,…,qxN1,uxN1+1,…,uxNT.Therefore, at the boundary node, the following theorem exists between the exact solution and the numerical solution of the IBEFM.Theorem 1.
When the boundary is sufficiently smooth and the radius of the influence domain of the boundary node is small enough, there exists a constantC independent with ρ such that(25)Xe−X^∞≤Cku,qcondAf∞ρm+1,where ⋅∞ represents the maximum norm of a vector (matrix) and(26)ku,q=maxqΓ,m+1,uΓ,m+1maxqΓ,∞,uΓ,∞.Proof.
Let(27)uxs=uhxs+r1s,qxs=qhxs+r2s,where uhxs and qhxs are the approximation functions obtained by the IMLS method.
Then, from the error estimates of the IMLS method [28, 44, 45], there exist(28)r1sΓ,0≤Cρm+1uΓ,m+1,(29)r2sΓ,0≤Cρm+1qΓ,m+1.
It follows from (10) that(30)Cξuξ=∫Γqhx+r2su∗ξ,xdΓ−∫Γuhx+r1sq∗ξ,xdΓ−∫Ωbxu∗ξ,xdΩ,where ξ is the source point on the boundary Γ.
Then,(31)CξIuξI=∫Γ∑J=1nIΦJsqxJ+r2su∗ξI,xdΓ−∫Γ∑J=1nIΦJsuxJ+r1sq∗ξI,xdΓ−∫Ωbxu∗ξI,xdΩ,where ξI is a boundary node.
The same numerical integration is used for (16) and (31). Also, using the interpolation properties of the IMLS method, we have(32)B+HUe=GQe+ε,where(33)Ue=ux1,ux2,…,uxNT,Qe=qx1,qx2,…,qxNT,and ε is an N×1 column vector with the element given by(34)εI=∫Γr2su∗ξI,xdΓ−∫Γr1sq∗ξI,xdΓ.
By adopting the same shifting arrangement of (21) for (32), there exists(35)AXe=f+ε,where(36)Xe=qx1,…,qxN1,uxN1+1,…,uxNT.
According to matrix theory, it follows from (23) and (35) that(37)Xe−X^∞≤condAXe∞⋅ε∞f+ε∞,where ⋅∞ represents the maximum norm of the vector, e.g., ε∞=maxIεI.
From (28), (29), and (34), we have(38)εI≤∫Γr2s2dΓ1/2∫Γu∗ξI,x2dΓ1/2+∫Γr1s2dΓ1/2∫Γq∗ξI,x2dΓ1/2≤CmaxqΓ,m+1,uΓ,m+1ρm+1.
Then,(39)Xe−X^∞≤Cku,qcondAf+ε∞ρm+1,where(40)ku,q=maxqΓ,m+1,uΓ,m+1maxqΓ,∞,uΓ,∞.
When the Gaussian point and the source point do not coincide, it can be seen from equations (6) and (7) that the basic solutions are bounded. Also, from the errors of the IMLS method [28, 45], there exists constant C independent of ρ such that(41)ΦJs≤C.
Then, on the boundary, whenqJ or uJ is given, the expression u∗xI,x∑J=1nIΦJsqJ or q∗xI,x∑J=1nIΦJsuJ must be bounded. Thus, it can be seen from (34)–(36) that the element of the vector f must be a bounded value that is not a constant equal to 0. From (38), when the radius of the influence domain is small enough, the element value of ε is much smaller than that of f, which is almost negligible. Then, the theorem can be directly obtained from (39).
When the error caused by numerical integration is ignored, the following error estimate can be obtained.Theorem 2.
Supposeuξ is the solution of (1) determined by (9) and u^ξ is the numerical solution of the IBEFM solved from (19). When the boundary is sufficiently smooth and the radius ρ is small enough, there is a constant C independent of ρ such that(42)uξ−u^ξ≤Cku,qcondAf∞ρm+1.Proof.
From (19) and (9), we have(43)uξ−u^ξ=∫Γqx−∑J=1nIΦJsq^Ju∗ξ,xdΓ−∫Γux−∑J=1nIΦJsu^Jq∗ξ,xdΓ.
Then,(44)uξ−u^ξ=∫Γqx−qhx+qhx−∑J=1nIΦJsq^Ju∗ξ,xdΓ−∫Γux−uhx+uhx−∑J=1nIΦJsu^Jq∗ξ,xdΓ,=∫Γqx−qhx+∑J=1nIΦJsqJ−q^Ju∗ξ,xdΓ−∫Γux−uhx+∑J=1nIΦJsuJ−u^Jq∗ξ,xdΓ.
From the errors of the IMLS method [28, 45], there exists constant C independent of ρ such that(45)ΦJs≤C.
From (25), it follows that(46)qJ−q^J≤Cku,qcondAf∞ρm+1,(47)uJ−u^J≤Cku,qcondAf∞ρm+1.
Then, using trigonometric inequality and Hölder’s inequality and substituting (28), (29), (46), and (47) into (44), the theorem holds.
The value off∞ contained in Theorem 2 is determined by the value of the potential function and its potential gradient on the boundary. Therefore, when the solution of the problem (1) is sufficiently smooth, there are three main aspects that affect the errors of the numerical solution: the size of the influence domain of nodes, that is, the density of nodes; the order of the complete polynomials basis function of the IMLS method; and the condition number of the coefficient matrix A.


## 4. Numerical Examples

Two numerical examples will be given to verify the error estimation of the interpolated boundary element-free method in this section. The linear basis function is used in all calculations. The radius of influence domain isρ=2h, where h is the arc length between two adjacent nodes.Example 1.
Consider the following Laplace equation:(48)∇2T=0,x1,x2∈Ω,where Ω is a circular domain with boundary Γ: x12+x22=1. The essential boundary condition is(49)Tx1,x2=T¯,x1,x2∈Γ,x2≥0,and the derivative boundary condition is(50)∂∂nTx1,x2=q¯,x1,x2∈Γ,x2<0.
T¯andq¯ are the known functions determined by the following exact solution:(51)Tx1,x2=ex2sinx1.
When 30 regular nodes are used (as shown in Figure1), in Figure 2, the exact and numerical results on the inner circle with radii of 0.5, 0.7, and 0.9 are given. From the figure, it can be seen that the numerical solution obtained by the IBEFM has high accuracy. Figure 3 shows the absolute errors of the solutions obtained by the BEFM and IBEFM on the half circle with a radius 0.5, respectively, in which the BEFM and IBEFM adopt the same boundary integral grid. It can be seen that the solutions of the IBEFM have better accuracy, which also shows that the IBEFM has better advantages than the BEFM. The MATLAB programs of the BEFM and IBEFM run 100 times, respectively, and their average CPU time is 2.2403 s and 1.9519 s. This shows that the IBEFM has high computational efficiency.
To study the convergence of the IFEFM with respect to the radiusρ, when 10, 20, 30, 40, and 50 regular nodes are used, the corresponding absolute errors between the analytical and numerical solutions at points 0.7,0.5 and 0.3,0.9 are as shown in Figures 4 and 5, respectively. As can be seen from these figures, the solution of the IBEFM has a good convergence order. It has a better convergence order than the error estimate given by this paper, which also shows the correctness of the theorem.
To study the influence of the condition number of the coefficient matrix on the error, we apply the random function “rand” of MATLAB software to randomly generate 50 kinds of irregular node distribution with 30 nodes. Then, under each irregular node distribution, the IBEFM is used to solve the numerical solutions at internal points 0.7,0.5 and 0.3,0.9.
Figures6 and 7 show, respectively, the changes of the absolute errors at 0.7,0.5 and 0.3,0.9 for the different matrix condition numbers. As can be seen from the figure, although the error is not strictly worse as the condition number becomes more massive, on the whole, there is a great correlation between the condition number and the error. When the value of the condition number becomes significantly larger, the error will deteriorate. This means that when the condition number of the coefficient matrix becomes considerably larger and the error will become larger, which is consistent with the theoretical results in the paper.Figure 1
Nodes distribution for Example1.Figure 2
The exact and numerical results on the inner circle with radii of 0.5, 0.7, and 0.9.Figure 3
The absolute errors of the solutions obtained by the BEFM and IBEFM, respectively, on the half circle with 0.5 radius.Figure 4
The variation of the absolute error for different influence radii at0.7,0.5.Figure 5
The variation of the absolute error for different influence radii at0.3,0.9.Figure 6
The variation of the absolute error for differentcondA at 0.7,0.5.Figure 7
The variation of the absolute error for differentcondA at 0.3,0.9.Example 2.
Consider the Poisson equation with a square temperature field:(52)∇2Tx1,x2=2x1+x2−x12−x22,0≤x1≤1,0≤x2≤1.
The Dirichlet and Neumann boundary conditions are, respectively,(53)Tx1,x2=T¯,x1,x2∈Γ,x2=0orx1=1,∂∂nTx1,x2=q¯,x1,x2∈Γ,x2=1orx1=0,where T¯andq¯ are known functions obtained from the exact solution:(54)Tx1,x2=x1−x12x2−x22.
When 40 regular nodes are used (as shown in Figure8), Figure 9 shows the exact and numerical solutions on the lines x1=0.5 and x1=0.9. From this figure, it can be seen that the numerical solutions are very consistent with the analytical results, which also shows that the solution of the IBEFM has good accuracy.
In Figure10, the absolute errors of the numerical solutions obtained by BEFM and IBEFM on the line x1=0.9 are given. It can be seen that the solution of the IBEFM still has good accuracy, which also shows the superiority of the IBEFM. The MATLAB programs of the BEFM and IBEFM are repeatedly run 100 times, and their corresponding average CPU time is 3.0335 s and 3.0170 s, which also shows that the IBEFM has high computational efficiency.
In order to study the convergence of the IFEFM with respect to the radiusρ, when the boundary nodes distributions of 16, 24, 40, and 60 nodes are, respectively, used, the corresponding absolute errors for different influence domain radii at the interior point 0.4,0.7 and 0.5,0.1 are given in Figures 11 and 12, respectively. It can be seen from these figures that the solution of the IBEFM has a good convergence order, which is consistent with the theoretical results.Figure 8
Nodes distribution for Example2.Figure 9
The exact and numerical solutions on the linesx1=0.5 and x1=0.9.Figure 10
The absolute errors of the numerical solutions obtained by the BEFM and IBEFM on the linex1=0.9.Figure 11
The absolute errors for different influence radii at the interior point0.4,0.7.Figure 12
The absolute errors for different influence domain radii at the point0.5,0.1.


## 5. Conclusions

The IBEFM has both the advantages of the high accuracy of the boundary element method and meshless method, and it can also directly apply the boundary conditions. Based on the error estimation of the IMLS method and related matrix theory, this paper studies the error estimation of the IBEFM for two-dimensional potential problems. The relationship between the error and the influence radius and the condition number of the coefficient matrix is studied. The theoretical results show that the upper error limit of the IBEFM solution will decrease with the decrease in the influence domain radius of the boundary nodes and rise with the increase in the condition number of the coefficient matrix. Two numerical examples are provided. The numerical results are in good agreement with the theoretical results in this paper, which verifies the correctness of the theoretical results in the article.



---

*Source: 6378745-2020-07-25.xml*