A fourth order block-hexagonal grid approximation for the solution of Laplace’s equation with singularities

The hexagonal grid version of the block-grid method, which is a difference-analytical method, has been applied for the solution of Laplace’s equation with Dirichlet boundary conditions, in a special type of polygon with corner singularities. It has been justified that in this polygon, when the boundary functions away from the singular corners are from the Hölder classes C4,λ, 0 < λ < 1, the uniform error is of order O(h4), h is the step size, when the hexagonal grid is applied in the ‘nonsingular’ part of the domain. Moreover, in each of the finite neighborhoods of the singular corners (‘singular’ parts), the approximate solution is defined as a quadrature approximation of the integral representation of the harmonic function, and the errors of any order derivatives are estimated. Numerical results are presented in order to demonstrate the theoretical results obtained.


Introduction
It is well known that angular singularities arise in many applied problems when the solution of Laplace's equation is considered, and that finite-difference and finite-element methods may become less accurate when singularities are not taken into account. In the last two decades, for the solution of singularity problems, various combined and highly accurate methods have been proposed (see [-], and references therein).
Among these methods the block-grid method (BGM), presented in [-], on polygons with interior angles α j π , j = , , . . . , N , where α j ∈ {   , ,   , } (staircase polygons), requires the finite neighborhood of the singular vertices to be covered by sectors (blocks), and the remaining part of the domain by overlapping rectangles ('nonsingular' part). The finitedifference method with square grids is used for the approximate solution in the 'nonsingular' part, and in the blocks the integral representations of the harmonic function are approximated by the exponentially convergent mid-point quadrature rule (see []). Finally these subsystems are connected together by constructing an appropriate order matching operator. BGM is a highly accurate method not only for the approximation of the solution, but also for the approximation of its derivatives around singular points.
In this paper, the fourth order BGM is extended and justified for the Dirichlet problem of Laplace's equation on polygons with interior angles α j π , where α j ∈ {   ,   , , } (non-staircase), by gluing with the matching operator the -point approximation on a hexagonal grid in the 'nonsingular' part and the approximation of the integral representations around the singular points (on 'singular' parts). An advantage of using the hexagonal grid version of BGM in this domain is that a highly accurate approximation on the irregular grids is not required as in []. Thus the realization of the total system of algebraic equations becomes simpler. This may not be the case for this type of domain when square or rectangular grids are applied.
Furthermore it is justified that, when the boundary functions on the sides except the adjacent sides of the singular vertices are given in C ,λ ,  < λ < , the proposed hexagonal grid version of BGM has an accuracy of O(h  ), h is the mesh step. The same order of accuracy is obtained for the -point scheme on a square grid (see [, ]).
Finally in the last section of the paper, numerical experiments are demonstrated to support the theoretical results obtained.

Boundary value problem on a special type of polygon
Let D be an open simply connected polygon, γ j , j = , , . . . , N , be its sides, including the ends, enumerated counterclockwise (γ  ≡ γ N , γ  ≡ γ N+ ), and let α j π , α j ∈ {   ,   , , }, be the interior angles formed by the sides γ j- and γ j . Furthermore, letγ j = γ j- ∩ γ j be the jth vertex of D, γ = N j= γ j be the boundary of D; s is the arclength measured along the boundary of D in the positive direction, and s j is the value of s atγ j . We denote by r j , θ j the polar system of coordinates with poles inγ j and the angle θ j is taken counterclockwise from the side γ j .
Consider the boundary value problem where ≡ ∂  /∂x  + ∂  /∂y  , ϕ j , j = , , . . . , N , are given functions and In addition, at the verticesγ j , for α j = /, the following conjugation conditions are satisfied: No compatibility conditions are required at the vertices for α j = /. Moreover, it is required that when α j = /, the boundary functions on γ j- and γ j are given as algebraic polynomials of the arclength s measured along γ . Let E = {j : α j = /,  ≤ j ≤ N}. We construct two fixed block sectors in the neighborhood ofγ j , j ∈ E, denoted by T i j = T j (r ji ) ⊂ D, i = , , where  < r j < r j < min{s j+s j , s js j- }, T j (r) = {(r j , θ j ) :  < r j < r,  < θ j < α j π}. On the closed sector T  j , j ∈ E, we consider the carrier function Q j (r j , θ j ) in the form given in [], which satisfies the boundary conditions We set (see []) is the kernel of the Poisson integral for a unit circle.
The following lemma acts as a basis for the approximation of the solution around the verticesγ j with angle α j π , α j = /.

Lemma . ([]) The solution u of problem ()-() can be represented on T
where V j is the curvilinear part of the boundary of the sector T  j .
For the approximation of problem (), () in the domain D, we apply the hexagonal grid version of the block-grid method (see [-]). The application of this method first of all requires the construction of two more sectors T  j and T  j , where  < r j < r j < r j . Let D T = D\( j∈E T  j ). The following steps are taken for the realization: () We blockade the singular cornersγ j , j ∈ E, by the double sectors T i j (r ji ), i = , , with T  k ∩ T  l = ∅, k = l, k, l ∈ E, and cover the polygon D by overlapping parallelograms denoted by D l , l = , , . . . , M, and sectors T  j , j ∈ E, such that the distance from D l tȯ γ j is not less that r j for all l = , , . . . , M.
For the arrangement of the parallelograms D l , l = , , . . . , M, it is required that any point P lying on η l ∩ D T ,  ≤ l ≤ M, or lying on V j ∩ D, j ∈ E, falls inside at least on of the parallelograms D l(P) ,  ≤ l(P) ≤ M, depending on P, where the distance from P to D T ∩ η l(P) is not less than some constant κ  independent of P.
Let h ∈ (, κ  /] be a parameter, and define a hexagonal grid on D l ,  ≤ l ≤ M, with maximal positive step h l ≤ h, such that the boundary η l lies entirely on the grid lines. Let D lh be the set of grid nodes on D l , η h l be the set of nodes on η l , and let D lh = D lh ∪ η h l . Furthermore, η h l denotes the set of nodes on η l ∩ D T , η h l = η h l \η h l , and t h j denotes the set of nodes on t j . We also specify a natural number n ≥ [ln +χ h - ] + , where χ >  is a fixed number, and the quantities n(j) = max{, [α j n]}, β j = α j π/n(j), and θ m j = (m -/)β j , j ∈ E,  ≤ m ≤ n(j). On the arc V j we choose the points (r j , θ m j ),  ≤ m ≤ n(j), and denote the set of these points by V n j . Finally, we have For the approximation of the solution at the points of the set ω h,n we use the fourth order linear matching operator S  constructed in [], which can be represented as follows: Consider the system of difference equations where  ≤ l ≤ M, j ∈ E, and the operator S is defined as The solution of this system is the approximation of the solution of problem (), () on D h,n * .
Theorem . There is a natural number n  such that for all n ≥ n  the system of equations ()-() has a unique solution.
Proof Taking into account the corresponding homogeneous system of system ()-(), the proof follows by analogy to Lemma  in [].

Now consider the sector T
Let u h be the solution of the system of equations ()-(). The function

Error analysis of the 7-point approximation on the special parallelogram
Let D be one of the parallelograms covering the 'nonsingular' part of the polygon D defined in Section . The boundaries of the parallelogram D are denoted by γ j , enumerated counterclockwise starting from left, including the ends, We consider the system of finite-difference equations: where Since () has nonnegative coefficients and their sum is equal to , the solution of system (), () exists and is unique (see []).
Everywhere below we will denote constants which are independent of h and of the cofactors on their right by c, c  , c  , . . . , generally using the same notation for different constants for simplicity. Let D h,k be the set of nodes whose distance from the point P ∈ D h to γ h is

Lemma . Let
and z k h = const. be the solution of the system of equations Proof The proof follows analogously to the proof of the comparison theorem given in [].

Lemma . Let v be the trace of the solution of problem (), () on D h , and v h be the solution of system (), ().
If on the vertices with an interior angle of α j π = π/, j = , , , , then Let D h contain the set of nodes whose distance from the boundary γ is √ h  , and hence for (x, y) ∈ D h , (x + sH, y Moreover, let We rewrite problem (), () as In order to obtain an estimation for Svv on D h , we use Taylor's formula. On the basis of Lemma ., we have Since at least two values of  h in S  h are lying on the boundary γ h , on which  h = , from (), (), and the maximum principle (see []), we obtain where c  = c  . Next, we consider the estimation of  h . Let D h,k be the set of nodes whose distance from the point P ∈ D h to γ h is denotes the integer part of c, and d t is the minimum of the half-lengths of the sides of the parallelogram. Furthermore, D h, ≡ D h and D h, ≡ γ h . Since the vertices with α j =   of the parallelogram D are never used as a node of the hexagonal grid for the estimation of |Sv -v| on D h,k ,  ≤ k ≤ a * , we use the inequality for the sixth order derivatives, where ρ is the distance from (x, y) ∈ D to γ m . Hence, we obtain Consider a majorant function of the form Hence Y k is a solution of the finite-difference problem We represent the solution of system () as the sum of the solution of the following subsystems: where  ≤ k ≤ a * , μ k =  when k =  and |μ k | ≤ c  h +λ k -λ when k = , , . . . , a * . By (), (), and Lemma ., it follows that Hence, by taking () and () into consideration, we have On the basis of (), (), and (), we have estimation ().

Error analysis of the hexagonal block-grid equations
where u h is the solution of system ()-(), and u is the trace of the solution of problem (), () on D h,n * . It is easy to show that () satisfies the system of equations where  ≤ m ≤ N ,  ≤ l ≤ M, j ∈ E, and where j ∈ E,  ≤ q ≤ p, and c p , p = , , . . . , are constants independent of r j , θ j , and h.
Proof By taking estimation () into account, the proof follows by analogy to the proof of Theorem . in [].

Numerical results
Two examples have been solved in order to test the effectiveness of the proposed method. In Example ., it is assumed that there is a slit in the domain D, thus causing a strong singularity at the origin. The vertexγ  containing the singularity has an interior angle of α  π = π . The exact solution of this problem is assumed to be known. In Example ., we consider a problem with two singularities. The vertices which contain the singularities have interior angles of α j π =   π , j = , . In this example, the exact solution is not known. After separating the 'singular' part in Example ., the remaining part of the domain is covered by  overlapping parallelograms, whereas in Example ., the 'nonsingular' part of the domain is covered by only two parallelograms. For the solution of the block-grid equations, Schwarz's alternating method is used. In each Schwarz iteration the system of equations on the parallelograms are solved by the block Gauss-Seidel method. The carrier function is constructed for each example, taking into consideration the boundary conditions given on the adjacent sides of the vertices in the 'singular' parts. Furthermore, the derivatives are approximated in the 'singular' parts for both of the examples.
The results are provided in Tables - We assume that there is a slit along the straight line y = ,  ≤ x ≤ . Let γ j , j =      , , . . . , , be the sides of D, including the ends, enumerated counterclockwise starting from the upper side of the slit (γ  ≡ γ  ), γ =  j= γ j , andγ j = γ j ∩ γ j- be the vertices of D. Let (r, θ ) ≡ (r  , θ  ) be a polar system of coordinates with pole inγ  , where the angle θ is taken counterclockwise from the side γ  .
We consider the boundary value problem u =  on D, u = ϕ j on γ j , j = , , . . . , , where ϕ j is the value of the function v(r, θ ) = .r / sin θ  + .r / sin θ  + r  cos θ + .r  cos θ + θ on γ j .     As ϕ  = x  + .x  + π and ϕ  = x  + .x  , we obtain the carrier function in the form where ξ  (r, θ ) = r  ((πθ ) cos (πθ ) + ln r sin (πθ ))/π and ξ  (r, θ ) = r  ((πθ ) cos (πθ ) + ln r sin (πθ ))/π . The following notation is used in the table of results. Let D l , l = , , . . . , , be the open overlapping parallelograms, D NS =  l= D l be the 'nonsingular' part, and D S = D\D NS denote the 'singular' part of D (see Figure ). In Table , the values are obtained in the maximum norm of the difference between the exact and the approximate solutions, for the values of h =  -k , k = , , , , and n, which is the number of quadrature nodes on V j . The order of convergence, has also been included. We also present the error obtained between the derivatives of the exact and the block-grid solutions ( , in the maximum norm, in Tables  and , respectively. Figures  and  illustrate the shapes of the derivatives ∂u ∂x and ∂  u ∂x  of the approximate (a) and the exact (b) solutions. These figures also demonstrate the highly accurate approximation of the derivatives.

Example . Let P be the open parallelogram
 }, let γ j , j = , , , , be the sides of P, including the ends, enumerated counterclockwise starting from left (γ  ≡ γ  , γ  ≡ γ  ), γ =  j= γ j , andγ j = γ j ∩ γ j- be the vertices of P. We look at a problem with two corner singularities at the verticesγ  andγ  , where α j π =   π , j = , . The two 'singular' corners of P are covered by sectors and these areas are denoted by P i S , i = , , and two overlapping parallelograms cover the 'nonsingular' part of the domain, denoted by P i NS , i = ,  (see Figure ). We consider the boundary value problem u =  on P, u =  on γ j , j = , , () u =  on γ j , j = , .
The carrier functions constructed for each singularity are Q  (r  , θ  ) =  -θ  π and Q  (r  , θ  ) = θ  π . We have checked the accuracy of the obtained approximate results u h by looking at the order of convergence using the formula R m P = u  -m -u  -m+ P u  -m- -u  -m P , which corresponds to   , for the pairs (h, n) = ( - , ), ( - , ), ( - , ), ( - , ). The results are presented in Table . Moreover, ∂  u ∂x  has been approximated in the 'singular' part, where u is the unknown exact solution of problem (). The results are presented in Table  and illustrated further in Figure .

Conclusion
A fourth order square and hexagonal grid version of the block-grid method, for the solution of the boundary value problem of Laplace's equation on staircase polygons, with interior angles α j ∈ {   , ,   , }, is extended for the polygons with interior angles α j π , α j ∈ {   ,   , , }, by constructing and justifying the block-hexagonal grid method. Moreover, the smoothness requirement on the boundary functions away from the singular vertices (outside of the 'singular' parts) is lowered down from the Hölder classes C ,λ ,  < λ < , as in [], to C ,λ ,  < λ < , which was proved for the -point scheme on square grids (see [, ]).
The proposed version of the BGM can be applied for the mixed boundary value problem of Laplace's equation on the above mentioned polygons. Furthermore, by this method any order derivatives of the solution can be highly approximated on the 'singular' parts, which are difficult to obtain in other numerical methods.
This method can also be used for the solution of the biharmonic equation by representing the problem with two problems for the Laplace and Poisson equations.