Hexagonal grid approximation of the solution of the heat equation on special polygons

*Correspondence: suzan.buranay@emu.edu.tr 1Department of Mathematics, Eastern Mediterranean University, Famagusta, TRNC, Turkey Abstract We consider the first type boundary value problem of the heat equation in two space dimensions on special polygons with interior angles αjπ , j = 1, 2, . . . ,M, where αj ∈ { 2 , 3 , 3 }. To approximate the solution we develop two difference problems on hexagonal grids using two layers with 14 points. It is proved that the given implicit schemes in both difference problems are unconditionally stable. It is also shown that the solutions of the constructed Difference Problem 1 and Difference Problem 2 converge to the exact solution on the grids of order O(h2 + τ 2) and O(h4 + τ ) respectively, where h and √ 3 2 h are the step sizes in space variables x1 and x2 respectively and τ is the step size in time. Furthermore, theoretical results are justified by numerical examples on a rectangle, trapezoid and parallelogram.


Introduction
The use of differential equations in modeling of physical phenomenon is inevitable. Some modeling examples include but are not limited to the modeling of musical instruments such as string and wind instruments using digital waveguides, addressed by Smith [1]. For modeling of infectious diseases see Ahmed et al. [2], and for auto-catalytic chemical reactions see Iqbal et al. [3]. More examples may be given to atmospheric and oceanographic models (see Sadourny et al. [4]-Randall et al. [5]). Analytical techniques cannot solve most of the model problems appearing in practice. Therefore, numerical simulations gain importance in getting an understanding of dynamical systems.
One of the most important issues in numerical methods for the solution of dynamical problems is the well-founded choice of stable and economical computational algorithms. In numerical calculations of nuclear reactors it has been found worthwhile to use the implicit schemes for the solution of a two dimensional age-diffusion equation by Ehrlich and Hurwitz [6]. Most recently, by Ahmed et al. [2] a novel and time efficient positivity preserving numerical scheme to find the solution of an epidemic model involving a reaction-diffusion system in three dimensions has been designed. By Iqbal et al. [3] an unconditionally stable and structure preserving computational technique for fractional order Schnakenberg model has been given.
High accurate implicit schemes on triangular nets whose meshes are equilateral triangles for the two dimensional diffusion equation where, ω > 0 is the thermometric conductivity or the diffusivity were studied by Richtmyer and Morton [7]. The analogue of O(h 2 + τ 2 ) accurate, unconditionally stable three layer scheme with 9 points on hexagonal grids and a three layer scheme with 21 points, a two layer scheme with 14 points both converging with order O(h 4 + τ 2 ) to the exact solution on hexagonal grids were given in Richtmyer and Morton [7]. Therein, the diffusion problem with heat source on a rectangle that irregular grids have centers h 2 units away from the sides of the rectangle at any time moment t with neighboring points emerging through these sides was not considered.
In the last few decades, hexagonal grid methods became of interest in many applied problems such as dynamical meteorology, dynamical oceanography and in the simulation of electrical wave phenomena. Sadourny et al. [4] developed a finite difference scheme for numerical integration of the nondivergent barotropic vorticity equation with an icosahedral-hexagonal grid covering the sphere. The given difference scheme for the advection of vorticity exactly conserved total vorticity, total square vorticity, and total kinetic energy. The authors of [4] showed that the high degree of the isotropy of the hexagonal grid made the formulation of an energy and vorticity conserving scheme simpler than for rectangular grids. Williamson [8] gave finite difference approximations for a nondivergent, barotropic model expressed in term of a streamfunction for an arbitrary triangular grid, which were applied to the spherical geodesic grid. Subsequently, Sadourney and Morel [9] demonstrated the applicability of the hexagonal grids to the primitive equations of fluid dynamics. Also, Sadourney [10], Thacker [11,12], Salmon and Talley [13], Ničkovič [14] considered various aspects of the finite differencing on hexagonal grids. Later Ničkovič et al. [15] showed that hexagonal lattices have some advantages over commonly used square grids. As the authors in [15] mentioned "having better isotropy, they provide more accurate dispersion of gravity waves than square grids do and therefore they can be more appropriate for simulation of smaller-scale divergent processes".
Recently, in the article Lee et al. [16] hexagonal grid finite difference methods were derived in a finite volume approach involving a standard Laplacian. The obtained schemes were used in the simulation of electrical wave phenomena propagated in two dimensional reserved-C type cardiac tissue, exhibiting both linear and spiral waves more efficiently than a similar computation carried out on rectangular finite volume schemes.
Dosiyev and Celiker [17] gave the approximation on the hexagonal grid of the Dirichlet problem for Laplace's equation. The fourth order matching operator on the hexagonal grid was constructed and applied to justify a hexagonal version of the combined Block-Grid method for the Dirichlet problem with corner singularity. Thus, they obtained O(h 4 ) order of accuracy where h is the step size when using the 7-point scheme on the hexagonal grid instead of using 9-point scheme on the rectangular grid giving the computational advantages such as memory space and computational cost. Further, Dosiyev and Celiker [18] investigated a fourth order block-hexagonal grid approximation for the solution of Laplace's equation on special polygons with singularities. It has been justified that in these polygons if the boundary functions away from the singular corners are from the Hölder classes C 4,λ , 0 < λ < 1, the uniform error is of order O(h 4 ) when the hexagon grid is applied in the "nonsingular" part of the domain.
Main contributions by this research are that we give two layer implicit schemes with 14 points by using the hexagonal grids for approximating the solution of first type boundary value problem of the heat equation in two space dimensions on special polygons Ω with interior angles α j π , j = 1, 2, . . . , M, where α j ∈ { 1 2 , 1 3 , 2 3 }. It is assumed that the heat source, the initial and boundary functions are given on x,t (Q T ), 0 < α < 1. Special difference schemes are proposed for the hexagonal grids that have centers h 2 units away from the orthogonal sides of these polygons at time moment t, which have a neighboring point in the pattern emerging through these sides. In Sect. 2, we consider the first type boundary value problem for a two dimensional heat equation on special polygons. In Sect. 3, a two layer implicit difference scheme with 14 points on hexagonal grids is proposed and it is proved that this scheme is unconditionally stable and the solution of the constructed Difference Problem 1 converges to the exact solution on the grids with O(h 2 + τ 2 ) order of accuracy. In Sect. 4, we give a two layer implicit unconditionally stable scheme with 14 points on hexagonal grids and showed that the solution of the constructed Difference Problem 2 converges to the exact solution on the grids with O(h 4 + τ ) order of accuracy. Section 5 is devoted to numerical experiments to justify the obtained theoretical results.

First type boundary value problem of the heat equation on special polygons
Let x = (x 1 , x 2 ), and Ω be an open simply connected polygon and υ j , j = 1, 2, . . . , M, be its sides, including the ends, enumerated counterclockwise (υ 0 = υ M , υ 1 = υ M+1 ), and also let α j π , j = 1, 2, . . . , M, where α j ∈ { 1 2 , 1 3 , 2 3 } be the interior angles formed by the sides υ j-1 and υ j . Furthermore, let S = M j=1 υ j be the boundary of Ω and denote by Ω = Ω ∪ S the closure of Ω. Let Q T = Ω × [0, T), with lateral surface S T , more precisely, the set of points (x, t), x ∈ S and t ∈ [0, T]; also Q T shows the closure of Q T . Let l be a noninteger positive number, x,t (Q T ) be the Banach space of functions u(x, t) that are continuous in Q T together with all derivatives of the form with bounded norm where u (j) and the quantities u α x , u β t for α, β ∈ (0, 1) are defined as We consider the first type boundary value problem for two space dimensional heat equation.

Boundary Value Problem (BVP)
where x = (x 1 , x 2 ) and L ≡ ∂ ∂tω( ∂ 2 and ω > 0 is constant. The differentiability properties of solutions of boundary value problems for the Laplace equation on polygons were given by Volkov [19]. For elliptic equations, the behaviour of solutions near singularities of the boundary of the domain had been treated by Kondrat'ev [20]. For the differentiability properties of solutions of the parabolic equations on cylindrical domains with smooth boundary, see Ladyženskaja et al. [21], and Friedman [22]. The smoothness of solutions of parabolic equations in regions with edges was studied by Azzam and Kreyszig for the Dirichlet [23] and for the mixed boundary value problems in [24]. Hence, in this paper the obtained subsequent theoretical and numerical results are given under the assumption that the heat source function f (x, t) and the initial and boundary functions ϕ(x) and φ(x, t) respectively, are given such that the BVP has unique solution u belonging to the class C Remark 1 It is known that the use of classical finite difference method to solve the boundary value problems with singularities is ineffective. Therefore, a special construction is usually needed for the numerical scheme near the singularities in such a way that the order of convergence is the same as in the case of a smooth solution (see Dosiyev and Celiker [17], Dosiyev [25] and Dosiyev et al. [26]).

Second order accurate difference problem
Let h > 0, we assign Ω h a hexagonal grid on Ω, with step size h, defined as the set of nodes as shown in Fig. 1, Fig. 2 and Fig. 3 when Ω is a rectangle, parallelogram and trapezoid respectively. Let υ h j , j = 1, 2, . . . , M be the set of nodes on the interior of υ j and let υ h We assume that the lengths of the sides of the polygon are given such that irregular hexagon grids only have a right neighboring point or a left neighboring point emerging through the side of Ω when the center of the hexagon is h 2 units away from this side. Accordingly, we shall call these points right ghost points and left ghost points. Further, let Ω * lh , Ω * rh denote the set of interior nodes whose distance from the boundary is h 2 and the hexagon has a left ghost point as shown in Fig. 4 or a right ghost point given in Fig. 5, respectively. We also denote by Ω * h = Ω * lh ∪ Ω * rh and Ω 0h = Ω h \Ω * h . Next, let Let We use the following notations: P 0 denotes the center of the hexagon. Patt(P 0 ) is the pattern of the hexagon consisting the neighboring points P i , i = 1, . . . , 6. Inc denotes the incidence matrix related to the neighborhood topology of all the hexagon centers. λ s (K) present the eigenvalues λ s , s = 1, 2, . . . , N of real matrix K ∈ R N×N .
Also u k+1 ∈ Ω * rh γ τ as also given in (29). Analogously, the values u k P i i = 0, . . . , 6 and u k P A present the exact solution at the same space coordinates of P i i = 0, . . . , 6 and P A respectively, but at time level t. Also we use the notations u k+1 h,τ ,P i , i = 0, . . . , 6, u k+1 h,τ ,P A , and u k h,τ ,P i , i = 0, . . . , 6, u k h,τ ,P A to present the numerical solution at the same space coordinates of P i , i = 0, . . . , 6 and P A for time moments t + τ and t = kτ , respectively. Further, f where p is as defined in (29).
For the numerical solution of the BVP we propose the following difference problem.

Theorem 2 The order of approximation of the implicit scheme in Difference Problem
Proof Let (x 1 , x 2 , t + τ ) ∈ Ω * γ τ be the center of the hexagon (P 0 ) at time moment t + τ . By the continuity of the solution u, the heat equation (10) is also satisfied at the boundary points ( p, x 2 , t + τ 2 ) denoted by P A . Therefore, we give the O(h 2 + τ 2 ) difference approximation of the heat equation (10) Using Eq. (30) the following equation can be derived for the sum of left ghost points P 2 at time moments t + τ , and t = kτ , Analogously, we use the following difference approximation of the heat equation (10) at the boundary points ( p, h,τ ,P 5 with the order of approximation O(h 2 + τ 2 ). The sum of right ghost points P 5 at time moments t + τ , and t = kτ is obtained from (32), Using (31), (33) and (18) we obtain the scheme (19). Let the error function be ε h,τ = u h,τu.
Then ε h,τ satisfies the following difference problem: where, Let Θ h,τ be the operator that coincides with Θ h,τ for the points in Ω 0h γ τ , and coincides with Θ * h,τ for the points in Ω * h γ τ . Analogously, let Λ h,τ be the operator that coincides with Λ h,τ for the points in Ω 0h γ τ and coincides with Λ * h,τ for the points in Ω * h γ τ . Further Ψ k denotes the truncation error Ψ k 1 and Ψ k 2 for the points belonging to Ω 0h γ τ and Ω * h γ τ respectively. Then the system (34)-(37) can be given as Using Taylor's expansion around the point (x 1 , x 2 , t + τ 2 ) and from the assumption that order of the approximation. Thus, the order of the approximation of the implicit scheme in Difference Problem 1 is O(h 2 + τ 2 ).
Next we analyze the stability for Difference Problem 1 by using spectral method. Let us label the interior grid points in Ω h γ τ by Q j , j = 1, 2, . . . , N at every time level along spatial variable x 1 (lexicographical ordering). The neighboring topology of all hexagon centers can be given by the set and shows the sparsity pattern of the incidence matrix Inc. Thus, the entries of the matrix Inc ∈ R N×N are The algebraic linear system of equations obtained by the Difference Problem 1 can be given in matrix form: where K, S ∈ R N×N are given as and Here U k+1 , U k , F k+ 1 2 and G k * are vectors of order N and F k+ 1 2 and G k * are obtained by evaluating the heat source function in (22), (23) at time level (k + 1 2 )τ and the boundary and initial functions in Difference Problem 1 (18)-(21), respectively. Also k * denotes that values at time moments t + τ , and t = kτ are used, and I is the identity matrix. D 1 is a diagonal matrix with entries 7 3 if Q j ∈ Ω * h γ τ , j = 1, 2, . . . , N. Proof (a) Using (43) if Q i ∈ Patt(Q j ) for i = j, 1 ≤ i, j ≤ N this implies that Q j ∈ Patt(Q i ) giving Inc T = Inc. Thus, B is symmetric and since a hexagonal grid is a connected grid in a simply connected polygon Ω, using (47) one can easily show that the matrix B is an irreducibly diagonally dominant matrix with b ii > 0, i = 1, . . . , N . Hence by Theorem 4, B is a positive definite matrix.
(b) The main diagonal entries k ii > 0, i = 1, . . . , N , and k ij ≤ 0 for i = j, and 1 ≤ i, j ≤ N . Since k ii > N j=1,j =i |k ij | for all i = 1, 2, . . . , N , the matrix K = I + ωτ h 2 D 1 -ωτ 3h 2 Inc is strictly diagonally dominant matrix and Lemma 3 implies that K is an M-matrix and its inverse K -1 ≥ 0. Also K T = (I + ωτ h 2 B) T = K and K is symmetric real matrix. Therefore, all eigenvalues λ s of K are real. Using Theorem 4 we find that Re(λ s (K)) = λ s (K) > 0 for all eigenvalues λ s of K ; thus, K is a symmetric positive definite matrix.

Theorem 6
The implicit scheme of the Difference Problem 1 is unconditionally stable and the solution u h,τ converges to the exact solution u with order of accuracy (h 2 + τ 2 ).
Proof On the basis of Lemma 5, the matrices B and K are symmetric and positive definite matrices. Since K is symmetric, K -1 is also symmetric and the eigenvalues of K -1 satisfy 0 < λ s (K -1 ) ≤ 1 1+ ωτ h 2 min 1≤s≤N (λ s (B)) < 1 for ωτ h 2 > 0 and we get K -1 2 < 1. Also where det(I + ωτ h 2 B) and Adj(I + ωτ h 2 B) are the determinant and the adjoint matrix of K = I + ωτ h 2 B, respectively. Thus K -1 S is a real symmetric matrix; therefore, the eigenvalues λ s (K -1 S) are real. Using the Gershgorin theorem for the estimation of the spectrum of the matrix B gives Since B is a symmetric positive definite matrix, it follows that B = P T DP with P orthogonal and D diagonal matrix of eigenvalues λ s (B). Then K = I + ωτ h 2 B = P T (I + ωτ h 2 D)P and K -1 = P T (I + ωτ h 2 D) -1 P and that is, the matrix K -1 S is similar to (I + ωτ h 2 D) -1 (I -ωτ h 2 D) so from (49) we get Using (52) and by induction follows that Since K -1 S is a real symmetric matrix it is also a normal matrix; hence, the von Neumann condition for stability is sufficient as well as necessary for stability (see Lax and Richtmyer [28]). Therefore, Eq. (53) shows that the implicit scheme in Difference Problem 1 is unconditionally stable. The error function ε h,τ satisfying (40)-(42) can be given in the matrix form (45) at time level t = (k + 1)τ as where, E is vector of order N . Thus, from Theorem 2 and (53), (54) we have where, c 1 is positive constant independent from h and τ and depends on the bounded derivatives of the solution u of the form (2) of at most sixth order in the truncation error Ψ k as given in (38) and (39). Let ε k+1 h,τ C = Ω h γ τ {t=(k+1)τ } max|ε k+1 h,τ | = E k+1 ∞ , then on the basis of norm concordance and using (55) we get Therefore, the solution u h,τ converges to the exact solution u with order of accuracy (h 2 + τ 2 ).

Fourth order accurate implicit difference problem
| (x 1 ,x 2 ,t+τ ) . We give the following difference problem for the solution of the given BVP.

Difference Problem 2
and p, p, η are as defined in (29). We remark that ψ in (61) can be taken as (see also Lee et al. [16])

Theorem 7 The order of approximation of the implicit scheme in Difference Problem 2 is O(h 4 + τ ).
Proof Let (x 1 , x 2 , t + τ ) ∈ Ω * lh γ τ be the center of the hexagon (P 0 ) at time level t + τ . The heat equation (10) is also satisfied at the boundary points ( p, x 2 , t + τ ) and ( p, x 2 , t) where, p = x 1 -h 2 . We give the difference approximation of the heat equation (10) at these boundary points (P A ), h,τ ,P 2 and from (70) results Analogously, when (x 1 , x 2 , t + τ ) ∈ Ω * rh γ τ for the right ghost points approximation we obtain Using (71)-(74) and (57) we obtain the scheme (58). Let the error function be ε h,τ = u h,τ -u. Then ε h,τ satisfies the following difference problem: where, Let Θ h,τ be the operator that coincides with Θ h,τ for the points in Ω 0h γ τ , and coincides with Θ * h,τ for the points in Ω * h γ τ . Analogously, let Λ h,τ be the operator that coincides with Λ h,τ for the points in Ω 0h γ τ and coincides with Λ * h,τ for the points in Ω * h γ τ . Also Ψ k denotes the truncation error Ψ k 1 and Ψ k 2 for the points belonging to Ω 0h γ τ and Ω * h γ τ , respectively. Then the system (75)-(78) can be given as Using Taylor's expansion at the point (x 1 , x 2 , t + τ ) and from the assumption that u ∈ C 6+α,3+ α where K, S ∈ R N×N are given as F k * and G k * are vectors of order N obtained by evaluating the heat source function f in (61), (62) and the boundary and initial function values in Difference Problem 2 (57)-(60), respectively. Also k * denotes that values from (k + 1)τ and kτ are used and D 1 , D 2 are diagonal matrices with entries 17 24 if Q j ∈ Ω * h γ τ , 14 3 if Q j ∈ Ω * h γ τ , respectively.

Lemma 8
The matrices K, S in (85) and B in (86) are symmetric positive definite matrices.
Proof The matrix S is nonnegative and strictly diagonally dominant matrix and B is irreducibly diagonally dominant matrices with positive main diagonal entries. Since Inc is symmetric we see that B, S are also symmetric. Thus, Theorem 4 implies that B, S have positive eigenvalues. Therefore B, S are symmetric positive definite matrices. Using K = S + ωτ h 2 B and that the sum of two symmetric positive definite matrices is symmetric positive definite shows that K is symmetric positive definite.
On the basis of Lemma 8, S is invertible and the algebraic linear system (84) can be rewritten as

Lemma 9
The matrix A in (90) is symmetric positive definite.
Proof S = I -1 16 B and B have the same basis vectors spanning their eigenspaces, therefore on the basis of Lemma 8 and using that every real symmetric matrix is orthogonal equivalent to a real diagonal matrix we have S = P T 1 P and B = P T 2 P where 1 and 2 are diagonal matrices and P is orthogonal matrix of which the columns are the normalized basis vectors spanning the eigenspaces of S and B. From -1 So S -1 B is symmetric, thus S -1 B commutes. Since the product of two symmetric positive definite matrices that commute is also symmetric positive definite (see [27] and [29]) we have λ s ( S -1 B) > 0 as that is, S -1 B is similar to the symmetric matrix B We have Proof Using (90) Lemma 9 implies that for ωτ h 2 > 0 and since A is symmetric matrix A -1 is also symmetric matrix, therefore A -1 2 < 1. Also using (85)-(90) and the Gershgorin theorem to estimate the spectrum of B we get 0 < λ s ( B) ≤ 8, and ( K) -1 2 = ( S A) -1 2 ≤ 2. Multiplying both sides of (89) by A -1 and taking the second norm and by using norm properties and induction give Hence, the inequality (94) gives the Von Neumann necessary condition for stability. Since A is real symmetric matrix it is also normal matrix and the condition (94) is as well as the sufficient condition for stability (see Lax and Richtmyer [28]). Therefore, equation (95) yields that the implicit scheme (18), (19) is unconditionally stable. The error function ε h,τ satisfying (81)-(83) can be given in the matrix form (89) at time level t = (k + 1)τ as where, E is vector of order N . Thus, from Theorem 2 and (95), (96) we have where, c 2 is positive constant independent from h and τ and depends on the bounded derivatives of the solution u of the form (2) of at most sixth order in the truncation error Ψ k as given in (79) and (80). Using norm concordance and the inequality (97) we get Hence, the solution u h,τ converges to the exact solution u on the hexagonal grids with order of accuracy (h 4 + τ ).

Numerical results
We consider the open polygon Ω as the rectangle Ω Rec = {(x 1 , x 2 ) : 0 < x 1 < 1, 0 < x 2 < ) is taken one. All the computations are performed using Mathematica in double precision on a personal computer with properties AMD Ryzen 7 1800X Eight Core Processor 3.60 GHz. Also we used conjugate gradient method to solve the obtained algebraic linear system of equations at each time level. All tables adopt the following notations: CT Exi , i = 1, 2, 3 present the total Central Processing Unit time in seconds per time level for the Example 1, Example 2 and Example 3, respectively. neg means that CT Exi , i = 1, 2, 3, is less than one milliseconds.
2 Exi 2 Ex1 for the considered domains respectively. Table 1, Table 2, Table 3 demonstrate the CT Ex1 , CT Ex2 and the maximum norm of the errors for h = 2 -μ , μ = 3, 4, 5, 6, 7, 8 when τ = 2 -λ , λ = 8, 9, 10, 11, 12, 13 and the order of convergence of u h,τ to the exact solution u with respect to h and τ obtained by using the constructed Difference Problem 1 for the Example 1 and Example 2 on a rectangle, trapezoid and parallelogram, respectively. These tables show that the proposed Difference Problem 1 has quadratic convergence order in both the spatial and time variables. Next we solve both examples by using the Difference Problem 2 and denote the obtained order of convergence of the approximate solution u h,τ to the exact solution u for the Ex-   4 Exi 4 Exi for the considered domains respectively. Table 4, Table 5, Table 6 show the CT Ex1 , CT Ex2 , maximum norm of the errors for h = 2 -μ , μ = 4, 5, 6, 7, 8 when τ = 2 -λ , λ = 6, 10, 14, 18, 22 and the order of convergence of u h,τ to the exact solution u with respect to h and τ obtained by using the constructed Difference Problem 2 for the Example 1 and Example 2 on rectangle, trapezoid and parallelogram respectively. These tables demonstrate that the approximate solution u h,τ of the proposed Difference Problem 2 converges to the exact solution u with fourth order in the spatial variables and linearly with respect to time variable t.       The exact solution of Example 3 is not given. Using the proposed Difference Problem 1 we obtain the approximate solution u 2 -μ, 2 -λ (x 1 , x 2 , t) at each time level for μ = 4, 5, 6, and λ = 9, 10, 11, respectively. Table 7 presents u 2 -μ, 2 -λ (x 1 , x 2 , t) at the grid points (0.125, Rec (P) = log 2 u 2 -4 ,2 -9 (P)u 2 -5 ,2 -10 (P) u 2 -5 ,2 -10 (P)u 2 -6 ,2 -11 (P) .
are shown in Table 8. By analyzing the values of (105) and (106) in the fifth columns of Table 7 and Table 8 respectively, we conclude that the order of convergence is quadratic in both the spatial and time variables on t = 1 when Difference Problem 1 is used, and it is of  Table 8 Solution at some points on t = 1, and the order of convergence by using Difference Problem 2 for the Example 3