Skip to main content

Approximation on the hexagonal grid of the Dirichlet problem for Laplace’s equation

Abstract

The fourth order matching operator on the hexagonal grid is constructed. Its application to the interpolation problem of the numerical solution obtained by hexagonal grid approximation of Laplace’s equation on a rectangular domain is investigated. Furthermore, the constructed matching operator is applied to justify a hexagonal version of the combined Block-Grid method for the Dirichlet problem with corner singularity. Numerical examples are illustrated to support the analysis made.

MSC:35A35, 35A40, 35C15, 65N06, 65N15, 65N22, 65N99.

1 Introduction

When solving PDEs, many approximate methods, such as overlapping versions of domain decomposition, composite grids and different types of combined methods, use the matching operator to connect the subsystems within them. Hence, the approximation of the solutions relies heavily on the order of accuracy of the matching operator, as well as on the order of accuracy of the subsystems. In [1] and [2], the second order matching operator is used to construct and justify the second order composite grid method for solving Laplace’s boundary value problems. In [3], the fourth order matching operator is constructed and used for the fourth order composite grids and in [4] and [5] it is used for the fourth order Block-Grid method. In [69], the sixth order matching operator is constructed for the Block-Grid method and it is used for the sixth order composite grids in [10].

In all of the above mentioned papers, the fourth and sixth order matching operators were constructed on the basis of the 9-point finite difference solution of Laplace’s equation on square grids. In this paper, the matching operator is constructed for the solution of the Dirichlet problem on a hexagonal grid. In order to approximate the given differential equations at each regular node P 0 on a hexagonal grid, the six equidistant nodes surrounding P 0 are used, and the truncation error obtained is O( h 4 ). Thus, we obtain the same order of accuracy when using the 7-point scheme on the hexagonal grid, as we do when using the 9-point scheme on the rectangular grid (see [11]). This has many computational advantages such as (i) the matrix of the system will contain seven diagonals rather than nine and will lead to less use of memory space, (ii) the calculations will require less computational effort and (iii) the algorithm will be easier to implement. Hexagonal grids are favored in many applied problems in dynamical meteorology and dynamical oceanography as well (see [1214]), due to the benefits a hexagonal grid provides, compared to a rectangular grid.

Even though using a hexagonal grid has the above mentioned advantages, it has not been used before in methods such as composite grids, domain decomposition, and combined methods, as the fourth order matching operator for connecting the subsystems together was not constructed. In this paper, in Section 2, the approximate solution in a hexagonal grid on a rectangular domain is analyzed. In Section 3 a fourth order matching operator is constructed and its application to find the fourth order accurate approximate solution on the closed domain is considered. Section 4 contains the justification of using a hexagonal grid for the solution of Laplace’s equation on a staircase polygon, with the use of the Block-Grid method. This method requires the application of the matching operator constructed in Section 3 and gives an overall fourth order accuracy. Numerical examples are illustrated in Section 5 to support the analysis made.

2 Approximation in the hexagonal grid of the Dirichlet problem on a rectangle

Let Π={(x,y):0<x<a,0<y<b} be a rectangle, γ j , j=1,2,3,4, be its sides, including the ends, enumerated counterclockwise starting from left ( γ 0 γ 4 , γ 1 γ 5 ), γ= j = 1 4 γ j be the boundary of Π, and let A j = γ j 1 γ j be the j th vertex. We consider the boundary value problem

Δu=0on Π,
(2.1)
u= φ j on  γ j ,j=1,2,3,4,
(2.2)

where Δ= 2 / x 2 + 2 / y 2 , φ j is a given function of arclength s taken along γ, and

φ j C 6 , λ ( γ j ),0<λ<1,j=1,2,3,4.
(2.3)

At the vertices s= s j ( s j is the beginning of γ j ), the conjugation conditions

φ j ( 2 q ) ( s j )= ( 1 ) q φ j 1 ( 2 q ) ( s j ),q=0,1,2,3,
(2.4)

are satisfied.

Let h>0, with a/h2, b/ 3 h2 integers. We assign Π h a hexagonal grid on Π, with step size h, defined as the set of nodes

Π h = { ( x , y ) Π : x = k j 2 h , y = 3 ( k + j ) 2 h , k = 1 , 2 , ; j = 0 ± 1 ± 2 , } .
(2.5)

Let γ j h be the set of nodes on the interior of γ j , and let γ ˙ j h = γ j γ j + 1 , γ h =( γ j h γ ˙ j h ), Π ¯ h = Π h γ h . Also let Π h denote the set of nodes whose distance from the boundary γ of Π ¯ is h 2 and Π 0 h = Π h Π h .

We consider the system of finite difference equations

u h =S u h on  Π 0 h ,
(2.6)
u h = S j u h + E j h ( φ j )on  Π h ,
(2.7)
u h = φ j on  γ j h ,j=1,2,3,4,
(2.8)

where

S u ( x , y ) = 1 6 ( u ( x + h , y ) + u ( x + h 2 , y + 3 2 h ) + u ( x h 2 , y + 3 2 h ) + u ( x h , y ) + u ( x h 2 , y 3 2 h ) + u ( x + h 2 , y 3 2 h ) ) ,
(2.9)
S j u ( x , y ) = 1 7 ( u ( x + h 2 , y 3 h 2 ) + u ( x + h , y ) + u ( x + h 2 , y + 3 h 2 ) ) ,
(2.10)
E j h ( φ j )= 1 21 ( 2 φ j ( y + 3 h 2 ) + 8 φ j ( y ) + 2 φ j ( y 3 h 2 ) ) .
(2.11)

From formulae (2.9) and (2.10) it follows that the coefficients of the expressions Su(x,y) and S j u(x,y) are non-negative, and their sums do not exceed one. Hence, on the basis of the maximum principle, it follows that the solution of system (2.6)-(2.8) exists and it is unique (see [11]).

Everywhere below we will denote constants which are independent of h and of the cofactors on their right by c, c 0 , c 1 , , generally using the same notation for different constants for simplicity.

Lemma 2.1 Let

v 1 = S v 1 + f h on  Π 0 h , v 1 = S j v 1 on  Π h , v 1 = 0 on  γ h ,

and

v 2 = S v 2 + f ¯ h on  Π 0 h , v 2 = S j v 2 + f ¯ h on  Π h , v 2 = η ¯ h on  γ h ,

where f h , f ¯ h , f ¯ h and η ¯ h are arbitrary grid functions. If the conditions

f ¯ h 0,| f h | f h ¯ ,and η ¯ h 0

are satisfied, then

| v 1 | v 2 .

Proof The proof of this lemma is similar to the proof of the comparison theorem (see Ch. 4 in [11]). □

Theorem 2.2 Let u be the solution of problem (2.1), (2.2) and u h be the solution of system (2.6)-(2.8), then

max Π ¯ h | u h u|c h 4 .
(2.12)

Proof Let

ϵ h = u h u,

where u is the trace of the solution of problem (2.1), (2.2) on Π ¯ h , and u h is the solution of system (2.6)-(2.8). Then, the error function ϵ h satisfies the following system:

ϵ h =S ϵ h + Ψ h on  Π 0 h ,
(2.13)
ϵ h = S j ϵ h + Ψ h on  Π h ,
(2.14)
ϵ h =0on  γ h ,
(2.15)

where

Ψ h =Suu,
(2.16)
Ψ h = S j uu+ E j h ( φ j )
(2.17)

are the truncation errors of equations (2.6) and (2.7), respectively.

On the basis of conditions (2.3) and (2.4), from Theorem 3.1 in [15] it follows that u C 6 , λ ( Π ¯ ), 0<λ<1. Then, by Taylor’s formula, we obtain (see [16])

max ( x , y ) Π ¯ | Ψ h (x,y)| c 1 h 6 M 6 ,
(2.18)

where

M j = sup ( x , y ) Π ¯ { | j u ( x , y ) x i y j i | , i = 0 , 1 , , j } .

We represent the solution of (2.13)-(2.15) as

ϵ h = ϵ h 1 + ϵ h 2 ,
(2.19)

where

ϵ h 1 = S ϵ h 1 + Ψ h on  Π 0 h , ϵ h 1 = S j ϵ h 1 on  Π h , ϵ h 1 = 0 on  γ h ,

and

ϵ h 2 = S ϵ h 2 on  Π 0 h , ϵ h 2 = S j ϵ h 2 + Ψ h on  Π h , ϵ h 2 = 0 on  γ h .

By taking the function as v 2 = h 4 c 1 M 6 ( a 2 + b 2 x 2 y 2 ) in Lemma 2.1, we obtain

max ( x , y ) Π h | ϵ h 1 | max ( x , y ) Π ¯ | v 2 | c 2 h 4 M 6 .
(2.20)

Using Taylor’s formula about each of the points ( h 2 ,y) Π h and from (2.17), we have

max ( x , y ) Π h | Ψ | c 3 M 4 h 4 .

On the basis of the maximum principle, we obtain

max ( x , y ) Π h | ϵ h 2 | 7 4 max ( x , y ) Π h | Ψ h | c 4 M 4 h 4 .
(2.21)

From (2.19), (2.20), and (2.21) it follows that

max ( x , y ) Π ¯ h | ϵ h |c h 4 .
(2.22)

 □

Remark 2.3 Estimation (2.12) remains true when E j h ( φ j ) in system (2.6)-(2.8) is replaced by

E j , h = 1 7 φ j ( y 3 h 2 ) + 1 7 φ j ( y + 3 h 2 ) + 2 7 φ j (y) h 2 2 ! 14 φ j ( 2 ) (y)+ h 4 4 ! 56 φ j ( 4 ) (y).

3 Construction of the fourth order matching operator

Let z=x+iy be a complex variable, and let Ω={z:|z|<1} be a unit circle. For a harmonic function u on Ω with u C 4 , 0 ( Ω ¯ ), by Taylor’s formula, any point (x,y)Ω can be represented as

u(x,y)= k = 0 3 a k Re z k + k = 1 3 b k Im z k +O ( r 4 ) ,
(3.1)

where r= x 2 + y 2 ,

a 0 = u ( 0 , 0 ) , a 1 = u ( 0 , 0 ) x , a 2 = 1 2 2 u ( 0 , 0 ) x 2 , a 3 = 1 3 ! 3 u ( 0 , 0 ) x 3 , b 1 = u ( 0 , 0 ) y , b 2 = 1 2 2 u ( 0 , 0 ) x y , b 3 = 1 3 ! 3 u ( 0 , 0 ) x 2 y .

By analogy with the idea used in [8], we construct the operator S 4 from the condition that the expression

S 4 u= ξ k u k ,

where u k =u( P k ), P k is a node of the hexagonal grid Π h , gives the exact value of any harmonic polynomial

F 3 (x,y)= k = 0 3 a k Re z k + k = 1 3 b k Im z k ,

at each point PΠ, and

ξ k 0, ξ k 1.

Let Π 0 denote the set of points PΠ such that all the nodes P k to determine the expression S 4 u belong to Π ¯ h , and Π 01 contain the points P, where some of the nodes P k emerge through the side γ j , j=1,2,3,4. We construct the fourth order matching operator S 4 by considering the cases when the point P belongs to one of the sets Π 0 or Π 01 .

Case 1. The point P Π 0 lies on the line connecting two neighboring grid nodes (a grid line).

We place the origin of the rectangular system of coordinates on the node P 0 and direct the positive axis of x along the grid line, so that P=P(δh,0), 0<δ1/2, and take the nodes.

P 0 ( 0 , 0 ) , P 1 ( h , 0 ) , P 2 ( h 2 , 3 h 2 ) , P 3 ( h 2 , 3 h 2 ) , P 4 ( h 2 , 3 h 2 ) , P 5 ( h 2 , 3 h 2 ) .

First, we find the coefficients λ j , j=0,1,2,3, such that the representation

u 0 = λ 0 u+ λ 1 u 1 + λ 2 u 2 + λ 3 u 3
(3.2)

is true for the harmonic polynomials Re z n , n=0,1,2,3, where u=u(P), u k =u( P k ), k=0,1,2,3, z=x+iy. This gives the system

λ 0 + λ 1 + λ 2 + λ 3 = 1 , δ λ 0 + λ 1 + 1 2 λ 2 1 2 λ 3 = 0 , δ 2 λ 0 + λ 1 1 2 λ 2 1 2 λ 3 = 0 , δ 3 λ 0 + λ 1 λ 2 + λ 3 = 0 .
(3.3)

By solving system (3.3) and rearranging (3.2) for u, we obtain the equation

u= u 0 λ 0 λ 1 λ 0 u 1 λ 2 λ 0 u 2 λ 3 λ 0 u 3 .
(3.4)

We now take into consideration the nodes P 4 ( h 2 , 3 h 2 ) and P 5 ( h 2 , 3 h 2 ) which are symmetric to the points P 2 and P 3 , respectively, with respect to the x-axis. Since Im z k =0, k=1,2,3 for y=0, and odd with respect to y, and Re z k , k=0,1,2,3, is even with respect to y, from (3.4) we have

u= u 0 λ 0 λ 1 λ 0 u 1 λ 2 2 λ 0 u 2 λ 3 2 λ 0 u 3 λ 2 2 λ 0 u 4 λ 3 2 λ 0 u 5 .

Hence, we obtain the matching operator S 4 , for which the expression

S 4 u= k = 0 5 λ k u k
(3.5)

gives the exact value of the harmonic polynomial F 3 (x,y) at the point P, where

λ 0 = ( 1 + δ ) ( 1 δ + δ 2 ) , λ 1 = 2 δ + δ 3 3 , λ 2 = λ 4 = ( 1 + δ ) δ 2 , λ 3 = λ 5 = ( 1 + δ ) ( δ + 2 δ 2 ) 6 .

It is easy to check that

λ 0 >0, λ j 0,j=1,2,3, for 0<δ1/2,
(3.6)

and

k = 0 5 λ k =1.
(3.7)

Remark 3.1 When 1/2<δ<1, the node P 1 , which is closest to P, is taken as the origin.

Case 2. The point P Π 0 lies inside a grid cell of the hexagonal grid.

Again, we place the origin of the rectangular system of coordinates at the node P 0 and direct the positive axis of x along the grid line, so that P has the coordinates P(δh, 3 h κ 2 ), where 0<δ,κ1/2. We form an artificial grid by taking the following points:

P 0 ( κ h 2 , 3 h κ 2 ) , P 1 ( h + κ h 2 , 3 h κ 2 ) , P 2 ( h 2 + κ h 2 , 3 h 2 + 3 h κ 2 ) , P 3 ( h 2 + κ h 2 , 3 h 2 + 3 h κ 2 ) , P 4 ( h 2 + κ h 2 , 3 h 2 + 3 h κ 2 ) , P 5 ( h 2 + κ h 2 , 3 h 2 + 3 h κ 2 ) .

Each of the nodes P k , k=0,1,,5, of the artificial grid falls on a grid line, and for the approximation of P the expression

S 4 u= k = 0 5 λ k u ( P k )

is used. As P k , k=0,1,,5, all lie on grid lines, each of these points needs to be approximated using the matching operator as follows:

S 4 u= k = 0 5 λ k S 4 u ( P k ) .

From the distribution of the nodes it becomes obvious that only 17 nodes are needed for this approximation (see Figure 2).

Hence, we form the matching operator as

S 4 u= k = 0 16 ξ k u( P k ),
(3.8)

where ξ k , k=0,,16, are defined by the coefficients obtained earlier and

ξ k 0, k = 0 16 ξ k =1.
(3.9)

For the approximation, it is also important to examine the structure of the hexagonal grid. There are two types of triangles in each hexagon, Type A and Type B, as shown in Figure 1.

Figure 1
figure 1

Shapes of triangles in a hexagon.

We consider triangles of Type A with 0<δ,κ1/2. The nodes used in S 4 u are shown in Figure 2.

Figure 2
figure 2

Type A triangle with 0<δ,κ1/2 .

In the case 1/2<δ<1, 0<κ1/2, the 17 nodes used have the same distribution as the reflection of the nodes in Figure 2 about the line x=0. In the cases 0<δ1/2, 1/2<κ<1 and 1/2<δ,κ<1, the nodes for S 4 are defined analogously.

In the case when P falls into a triangle of Type B, we rotate the artificial grids formed for Type A with an angle of 180, for all four cases of δ and κ specified earlier.

Case 3. P Π 01 , where u= φ j on the side γ j , j=1,2,3,4, and φ j C 4 , λ ( γ j ), 0<λ<1.

We position the origin of the rectangular system of coordinates on γ j so that the point P lies on the positive y axis, and the x axis is in the direction of the vertex A j + 1 along γ j . It is obvious that k = 1 3 b k Im z k =0 if y=0, where z=x+iy. Hence, when the function φ j C 4 , λ ( γ j ), 0<λ<1, is represented using Taylor’s formula about the point x=0 in the neighborhood |z|4h of the origin, we define a k , k=0,1,2,3, of (3.1) as

a k = 1 k ! d k φ j ( 0 ) d x k .

We let

u ˜ (x,y)=u(x,y) k = 0 3 a k Re z k = k = 1 3 b k Im z k +O ( h 4 )

for y>0, and keeping in mind that Im z k is odd extendable, we complete the definition with u ˜ (x,y)= u ˜ (x,y) for y<0. Clearly, in the given neighborhood, u ˜ (x,y) is equal to the harmonic polynomial k = 1 3 b k Im z k , with an accuracy of O( h 4 ). To form an expression for the matching operator S 4 u ˜ , we use

S 4 u ˜ = 0 j 16 μ j ( u k = 0 3 a k Re z k ) ( P j )
(3.10)

or

S 4 u ˜ = 0 j 5 ν j ( u k = 0 3 a k Re z k ) ( P j ),
(3.11)

where

μ j 0, 0 j 16 μ j 1; ν j 0, 0 j 5 ν j 1.
(3.12)

Hence using (3.10) or (3.11), with the addition of the term

( k = 0 3 a k Re z k ) (P),

we have the following representation for the solution u of problem (2.1), (2.2) at any P Π 01

u= S 4 u ˜ + ( k = 0 3 a k Re z k ) (P)+O ( h 4 ) .
(3.13)

Remark 3.2 We obtain the representation (3.13), with a less number of grid nodes P j in (3.10) or (3.11) for the points on the boundary γ of Π.

Let φ= { φ j } j = 1 4 . We express the matching operator S 4 as follows:

S 4 (u,φ)= { S 4 u on  Π 0 , S 4 ( u k = 0 3 a k Re z k ) + ( k = 0 3 a k Re z k ) ( P ) on  Π 01 γ .
(3.14)

Theorem 3.3 Let the boundary functions φ j , j=1,2,3,4, in problem (2.1), (2.2) satisfy the conditions

φ j C 4 , λ ( γ j ),0<λ<1,
(3.15)
φ j ( 2 q ) ( s j )= ( 1 ) q φ j 1 ( 2 q ) ( s j ),q=0,1,2.
(3.16)

Then

max Π ¯ | S 4 (u,φ)u| c 5 h 4 ,
(3.17)

where u is the exact solution of problem (2.1), (2.2).

Proof According to Theorem 3.1 in [15], from conditions (3.15) and (3.16) it follows that u C 4 , λ ( Π ¯ ). Then on the basis of (3.1), (3.5), (3.8), (3.13), and Remark 3.2, we obtain inequality (3.17). □

We define the function u ˆ h as follows:

u ˆ h = S 4 ( u h ,φ)on  Π ¯ ,
(3.18)

where u h is the solution of the finite difference problem (2.6)-(2.8).

Theorem 3.4 Let conditions (2.3) and (2.4) be satisfied. Then the function u ˜ h is continuous on Π ¯ , and

max ( x , y ) Π ¯ | u ˆ h u| c 6 h 4 ,
(3.19)

where u is the solution of problem (2.1), (2.2).

Proof From the construction of the expression S 4 ( u h ,φ) it follows that u ˆ h = u h on Π h , and u ˆ h = φ j on γ j h , j=1,2,3,4. The continuity of u ˆ h on Π follows from the continuity S 4 ( u h ,φ) on each closed triangle Type A and Type B, and from the equality u ˆ h = u h on Π h . By Remark 3.2 and from the condition u ˆ h = φ j on γ j h , j=1,2,3,4, the continuity of the function u ˆ h on the closed rectangle Π ¯ follows. By virtue of (2.3) and (2.4) it follows that u C 6 , λ ( Π ¯ ), 0<λ<1 (see Theorem 3.1 in [15]). Then, on the basis of (3.6), (3.7), (3.9), (3.11) Theorem 2.2, Theorem 3.3 and (3.18), we obtain

max ( x , y ) Π ¯ | u ˆ h u | max ( x , y ) Π ¯ | S 4 ( u , φ ) u | + max ( x , y ) Π ¯ | S 4 ( u h u , 0 ) | c 5 h 4 + k = 0 16 ξ k max ( x , y ) Π ¯ h | u h u | c 6 h 4 .

 □

4 An application of the matching operator in the Block-Grid method

Let G be an open simply connected staircase polygon, let γ j , j=1,2,,N, be its sides, including the ends, and let α j π, α j { 1 2 ,1, 3 2 ,2}, be the interior angle formed by the sides γ j 1 and γ j ( γ 0 = γ N ). Furthermore, let s be the arc length measured along the boundary of G in the positive direction and s j be the value of s at the vertex A j = γ j 1 γ j , ( r j , θ j ) be a polar system of coordinates with pole in A j and the angle θ j taken counterclockwise from the side γ j .

We consider the boundary value problem

Δu=0on G,u= φ j on  γ j ,j=1,2,,N,
(4.1)

where φ j are given functions, and

φ j C 6 , λ ( γ j ),0<λ<1,1jN.
(4.2)

Moreover, at the vertices A j for α j = 1 2 the conjugation conditions

φ j ( 2 q ) ( s j )= ( 1 ) q φ j 1 ( 2 q ) ( s j ),q=0,1,2,3,
(4.3)

are satisfied. At the vertices A j for α j 1 2 no compatibility conditions for boundary functions are required; in particular the values of φ j 1 and φ j at these vertices might be different. Additionally, it is required that when α j 1/2, the boundary functions on γ j 1 and γ j are given as algebraic polynomials of arclength s measured along γ.

Let E={j: α j 1/2,j=1,2,,N}. We call the vertices A j , jE, the singular vertices of the polygon G. We construct two fixed block sectors in the neighborhood of A j , jE, denoted by T j i = T j ( r j i )G, i=1,2, where 0< r j 2 < r j 1 <min{ s j + 1 s j , s j s j 1 }, T j (r)={( r j , θ j ):0< r j <r,0< θ j < α j π}. On the closed sector T ¯ j 1 , jE, we consider the function Q j ( r j , θ j ), which has the following properties:

  1. (i)

    Q j ( r j , θ j ) is harmonic and bounded on the open sector T j 1 ;

  2. (ii)

    continuous everywhere on T ¯ j 1 apart from the point A j , jE when φ j 1 φ j ;

  3. (iii)

    continuously differentiable on T ¯ j 1 A j ;

  4. (iv)

    satisfies the given boundary conditions on γ j 1 T ¯ j 1 and γ j T ¯ j 1 , jE.

The function Q j ( r j , θ j ) with properties (i)-(iv) is given in [17].

Let

R j ( r j , θ j ,η)= 1 α j k = 0 1 ( 1 ) k R ( ( r r j 2 ) 1 / α j , θ α j , ( 1 ) k η α j ) ,jE,
(4.4)

where

R(r,θ,η)= 1 r 2 2 π ( 1 2 r cos ( θ η ) + r 2 )
(4.5)

is the kernel of the Poisson integral for a unit circle.

The approximation of the integral representation given in the following lemma is used to construct an approximate solution of problem (4.1) around the singular vertices A j , jE.

Lemma 4.1 The solution u of problem (2.1), (2.2) can be represented on T ¯ j 2 V j , jE, in the form

u( r j , θ j )= Q j ( r j , θ j )+ 0 α j π ( u ( r j 2 , η ) Q j ( r j 2 , η ) ) R j ( r j , θ j ,η)dη,
(4.6)

where V j is the curvilinear part of the boundary of the sector T j 2 .

Proof The proof follows from Theorems 3.1 and 5.1 in [17]. □

We define the approximate solution in the whole polygon G by applying a version of the Block-Grid method introduced in [8] (see also [9]).

Let us consider, in addition to the sectors T j 1 , T j 2 , the sectors T j 3 and T j 4 , which are also in the neighborhood of each vertex A j , jE, of the polygon G, with 0< r j 4 < r j 3 < r j 2 , r j 3 =( r j 2 + r j 4 )/2 and T k 3 T l 3 =, kl, where k,lE. Furthermore, let G T =G( j E T j 4 ). We give the description of the Block-Grid method on a hexagonal grid:

  1. (i)

    All singular corners A j , jE, are separated by the double sectors T j i = T j ( r j i ), i=2,3, with r j 3 < r j 2 , T k 2 T l 3 =, kl and k,lE. The polygon is covered by overlapping rectangles Π k , k=1,2,,M, and sectors T j 3 , jE, such that the distance from Π k to a singular point A j is greater than r j 4 for all k=1,2,,M and jE.

  2. (ii)

    On each rectangle Π k , the seven point difference scheme for the approximation of Laplace’s equation on a hexagonal grid is used, with step size h k h, h is a parameter, and as an approximate solution on T ¯ j 3 , jE, the harmonic function (4.6) is used.

  3. (iii)

    We use the matching operator S 4 constructed in Section 3 to connect the subsystems.

For obtaining the numerical solution of the algebraic system of equations (2.1), (2.2), we outline the procedure: Let Π k G T , k=1,2,,M, be certain fixed open rectangles with sides a 1 k and a 2 k parallel to the sides of G, and G( k = 1 M Π k )( j E T j 3 )G. We use η k to denote the boundary of the rectangle Π k , V j is the curvilinear part of the boundary of the sector T j 2 and t j =( k = 1 M η k ) T ¯ j 3 .

The overlapping condition is imposed on the arrangement of the rectangles Π k , k=1,2,,M: any point P lying on η k G T , 1kM, or located on V j G, jE, falls inside at least one of the rectangles Π k ( P ) , 1k(P)M, where the distance from P to G T η k ( P ) is not less than some constant ϰ 0 independent of P. The quantity ϰ 0 is called the gluing depth of the rectangles Π k , k=1,2,,M.

We introduce the parameter h(0, ϰ 0 /4] and consider a hexagonal grid on Π k , k=1,2,,M, with maximal possible step h k min{h,min{ a 1 k , a 2 k }/4}. Let Π k h be the set of nodes on Π k , let η k h be the set of nodes on η k , and let Π ¯ k h = Π k h η k h . We denote the set of nodes on the closure of η k G T by η k 0 h , and the set of nodes on Π k h whose distance from the boundary η k G T of Π k is h 2 by η k 0 h . We also have Π k h denoting the set of nodes whose distance from the boundary η k 1 of Π k is h 2 and Π k 0 h = Π k h ( Π k h η k 0 h ). Let t j h be the set of nodes on t j , and let η k 1 h be the set of remaining nodes on η k . We also specify a natural number n[ ln 1 + ϰ h 1 ]+1, where ϰ>0 is a fixed number and the quantities n(j)=max{4,[ α j n]}, β j = α j π/n(j), and θ j m =(m1/2) β j , jE, 1mn(j). On the arc V j we choose the points ( r j 2 , θ j m ), 1mn(j), and denote the set of these points by V j n . Finally, let

ω h , n = ( k = 1 M η k 0 h ) ( k = 1 M η k 0 h ) ( j E V j n ) , G ¯ h , n = ω h , n ( k = 1 M Π ¯ k h ) .

Consider the system of equations

u h =S u h on  Π k 0 h ,
(4.7)
u h = S m u h + E m h ( φ m )on  Π k h , η k 1 h γ m ,
(4.8)
u h = φ m on  η k 1 h γ m ,
(4.9)
u h ( r j , θ j ) = Q j ( r j , θ j ) u h ( r j , θ j ) = + β j k = 1 n ( j ) R j ( r j , θ j , θ j k ) ( u h ( r j 2 , θ j k ) Q j ( r j 2 , θ j k ) ) on  t j h ,
(4.10)
u h = S 4 ( u h ,φ)on  ω h , n ,
(4.11)

where 1k,mM, jE, φ= { φ j } j = 1 N ; S u h , S m u h and E m h ( φ m ) are defined as equations (2.9), (2.10), and (2.11) in Section 2, respectively.

The solution of the system of equations (4.7)-(4.11) is a numerical solution of problem (2.1), (2.2) on G ¯ T (‘nonsingular’ part of the polygon G).

Theorem 4.2 There is a natural number n 0 such that for all n n 0 and h(0, ϰ 0 4 ], where ϰ 0 is the gluing depth, the system of equations (4.7)-(4.11) has a unique solution.

Proof Let v h be a solution of the system of equations

u h = S u h on  Π k 0 h , u h = S m u h on  Π k h , η k 1 h γ m , u h = 0 on  η k 1 h γ m ,
(4.12)
u h ( r j , θ j ) = β j k = 1 n ( j ) R j ( r j , θ j , θ j k ) u h ( r j 2 , θ j k ) on  t j h , u h = S 4 u h on  ω h , n ,
(4.13)

where 1k,mM, jE. To prove the given theorem, we show that max G ¯ h , n | v h |=0. On the basis of the structure of operators S and S j , and the forms (3.5), (3.6), (3.7), (3.8), (3.9), and (3.10)-(3.11) of the matching operator S 4 and by the maximum principle (see Ch. 4, [11]) it follows that the nonzero maximum value of the function v h can be at the points on j E t j h . From estimation (2.29) in [18] the existence of the positive constants n 0 and σ>0 such that for n n 0

max ( r j , θ j ) T j ¯ 3 β j q = 1 n ( j ) R j ( r j , θ j , θ j q ) σ<1
(4.14)

follows. However, taking (4.14) into account in (4.13) we have that the nonzero maximum value can not be at the points on j E t j h either. Since the set G ¯ h , n is connected, from equation (4.12) it follows that max G ¯ h , n | v h |=0. □

Let u h be the solution of the system of equations (4.7)-(4.11). The function

U h ( r j , θ j )= Q j ( r j , θ j )+ β j q = 1 n ( j ) R j ( r j , θ j , θ j q ) ( u h ( r j 2 , θ j q ) Q j ( r j 2 , θ j q ) )
(4.15)

is the approximation of the integral representation (4.6) with the use of the composite mid-point rule. We use the function U h ( r j , θ j ) as an approximate solution of problem (2.1), (2.2) on the closed block T ¯ j 3 , jE (‘singular’ parts of the polygon G).

Let

ϵ h = u h u,
(4.16)

where u h is the solution of system (4.7)-(4.11) and u is the trace of the solution of (2.1), (2.2) on G ¯ h , n . On the basis of (2.1), (2.2), (4.7)-(4.11), and (4.16), ϵ h satisfies the following difference equations:

ϵ h =S ϵ h + r h 1 on  Π k 0 h ,
(4.17)
ϵ h = S m ϵ h + r h 2 on  Π k h , η k 1 h γ m , ϵ h =0on  η k 1 h γ m ,
(4.18)
ϵ h ( r j , θ j )= β j k = 1 n ( j ) R j ( r j , θ j , θ j k ) ϵ h ( r j 2 , θ j k ) + r j h 3 ,( r j , θ j ) t k j h ,
(4.19)
ϵ h = S 4 ϵ h + r h 4 on  ω h , n ,
(4.20)

where 1k,mM, jE and

r h 1 =Suuon  k = 1 M Π k 0 h , r h 2 = S m u+ E m h ( φ m )uon  1 k M Π k h ,
(4.21)
r j h 3 = β j k = 1 n ( j ) R j ( r j , θ j , θ j k ) ( u ( r j 2 , θ j k ) Q j ( r j 2 , θ j k ) ) ( u ( r j , θ j ) Q j ( r j , θ j ) ) on  j E t j h , r h 4 = S 4 ( u , φ ) u on  ω h , n .
(4.22)

Since all the rectangles Π k , k=1,2,,M are located away from the singular vertices A j , jE of the polygon G, at a distance greater than r j 4 >0 independent of h, by virtue of the conditions (2.3) and (2.4), up to sixth order derivatives of the solution of problem (2.1), (2.2) are bounded on k = 1 M Π k . Then, by Taylor’s formula, from (4.21) we obtain

max k = 1 M Π k 0 h | r h 1 | c 1 h 6 , max k = 1 M Π k h | r h 2 | c 2 h 4 .
(4.23)

Furthermore, as ω h , n k = 1 M Π k from (4.22) and Theorem 3.3, we have

max ω h , n | r h 4 | c 3 h 4 .
(4.24)

By analogy to the proof of Lemma 6.2 in [9], it is shown that there exists a natural number n 0 , such that for all nmax{ n 0 ,[ ln 1 + ϰ h 1 ]+1}, ϰ>0 being a fixed number,

max j E | r j h 3 | c 4 h 4 .
(4.25)

Theorem 4.3 There exists a natural number n 0 such that for all nmax{ n 0 ,[ ln 1 + ϰ h 1 ]}, ϰ>0 being a fixed number,

max G ¯ h , n | u h u|c h 4 .
(4.26)

Proof The proof follows from estimations (4.23)-(4.25) and the principle of maximum by analogy to the proof of Theorem 6.1 in [9]. □

Theorem 4.4 Let u h be the solution of the system of equations (4.7)-(4.11), and let an approximate solution of problem (2.1), (2.2) be found on blocks T ¯ j 3 , jE, by (4.15). There is a natural number n 0 such that for all nmax{ n 0 ,[ ln 1 + ϰ h 1 ]}, ϰ>0 being a fixed number, the following estimations hold:

| U h ( r j , θ j )u( r j , θ j )| c 0 h 4 on  T ¯ j 3 ,jE,
(4.27)
| p x p q y q ( U h ( r j , θ j ) u ( r j , θ j ) ) | c p h 4 / r j p 1 / α j on  T ¯ j 3 A j ,jE,
(4.28)

where 0qp, p=0,1, .

Proof Estimation (4.27) is obtained from the integral representation (4.6) and formula (4.15) by using estimations (4.25) and (4.26). Estimation (4.28) for p=0,1, is obtained by using inequality (4.27) and Lemma 6.12 in [17]. □

5 Numerical results and discussion

To support the theoretical results, numerical examples have been solved in two different domains.

Example 5.1 Approximation in a rectangular domain.

Consider the rectangular domain

Π= { ( x , y ) D : 0 < x < 1 , 0 < y < 3 2 } ,

with the boundary γ. The hexagonal grid (2.5), denoted Π h , is assigned to the grid Π, where γ h denotes the set of nodes on the boundary γ.

We consider the problem

Δ u = 0 on  Π h , u = v ( x , y ) on  γ h ,

where

v(x,y)= e y sinx
(5.1)

is the exact solution in the rectangular domain.

This example is solved using the incomplete LU-decomposition method (see [19], Ch. 5), and all the calculations are carried out in double precision. As a convergence test, we request the maximum residual error to be 10−12, and as a starting point v h =0 is used.

Table 1 gives the values obtained in the maximum norm of the difference between the exact and the approximate solutions, for the values of h= 2 k , k=3,4,5,6,7, i.e., ϵ h Π ¯ h = max Π ¯ h |v v h |. The ratios R Π ¯ h m = v v 2 m Π ¯ h v v 2 ( m + 1 ) Π ¯ h have also been included, where O( h 4 ) order of accuracy corresponds to 24 of the value R Π ¯ h m .

Table 1 Approximation in a rectangle with smooth exact solution

Example 5.2 Less smooth function.

We consider the same problem as in Example 5.1 with the exact solution

v(x,y)= 1 2 ln ( x 2 + y 2 ) Re z 7 tan 1 ( y x ) Im z 7 ,
(5.2)

which is less smooth than (5.1). The results obtained are consistent with the theoretical results and are summarized in Table 2.

Table 2 Approximation in a rectangle with less smooth exact solution

Example 5.3 The matching operator.

Examples of the matching operator have also been considered in the domain Π. The coordinate P 1 (0.55,0.4387) is chosen, where P 1 Π 0 . The harmonic function

u(x,y)= e x cosy
(5.3)

is assumed to be the exact solution. The result in Table 3 is obtained using h= 2 4 and demonstrates high accuracy of the above constructed matching operator.

Table 3 Results for approximation of inner points with the matching operator

The second coordinate considered demonstrates the accuracy of the approximation of near-boundary points. The point chosen is P 2 (0.195938,0.02), where P 2 Π 01 , and equation (3.13) is used for approximation. Again, the harmonic function (5.3) is used as the exact solution. Lastly, a point near one of the corners of the domain P 3 (0.005,0.005) has been considered, where the nodes of evaluation emerge outside of the domain from both adjacent sides of the corner. The function

u(x,y)= e y cosx

is used as the exact solution. The results obtained are summarized in Table 4.

Table 4 Results for approximation of near boundary points with the matching operator

Example 5.4 Approximation in an L-shaped domain.

The final example is solved in an L-shaped domain with an angle singularity at the origin, where α 1 π=3π/2. The domain is defined by

Ω= { ( x , y ) : 1 x 1 , 3 2 y 3 2 } Ω 1 ,

where Ω 1 ={(x,y):0x1, 3 2 y0}, and is covered by four overlapping rectangles and a sector. The singular part is defined to be the region

Ω S = { ( x , y ) : 1 2 x 1 2 , 3 4 y 3 4 } Ω 1 s ,

where Ω 1 s ={(x,y):0x 1 2 , 3 4 y0}, and the nonsingular part is Ω N S =Ω/ Ω S . The system of Block-Grid equations is solved by Schwarz’s alternating method. The quadrature nodes on the circular arc, whose radius is taken as 0.75, and the overlapping boundaries of the rectangles are renewed after each Schwarz’s iteration. The nodes on the circular arc, the inner boundaries of the overlapping rectangles and the nodes in the set k = 1 4 η k 0 h are renewed using the matching operator constructed above. Since the boundary functions are harmonic polynomials on the sides γ 1 and γ 0 γ 6 , the nodes whose neighbors emerge outside of the domain from these sides are approximated using the function u Q 1 . Finally, the solution on the singular part is approximated using the integral representation [3, 9].

The problem considered is

Δ u = 0 on  Ω h , u = v ( x , y ) on  γ h ,

where

v(x,y)=θ+ r 2 / 3 sin ( 2 3 θ ) +Re z 5 +Im z 5 ,

is the exact solution. Accordingly, the function Q 1 (x,y) used in the integral representation is constructed as

Q 1 (x,y)=θ+ r 5 ( cos ( 5 θ ) + sin ( 5 θ ) ) .

The results in Tables 5 and 6 show the solution for different pairs (h,N), where N is the number of quadrature nodes, h is the mesh size of the hexagonal grid.

Table 5 The order of convergence in ‘nonsingular’ part
Table 6 The order of convergence in ‘singular’ part

6 Conclusions

The fourth order classical 7-point scheme on a hexagonal grid is applied on a rectangular domain. This leads to some of the nodes emerging from the two parallel sides while approximating points whose distance from the boundary is h 2 . This problem was overcome by devising an approximating equation for near-boundary nodes, which included the use of three inner nodes around the point of evaluation and three points lying on the boundary. The fourth order matching operator has been constructed on the hexagonal grid functions. It is applied to construct a fourth order accurate interpolating function, on the closed rectangle, for the numerical solution of Laplace’s equation on the hexagonal grids. Further, the matching operator and the hexagonal grid approximation in a rectangle are used to obtain and justify the Block-Grid method in solving the Dirichlet problem for Laplace’s equation on staircase polygons.

Numerical examples have been provided as an illustration of the theoretical results mentioned above.

The matching operator constructed can be applied to many other forms of domain decomposition or combined methods. It will also be an interesting study to extend the approximation of the methods to using mixed or Neumann boundary conditions.

References

  1. Volkov EA: Method of composite meshes for finite and infinite domain with a piecewise-smooth boundary. Tr. Mat. Inst. Steklova 1968, 96: 117-148.

    Google Scholar 

  2. Volkov EA: On the method of composite meshes for Laplace’s equation on polygons. Tr. Mat. Inst. Steklova 1976, 140: 68-102.

    Google Scholar 

  3. Dosiyev AA: A fourth order accurate composite grids method for solving Laplace’s boundary value problems with singularities. Comput. Math. Math. Phys. 2002, 42(6):832-849.

    MathSciNet  Google Scholar 

  4. Dosiyev AA, Buranay SC: A fourth order accurate difference-analytical method for solving Laplace’s boundary value problem with singularities. In Mathematical Methods in Engineers. Edited by: Tas K, Machado JAT, Baleanu D. Springer, Berlin; 2007:167-176.

    Chapter  Google Scholar 

  5. Dosiyev AA, Buranay SC:A fourth-order block-grid method for solving Laplace’s equation on a staircase polygon with boundary functions in C k , λ . Abstr. Appl. Anal. 2013., 2013: Article ID 864865 10.1155/2013/864865

    Google Scholar 

  6. Dosiyev AA, Buranay SC, Subasi D: The highly accurate block-grid method in solving Laplace’s equation for nonanalytic boundary condition with corner singularity. Comput. Math. Appl. 2012, 64: 612-632.

    Article  MathSciNet  Google Scholar 

  7. Dosiyev AA: A block-grid method for increasing accuracy in the solution of the Laplace equation on polygons. Russ. Acad. Sci. Dokl. Math. 1992, 45(2):396-399.

    MathSciNet  Google Scholar 

  8. Dosiyev AA: A block-grid method of increased accuracy for solving Dirichlet’s problem for Laplace’s equation on polygons. Comput. Math. Math. Phys. 1994, 34(5):591-604.

    MathSciNet  Google Scholar 

  9. Dosiyev AA: The high accurate block-grid method for solving Laplace’s boundary value problem with singularities. SIAM J. Numer. Anal. 2004, 42(1):153-178. 10.1137/S0036142900382715

    Article  MATH  MathSciNet  Google Scholar 

  10. Volkov EA, Dosiyev AA: A high accurate composite grid method for solving Laplace’s boundary value problems with singularities. Russ. J. Numer. Anal. Math. Model. 2007, 22(3):291-307.

    Article  MATH  MathSciNet  Google Scholar 

  11. Samarskii AA: The Theory of Difference Schemes. Dekker, New York; 2001.

    Book  MATH  Google Scholar 

  12. Ničkovic S, Gavrilov MB, Tošic IA: Geostrophic adjustment on hexagonal grids. Mon. Weather Rev. 2002, 130(3):668-683. 10.1175/1520-0493(2002)130<0668:GAOHG>2.0.CO;2

    Article  Google Scholar 

  13. Majewski D, Liermann D, Prohl P, Ritter B, Buchhold M, Hanisch T, Paul G, Wergen W, Baumgardner J: The operational global icosahedral-hexagonal gridpoint model GME: description and high-resolution tests. Mon. Weather Rev. 2002, 130(2):319-338. 10.1175/1520-0493(2002)130<0319:TOGIHG>2.0.CO;2

    Article  Google Scholar 

  14. Sadourny R, Morel P: A finite-difference approximation of the primitive equations for a hexagonal grid on a plane. Mon. Weather Rev. 1969, 97(6):439-445. 10.1175/1520-0493(1969)097<0439:AFAOTP>2.3.CO;2

    Article  Google Scholar 

  15. Volkov EA: On differentiability properties of solutions of boundary value problems for the Laplace and Poisson equations on a rectangle. Proc. Steklov Inst. Math. 1965, 77: 101-126.

    MATH  Google Scholar 

  16. Kantorovich LV, Krylov VI: Approximate Methods of Higher Analysis. Noordhoff, Groningen; 1958.

    MATH  Google Scholar 

  17. Volkov EA: Block Method for Solving the Laplace Equation and Constructing Conformal Mappings. CRC Press, Boca Raton; 1994.

    MATH  Google Scholar 

  18. Volkov EA: An exponentially converging method for solving Laplace’s equation on polygons. Math. USSR Sb. 1980, 37(3):295-325. 10.1070/SM1980v037n03ABEH001954

    Article  MATH  Google Scholar 

  19. Smith GD: Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford University Press, Oxford; 1985.

    MATH  Google Scholar 

Download references

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Adiguzel A. Dosiyev.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

All authors contributed equally to the writing of this paper. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 2.0 International License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Dosiyev, A.A., Celiker, E. Approximation on the hexagonal grid of the Dirichlet problem for Laplace’s equation. Bound Value Probl 2014, 73 (2014). https://doi.org/10.1186/1687-2770-2014-73

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1687-2770-2014-73

Keywords