Localized Boundary Knot Method for Solving Two-Dimensional Laplace and Bi-Harmonic Equations

In this paper, a localized boundary knot method is proposed, based on the local concept in the localized method of fundamental solutions. The localized boundary knot method is formed by combining the classical boundary knot method and the localization approach. The localized boundary knot method is truly free from mesh and numerical quadrature, so it has great potential for solving complicated engineering applications, such as multiply connected problems. In the proposed localized boundary knot method, both of the boundary nodes and interior nodes are required, and the algebraic equations at each node represent the satisfaction of the boundary condition or governing equation, which can be derived by using the boundary knot method at every subdomain. A sparse system of linear algebraic equations can be yielded using the proposed localized boundary knot method, which can greatly reduce the computer time and memory required in computer calculations. In this paper, several cases of simply connected domains and multi-connected domains of the Laplace equation and bi-harmonic equation are demonstrated to evidently verify the accuracy, convergence and stability of this proposed meshless method.


Introduction
With the emergence of various engineering problems, an increasing number of numerical methods have been proposed in the past decades. To date, there are some well-known numerical methods, such as the finite difference method (FDM) [1,2], the finite element method (FEM) [3] and the boundary element method (BEM) [4]. These three methods are relatively mature, but they still have some obvious shortcomings and limitations. For example, with these methods, grid generation is difficult, especially for problems in complicated domains. They need to constantly remesh when calculating specific problems, which may result in high computational costs. On the contrary, there are some so-called meshless methods, which can solve problems without time-consuming mesh generation. From the present studies of meshless methods, it is evidently demonstrated that the meshless methods, proposed recently, are truly free from mesh and numerical quadrature, so they can be adopted to efficiently analyze physical problems and engineering applications with complex geometries.
In recent years, various meshless or meshfree numerical methods have been proposed to eliminate the time-consuming task of meshing. Meshless methods usually require node information only and do not need any orthogonal grid or unstructured mesh to acquire convergent, accurate and stable numerical solutions. These newly developed meshless methods have attracted the attention of many scholars [5]. Among them, the Radial basis function (RBF)-related method is the most popular, which

Mathematical Formulation of Laplace and Bi-Harmonic Equations
In this study, the numerical procedures of the localized BKM was proposed to solve the two-dimensional boundary value problem, which are governed by Laplace and bi-harmonic equations, respectively. The governing equation of Laplace is demonstrated as follows: where ∇ 2 = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 is the two-dimensional Laplacian, u(x, y) is the unknown variable in the computational domain Ω, Γ D and Γ N denote the boundary segments with Dirichlet boundary conditions and the Neumann boundary condition: where n = n x , n y is the unit outward normal vector on the boundary and h = (x, y) and g = (x, y) are the given boundary conditions, respectively. The other boundary value problem is the bi-harmonic equation: The Neumann boundary condition, Equation (2) and the Dirichlet boundary condition, Equation (3), should be simultaneously applied to all of the boundaries, since the bi-harmonic equation is a fourth-order PDE. In the next section, the numerical procedures of the proposed localized BKM will be clearly described.

Laplace Equation
To describe the numerical procedures of the proposed localized BKM, for a two-dimensional Laplace equation, we required N = n i + n b1 + n b2 computational nodes, where n i is the number of the nodes inside the computational domain, n b1 and n b2 are the points along the whole boundary, Γ D and Γ N . The schematic diagram of the computational nodes is displayed in Figure 1a.
where ( ) The Neumann boundary condition, Equation (2) and the Dirichlet boundary condition, Equation (3), should be simultaneously applied to all of the boundaries, since the bi-harmonic equation is a fourth-order PDE. In the next section, the numerical procedures of the proposed localized BKM will be clearly described.

Laplace Equation
To describe the numerical procedures of the proposed localized BKM, for a two-dimensional Laplace equation, we required  To solve the two-dimensional boundary value problems, a subdomain, which is depicted in Figure 1b, will be formed for each node and the numerical solutions within each subdomain can be approximated by  To solve the two-dimensional boundary value problems, a subdomain, which is depicted in Figure 1b, will be formed for each node and the numerical solutions within each subdomain can be approximated by where N s is the number of nodes in each subdomain. The determination of subdomain is identical to the LRBFCM, α s j are the unknown coefficients, r = x 2 + y 2 , x = x − x s , y = y − y s , is the Euclidean distance and x s represents the sth local node in the neighborhood of x. The nonsingular general solution ξ(x, y) = e (−c(x 2 −y 2 )) cos(2cxy) satisfies the two-dimensional Laplace equation and is regular inside the computational domain. c is the shape parameter, which is adopted in the range 0 < c ≤ 3 in this paper.
From Equation (5)-and similarly to the LMFS [18][19][20][21]-the unknown coefficients α j s can be transformed as follows: where By substituting Equation (6) into Equation (1), the following expression can be derived: where L is the differential operator , which is a vector with size 1 × N s and we use the MATLAB built-in command pinv to calculate ξ −1 is an ill-conditioned matrix. Furthermore, from Equation (2), we have In order to derive the expression for the Neumann boundary conditions, we can substitute Equations (8) and (9) into Equation (2): where are the interpolation vectors of the discrete point satisfying the Neumann boundary conditions ψ N s = h x( j) s n x + h y( j) s n y and B are the differential operators in Γ N . Equations (7) and (10) are the local forms, which is different from the global method. Now, we need to globalize this local form by inserting zero elements into the calculation. Therefore, Equations (7) and (10) can be transformed as follows: Bu(x, y) = Bψu = h(x, y), (x, y) ∈ Γ N , (12) u(x, y) = ψu = g(x, y), (x, y) ∈ Γ D , where ψ(x, y) = ψ 1 (x, y), ψ 2 (x, y), . . . , ψ N (x, y) is a sparse vector with the size of 1 × N, and N is the total number of nodes. All the values are zeros, except for the values obtained by interpolation at the corresponding positions of each discrete point in the local domain. Equation (13) represents the satisfaction of the Dirichlet boundary condition. The discretizations of the Laplace equation, the Neumann boundary condition and the Dirichlet boundary condition are combined to form a matrix as follows: where K N×N is the sparse coefficient matrix that avoids the ill-conditioned matrix, u = [u 1 , u 2 , . . . , u N ] T is the unknown field quantity and f is the given condition. Finally, u can be efficiently computed from Equation (14). From the descriptions of the numerical procedures, the proposed localized BKM, which is the combination of the conventional BKM and the localization concept from the LMFS and the LRBFCM, is very simple and easy to program. Furthermore, it can be expected that the localized BKM can be applied to engineering applications in complicated domains owing to the formation of a sparse system of linear algebraic equations.

Bi-Harmonic Equation
The numerical procedures of the proposed localized LBKM for solving the two-dimensional bi-harmonic equation are presented in this subsection. The numerical procedures of the localized BKM for solving the bi-harmonic equation is similar to the procedures for the Laplace equation described above. As compared with the Laplace equation, the bi-harmonic equation is a fourth-order PDE. In addition to the bi-harmonic equation, there are two boundary conditions, which should be imposed along the whole boundary at the same time. The process of the calculation is basically the same as that for the Laplace equation. The numerical solution within the local domain of the jth node can be represented as follows: where α j s and β j s are the unknown coefficients, and N s is the number of nodes in the subdomain. For simplicity, Equation (15) can be simply transformed as follows: where C is the coefficient matrix based on the jth point and the neighboring N s points through interpolation. ω = [α 1 , α 2 , . . . , α N s , β 1 , β 2 , . . . , β N s ] T is the vector of unknown coefficients. Similarly, from Equations (7)-(9), we can reformulate the governing equation and the boundary conditions as follows: where η ( j) n s , η x( j) n s and η y( j) n s are the weighting coefficients. The satisfaction of the governing equation is imposed at every interior node while the satisfactions of both boundary conditions, Equations (17)- (19), are enforced at each boundary node.
To combine Equations (20)- (22), we can acquire a matrix form of the linear algebraic equation where K is the matrix with the size (n i + 2n b ) × (n i + n b ). Once the above system is accurately solved, the numerical solution of the proposed localized BKM can be obtained. However, it is also known that the numerical solutions may not be accurate enough as a result of the overdetermined system. In order to acquire a highly accurate solution, we adopted the ghost nodes, which is shown in Figure 2. As shown in Figure 2, we place a circle of ghost nodes around the boundary nodes and the distance between the ghost nodes and the boundary nodes is 1/2 the distance between the adjacent boundary nodes. When the base nodes are interpolated with its adjacent n s , the adjacent ghost nodes are also selected to participate in the interpolation calculation. When K is calculated in this way, we can obtain a square coefficient matrix and an extremely accurate solution.   (17)- (19), are enforced at each boundary node.
To combine Equations (20)- (22), we can acquire a matrix form of the linear algebraic equation , . Once the above system is accurately solved, the numerical solution of the proposed localized BKM can be obtained. However, it is also known that the numerical solutions may not be accurate enough as a result of the overdetermined system. In order to acquire a highly accurate solution, we adopted the ghost nodes, which is shown in Figure  2. As shown in Figure 2, we place a circle of ghost nodes around the boundary nodes and the distance between the ghost nodes and the boundary nodes is 1/2 the distance between the adjacent boundary nodes. When the base nodes are interpolated with its adjacent ns, the adjacent ghost nodes are also selected to participate in the interpolation calculation. When K is calculated in this way, we can obtain a square coefficient matrix and an extremely accurate solution.

Numerical Results and Comparisons
In this section, six numerical cases are given and the numerical results of the LBKM are compared. In order to prove the accuracy of the proposed algorithm, the simply connected domain and multi-connected domain are calculated and analyzed by changing the parameters, such as the shape parameters, the total number of nodes or the number of subdomains. Finally, the calculated numerical solution u n is compared with the analytical solution u e , and the maximum absolute error is compared as an index. The calculation formula is as follows:

Case 1
In this case, we test a square domain, whose four boundaries are applied to the Dirichlet boundary conditions. A detailed diagram of the computational domain and the applied boundary conditions are shown in Figure 3. The analytical solution of the applied boundary condition is as follows: where N is the total nodes, n b is the boundary nodes, n s is the number of subdomains, and c is the shape parameter. The maximum absolute error is 5.3892 × 10 −5 , which can preliminarily prove the accuracy of LBKM. The distribution of the numerical solution is shown in Figure 4.

Numerical Results and Comparisons
In this section, six numerical cases are given and the numerical results of the LBKM are compared. In order to prove the accuracy of the proposed algorithm, the simply connected domain and multi-connected domain are calculated and analyzed by changing the parameters, such as the shape parameters, the total number of nodes or the number of subdomains. Finally, the calculated numerical solution n u is compared with the analytical solution e u , and the maximum absolute error is compared as an index. The calculation formula is as follows: In this case, we test a square domain, whose four boundaries are applied to the Dirichlet boundary conditions. A detailed diagram of the computational domain and the applied boundary conditions are shown in Figure 3. The analytical solution of the applied boundary condition is as follows:   To verify the stability of the LBKM, we changed the number of boundary nodes and interior nodes at the same time. As can be seen in the results shown in Table 1, the results calculated using different parameters are considerable, the shape parameter c remains unchanged, and the maximum  To verify the stability of the LBKM, we changed the number of boundary nodes and interior nodes at the same time. As can be seen in the results shown in Table 1, the results calculated using different parameters are considerable, the shape parameter c remains unchanged, and the maximum absolute error is 8.7091 × 10 −6 when the total number of nodes reaches N = 62997.

Case 2
In order to further verify the impact of geometric complexity on the precision of the algorithm. In this case, a multiply connected domain is adopted as the computational domain. As shown in Figure 5, we apply Dirichlet boundary conditions and Neumann boundary conditions on each boundary, respectively. Here, the analytical solutions of the applied Dirichlet boundary conditions are as follows: u e (x, y) = x 2 − y 2 + xy + 2. In this case, the shape parameter is 0.2 c = . The maximum absolute error can reach  Figure 6 shows the distribution of the numerical solution. In addition, the results computed by different parameters are shown in Table 2. It can be seen that when the number of nodes in the computational domain increases slowly and the number of nodes in the subdomain changes, the results of the calculation are convergent. Therefore, the error of LBKM can reach the accuracy requirement when calculating the multi-domain. In this case, the shape parameter is c = 0.2. The maximum absolute error can reach 5.1431 × 10 −4 , when the total number of nodes is N = 6762, the boundary number is n b = 220, and the number of the nodes in the local domain is 29. Figure 6 shows the distribution of the numerical solution. In addition, the results computed by different parameters are shown in Table 2. It can be seen that when the number of nodes in the computational domain increases slowly and the number of nodes in the subdomain changes, the results of the calculation are convergent. Therefore, the error of LBKM can reach the accuracy requirement when calculating the multi-domain.  Figure 6 shows the distribution of the numerical solution. In addition, the results computed by different parameters are shown in Table 2. It can be seen that when the number of nodes in the computational domain increases slowly and the number of nodes in the subdomain changes, the results of the calculation are convergent. Therefore, the error of LBKM can reach the accuracy requirement when calculating the multi-domain.

Case 3
For the third case, we simulate a simple connected domain to verify the accuracy of the solutions for the Laplace problem with non-harmonic boundary conditions. Furthermore, the analytic solution that we have is u e (x, y) = x 2 y 3 + 1, which does not satisfy the Laplace equation. Figure 7 shows the schematic diagram of the computational domain. Table 3 shows the maximum absolute error for case 3. It can be seen from Table 3 that when the total number of nodes increases and the number of nodes in the local domain changes, the calculated result is convergent. When the total number of nodes is N = 25064, the number of boundary nodes is n b = 600, the number of subdomain nodes n s = 18, and the shape parameter is c = 1, the maximum absolute error is 5.2868 × 10 −5 .

Case 3
For the third case, we simulate a simple connected domain to verify the accuracy of the solutions for the Laplace problem with non-harmonic boundary conditions. Furthermore, the analytic solution that we have is 1 , e u x y x y = + which does not satisfy the Laplace equation. Figure 7 shows the schematic diagram of the computational domain. Table 3 shows the maximum absolute error for case 3. It can be seen from Table 3 that when the total number of nodes increases and the number of nodes in the local domain changes, the calculated result is convergent. When the total number of nodes is 25064 N = , the number of boundary nodes is     In the above three cases, the new LBKM is used to solve the boundary value problem with the interior nodes satisfying the Laplace equation. By changing the total number of nodes N, the number of nodes n s in local domain and the shape parameter c, the convergence of the results is obtained when calculating the simply connected domain and the multi-connected domain, respectively, which preliminarily verifies the accuracy of the algorithm. In the cases below, the bi-harmonic equations are tested to further verify the stability and accuracy of the method.

Case 4
In this case, we make a simulation of a governed bi-harmonic equation with two types of boundary conditions. There are two tests in this case. The calculation domain is shown in Figure 8.  Table 4 shows the results calculated with different parameters. The maximum absolute error is  Figure 9.  For the first test, the analytic solution we use is u e (x, y) = x 4 + y 4 − 6x 2 y 2 + 3, which satisfies both the Laplace equation and the bi-harmonic equation. Table 4 shows the results calculated with different parameters. The maximum absolute error is 1.7719 × 10 −4 when the total number of nodes is N = 2336, the number of boundary nodes is n b = 296, the number of local domain nodes is n s = 41, and the shape parameter is c = 3. The distribution of the numerical solution is shown in Figure 9.    The other test uses the same computational domain and the same boundary conditions. However, the boundary condition can be calculated from the analytical solution as follows: where r = x 2 + y 2 is the Euclidean distance. The origin of the cylindrical coordinate (r s , θ s ) is located at (x, y) = (8, −4); Table 5 shows the maximum absolute error for this test. The distribution of the numerical solution is shown in Figure 10 for when the total number of nodes N = 24896, the number of boundary nodes n b = 776, the number of nodes in the subdomain n s = 26, and the shape parameter is c = 1. The other test uses the same computational domain and the same boundary conditions. However, the boundary condition can be calculated from the analytical solution as follows:

Case 5
In the fifth case, a multi-connected domain is considered. The computational domain is a square domain with a small square domain dug out, as shown in Figure 11.

Case 5
In the fifth case, a multi-connected domain is considered. The computational domain is a square domain with a small square domain dug out, as shown in Figure 11. The analytical solution is as follows:  Table 6 shows the results for when we change the number of total nodes, N , and the number of different nodes in the subdomain, ns . We can see from Table 6 Figure 12 is the distribution of the numerical solution.   The analytical solution is as follows: which only satisfies the bi-harmonic equation, but not the Laplace equation. The shape parameter in this case is c = 0.5. When the parameters are taken as follows: N = 14552, n b = 720, n s = 49, the maximum absolute error is 1.0615 × 10 −3 . Table 6 shows the results for when we change the number of total nodes, N, and the number of different nodes in the subdomain, ns. We can see from Table 6 that the result converges and when the total number of nodes is N = 27712, the number of boundary nodes is n b = 1000 and n s = 65, the error can reach 9.7262 × 10 −4 . Figure 12 is the distribution of the numerical solution.  The analytical solution is as follows:  Table 6 shows the results for when we change the number of total nodes, N , and the number of different nodes in the subdomain, ns . We can see from Table 6 Figure 12 is the distribution of the numerical solution.

Case 6
In the last case, we use the non-harmonic boundary conditions, governed by the bi-harmonic equation. Figure 13 shows the computational domain and the applied boundary conditions. The analytical solution we use is as follows: u e (x, y) = x 3 y 4 + 1.

Case 6
In the last case, we use the non-harmonic boundary conditions, governed by the bi-harmonic equation. Figure 13 shows the computational domain and the applied boundary conditions. The analytical solution we use is as follows: 1 .
e u x y x y = + Figure 13. Schematic diagram of computational domain in case 6.
In this case, the shape parameter 0.3 c = is used. Since the non-harmonic is mainly affected by the boundary nodes, the number of nodes in the interior nodes and local domain remains unchanged, only the number of boundary nodes is changed. Table 7 shows the numerical results calculated at different boundary nodes.

Conclusions
In this paper, the localized BKM is proposed to solve two-dimensional Laplace and bi-harmonic equations accurately and efficiently. The proposed meshless method is the combination of the convectional BKM and the concept of localization from the LMFS and the LRBFCM. The proposed method is truly free from time-consuming mesh generation and numerical quadrature since only a set of randomly distributed nodes are required for the numerical simulation. Furthermore, the troublesome problems of determination of the fictitious boundary for sources can be simply avoided by using the nonsingular general solution. Since the resultant system of algebraic equations is sparse, it can be expected that the proposed localized BKM can be used to solve realistic engineering applications in complicated domains.
In this paper, six numerical examples are provided to demonstrate the merits of the proposed localized BKM. According to the numerical results and comparisons, the following conclusions can be drawn: In this case, the shape parameter c = 0.3 is used. Since the non-harmonic is mainly affected by the boundary nodes, the number of nodes in the interior nodes and local domain remains unchanged, only the number of boundary nodes is changed. Table 7 shows the numerical results calculated at different boundary nodes.

Conclusions
In this paper, the localized BKM is proposed to solve two-dimensional Laplace and bi-harmonic equations accurately and efficiently. The proposed meshless method is the combination of the convectional BKM and the concept of localization from the LMFS and the LRBFCM. The proposed method is truly free from time-consuming mesh generation and numerical quadrature since only a set of randomly distributed nodes are required for the numerical simulation. Furthermore, the troublesome problems of determination of the fictitious boundary for sources can be simply avoided by using the nonsingular general solution. Since the resultant system of algebraic equations is sparse, it can be expected that the proposed localized BKM can be used to solve realistic engineering applications in complicated domains.
In this paper, six numerical examples are provided to demonstrate the merits of the proposed localized BKM. According to the numerical results and comparisons, the following conclusions can be drawn: (1) In this paper, a novel meshless method, the localized BKM, is proposed to accurately and efficiently solve two-dimensional Laplace and bi-harmonic equations; (2) As a result of the use of localization, the resultant system of algebraic equations is sparse, so it is evident that the proposed localized BKM is capable of efficiently solving large-scale problems; (3) As compared with the MFS, the RBFCM and other meshless methods, the problems of ill-conditioned matrix and the determination of the fictitious boundary for sources are avoided in the proposed method; (4) In the examples provided, it is evident that the proposed localized BKM can accurately solve problems in simply connected and multiply connected domains. Furthermore, both the Laplace equation with a non-harmonic condition and the bi-harmonic equation with a non-bi-harmonic conditions can be stably analyzed by the proposed meshless method.
In the future, the proposed localized BKM will be improved and extended to analyze various physical and mathematical problems, such as temporal transient problems, inverse problems, the acoustic problem, the moving-boundary problem, three-dimensional problems, etc.