Two Stage Implicit Method on Hexagonal Grids for Approximating the First Derivatives of the Solution to the Heat Equation

: The ﬁrst type of boundary value problem for the heat equation on a rectangle is considered. We propose a two stage implicit method for the approximation of the ﬁrst order derivatives of the solution with respect to the spatial variables. To approximate the solution at the ﬁrst stage, the unconditionally stable two layer implicit method on hexagonal grids given by Buranay and Arshad in 2020 is used which converges with O (cid:0) h 2 + τ 2 (cid:1) of accuracy on the grids. Here, h and √ 32 h are the step sizes in space variables x 1 and x 2 , respectively and τ is the step size in time. At the second stage, we propose special difference boundary value problems on hexagonal grids for the approximation of ﬁrst derivatives with respect to spatial variables of which the boundary conditions are deﬁned by using the obtained solution from the ﬁrst stage. It is proved that the given schemes in the difference problems are unconditionally stable. Further, for r = ωτ h 2 ≤ 37 , uniform convergence of the solution of the constructed special difference boundary value problems to the corresponding exact derivatives on hexagonal grids with order O (cid:0) h 2 + τ 2 (cid:1) is shown. Finally, the method is applied on a test problem and the numerical results are presented through tables and ﬁgures.


Introduction
Numerical methods have gained considerable attention in many applications, since the exact solution of many problems arising in the models of chemistry, physics, biology, engineering, and many other fields of different sciences is an uphill task. Modeling of these problems leads us to consider a number of physical quantities, representing physical phenomena on a modeling domain. These physical quantities then occur in the model via functions or function derivatives of which for a considerable number of them the Newtonian concept of a derivative satisfies the complexity of the natural occurrences. However, "time's evolution and changes occurring in some systems do not happen in the same manner after a fixed or constant interval of time and do not follow the same routine as one would expect. For instance, a huge variation can occur in a fraction of a second, causing a major change that may affect the whole system's state forever" as stated in [1]. Indeed, it has turned out recently that many of phenomena involved in many branches of chemistry, engineering, biology, ecology, and numerous domains of applied sciences can be described very successfully by models using fractional order differential equations such as acoustic dissipation, viscoelastic systems, mathematical epidemiology, continuous time random walk, and biomedical engineering (see [1] and references therein). Analytical techniques can not solve most of these models with the conventional integer-order derivative, and models with fractional derivatives appearing in practice. Hence, various methods for the solution of these model problems have been developed and proposed in numerous works, in order to provide an improved description of the phenomenon under investigation.
Common numerical methods include finite difference method, finite element method, finite volume method, variational iteration method, adomian or homotopy analysis, wavelet method, etc.
Most recently, for the model problems of chaotic attractors, exhaustive studies were given for the solvability of chaotic fractional systems with 3D four-scroll attractors in [2], for the Proto-Lorenz system in its chaotic fractional and fractal structure in [3], and for a new auto-replication in systems of attractors with two and three merged basins of attraction via control in [4]. To model some symbiosis systems describing commmensalism and predator-prey processes, the Atangana-Baleanu derivative operator was applied in [5] and a numerical approximation technique was given. Model problems involving a derivative with fractional parameter β and the application to transport-convection differential equations were given in [1]. In addition, the use of a control parameter to control on processes related to stationary state system, and on relaxation and diffusion, was studied in [6]. Further, in [7] a comparative analysis between differential fractional operators including the Atangana-Baleanu derivative and the Caputo-Fabrizio derivative applied to the non-linear Kaup-Kupershmidt equation was given and methods of performing numerical approximations of the solutions were presented. Furthermore, for the numerical solution of fractional Volterra type model problems, recent studies include [8] that a class of system of nonlinear singular fractional Volterra integro-differential equations was solved by a proposed computational method. In addition, [9], in which delay-dependent stability switches in fractional differential equations were studied and obtained results were illustrated via a fractional Lotka-Volterra population model. Moreover, [10] as a biological fractional n-species delayed cooperation model of Lotka-Volterra type was presented. Examples to recent studies on numerical solutions of model problems in fractional structure with both stiff and nonstiff components and the leading-edge model problem can be given to [11,12], respectively. A second-order diagonally-implicit-explicit multi-stage integration method was given in [11] for the solution of problems with both stiff and nonstiff components. An implicit method for numerical solution of singular and stiff initial value problem was developed in [12]. For the epidemic models latest studies include [13] that the Crank Nicolson difference scheme and iteration method were used for finding the approximate solution of system of nonlinear observing epidemic model. In addition, [14], in which a novel and time efficient positivity preserving numerical scheme was designed to find the solution of epidemic model involving a reaction-diffusion system in three dimension.
Apart from rectangular grids, hexagonal grids have been also used to develop finite difference methods for the approximate solution of modeled problems in many applied sciences for more than the half century. These studies include the hexagonal grid methods given in meteorological and oceanographic applications by [15][16][17][18][19][20][21][22][23][24][25], of which favorable results were obtained compared with rectangular grids. Hexagonal grids were applied in reservoir simulation in [26] and it was shown that for seven-point floods, hexagonal grid method provides good numerical accuracy at substantially less computational work than rectangular grid method (five or nine point methods). Hexagonal grids were also used in the simulation of electrical wave phenomena propagated in two dimensional reserved-C type cardiac tissue in [27]. The exhibited linear and spiral waves were more efficient than similar computation carried out on rectangular finite volume schemes. Furthermore, hexagonal grids were applied to approximate the solution of the first type boundary value problem of the heat equation in [28][29][30], convection-diffusion equation in [31], and Dirichlet type boundary value problem of the two dimensional Laplace equation in [32]. In the recent study [29], the solution of first type boundary value problem of heat equation on special polygons with interior angles α j π, j = 1, 2, ..., M, for α j ∈ 1 2 , 1 3 , 2 3 where, ω > 0 and f is the heat source by using hexagonal grids has been given. Therein, two implicit methods named as Difference Problem 1 and Difference Problem 2 both on two layers with 14-points have been proposed. It was assumed that the heat source and the initial and boundary functions are given such that the exact solution belongs to the Hölder space C 6+α,3+ α 2 x,t , 0 < α < 1. Under this condition, it was proved that the given Difference Problem 1 and Difference Problem 2 converge to the exact solution on the grids with O h 2 + τ 2 and O h 4 + τ order of accuracy, respectively.
On the other hand, as well as the solution of the modeled problem, the first order partial derivatives of the solution are also essential to determine the rate of change in the solution and the gradients which determines important phenomena in that model. Such as in the electrostatics the first derivatives of electrostatic potential function define electric field. As the calculation of ray tracing in electrostatic fields by the interpolation methods require the specification at each mesh point not only the potential function Φ but also the gradients ∂Φ ∂x 1 , ∂Φ ∂x 2 and the mixed derivative ∂ 2 Φ ∂x 1 ∂x 2 . Motivated by this aim, in this study a second order accurate two stage implicit method for the approximation of the first order derivatives of the solution u(x 1 , x 2 , t) of (1) with respect to the spatial variables x 1 and x 2 on rectangle D is developed.
The research is organized as follows: In Section 2, we consider the first type boundary value problem for the heat equation in (1) on a rectangle D under the assumption that the heat source and the initial and boundary functions are given such that on T], and D is the closure of D. In addition, hexagonal grid structure and basic notations are given. Further, at the first stage, a two layer implicit method on hexagonal grids given in [29] with O h 2 + τ 2 order of accuracy, where h and √ 3 2 h are the step sizes in space variables x 1 and x 2 , respectively, and τ is the step size in time used to approximate the solution u(x 1 , x 2 , t). For the error function when r ≤ 3 7 , we provide a pointwise prior estimation depending on ρ(x 1 , x 2 , t), which is the distance from the current grid point to the surface of Q T . In Sections 3 and 4, the second stages of the two stage implicit method for the approximation to the first order derivatives of the solution u(x 1 , x 2 , t) with respect to the spatial variables x 1 and x 2 are proposed, respectively. It is proved that the constructed implicit schemes at the second stage are unconditionally stable (see Theorem 1 in [33] which gives the sufficient condition of stability). For r = ωτ h 2 ≤ 3 7 , priory error estimations in maximum norm between the exact derivatives ∂u ∂x 1 , ∂u ∂x 2 and the obtained corresponding approximate solutions are provided giving O h 2 + τ 2 order of accuracy on the hexagonal grids. In Section 5, a numerical example is constructed to support the theoretical results. We applied incomplete block preconditioning given in [34] (see also [35,36]) for the conjugate gradient method to solve the obtained algebraic systems of linear equations for various values of r. In Section 6, conclusions and some remarks are given. , with lateral surface S T more precisely the set of points (x, t), x = (x 1 , x 2 ) ∈ S and t ∈ [0, T] also Q T is the closure of Q T . We consider the first type boundary value problem for the two space dimensional heat equation:

BVP(1):
where ω is a positive constant. We assume that the heat source function f (x 1 , x 2 , t) and the initial and boundary functions ϕ(x 1 , x 2 ) and φ(x 1 , x 2 , t), respectively, are given such that the Problems (2)-(4) has a unique solution u belonging to the Hölder class C 7+α, 7+α 2 x,t Q T . For the smoothness of solutions of parabolic equations in regions with edges, see [37] for the Dirichlet and [38] for the mixed boundary value problems. Let h > 0, with h = a 1 /N 1 , where N 1 is positive integer and assign D h a hexagonal grid on D, with step size h, defined as the set of nodes Let γ h j , j = 1, ..., 4 be the set of nodes on the interior of γ j and let γ h Further, let D * lh , D * 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 Figure 1 or a right ghost point as presented in Figure 2, emerging through the left or right side of the rectangle, respectively. We also denote by D * h = D * lh ∪ D * rh and D 0h = D h \D * h .  Next, let and the set of internal nodes and lateral surface nodes be defined by respectively. Let D * lh γ τ = D * lh × γ τ ⊂ D h γ τ and D * rh γ τ = D * rh × γ τ ⊂ D h γ τ and D * h γ τ = D * lh γ τ ∪ D * rh γ τ . In addition, D 0h γ τ = D h γ τ \D * h γ τ and D h γ τ is the closure of D h γ τ . Let P 0 denote the center of the hexagon and Patt(P 0 ) denote the pattern of the hexagon consisting the neighboring points P i , i = 1, ..., 6. In addition, u k+1 P i denotes the exact solution at the point P i and u k+1 P A denotes the value at the boundary point for the time moment t + τ as follows: where the value of p = 0 if P 0 ∈ D * lh γ τ and p = a 1 if P 0 ∈ D * rh γ τ as also given in (21). 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 = kτ. Further, 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 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 and f For the numerical solution of the BVP(1) the following difference problem (named as Difference Problem 1) was given in [29] which we will consider as the Stage 1 H 2nd of the two stage implicit method: ψ and By numbering the interior grid points using standard ordering as L j , j = 1, 2, ..., N, the obtained algebraic linear system of equations in matrix form given in [29] is as follows: where A, B ∈ R N×N are and and u k , q k u ∈ R N . The matrix Inc is the incidence matrix of the neighboring topology with entries unity for the points in the pattern of the hexagon center. In addition, I is the identity matrix, D 1 is a diagonal matrix with entries (22) is a nonsingular M-matrix and is also a symmetric positive definite matrix.

Lemma 1. (a) The matrix A in
Proof. Proof is given in ( [29]).

Pointwise Priory Estimation For the Error Function (27)-(30)
Consider the following systems for k = 0, ..., M − 1, where g 1 , g 2 and g 1 , g 2 are given functions. The algebraic systems (33)- (36) and (37)-(40) can be written in matrix form respectively, for every time level k = 0, ..., M − 1 where A and B are the matrices given in (22) and q k , q k , g k , g k ∈ R N . We also use the partial order K 1 ≤ K 2 which means that

Theorem 1.
For the solution of the problem (27)-(30), the following inequality holds true where and u is the exact solution of BVP(1) and ρ(x 1 , x 2 , t) is the distance from the current grid point in D h γ τ to the surface F of Q T .
Proof. We consider the system and the majorant functions in which ε u i , i = 1, 2, 3 satisfy the difference boundary value problem Therefore, we write the difference problems (56)-(59) and (63)-(66) for fixed k ≥ 0 in matrix form respectively, where A and B are the matrices given in (22) and e u,k i , ε u,k i , i = 1, 2, 3 and ε u,k , e u,k ∈ R N . Using (52)  on D 0h γ τ , and 5 6

Difference Problem Approximating ∂u ∂x 1 on Hexagonal grids with O(h 2 + τ 2 ) Order of Accuracy
We use the notations ∂ x 1 f and ∂ x 1 f . Given the boundary value Problems (2)

Lemma 3.
The following inequality holds true for r = ωτ h 2 ≤ 3 7 , where u is the solution of the boundary value Problems (2)-(4) and u h,τ is the solution of the difference problem in Stage 1 H 2nd and Ω 1 (h, τ) and d are as given in (52) and (55), respectively.

Lemma 4. The following inequality
holds true for r = ωτ h 2 ≤ 3 7 where u h,τ is the solution of the difference problem in Stage 1 H 2nd and Ω 1 and d are as given in (52) (73) and (74) give the second order approximation of ∂u ∂x 1 , respectively. From the truncation error formula (see [39]) it follows that Analogously, Using Lemma 3 and the estimations (81) and (82) follows (80).

Theorem 2.
The implicit scheme given in Stage 2( ∂u ∂x 1 ) is unconditionally stable.
Proof. The obtained algebraic linear system of Equations (83)-(87) can be given in matrix form: k = 0, 1, ..., M − 1, where A and B are the matrices given in (22) and v k , q k v ∈ R N . On the basis of the assumption that the exact solution v of the BVP(2) belongs to C 6+α,3+ α 2 x,t Q T and using Lemma 1 and by induction we get v k+1 Thus, Lax and Richtmyer sufficient condition for stability given in Theorem 1 of [33] is satisfied and the scheme is unconditionally stable.

Second Order Approximation of ∂u ∂x 2 on Hexagonal grids
Let the BVP(1) be given, then we use the notations ∂ x 2 f and and denote q i = ∂u ∂x 2 on S T γ i , i = 1, 2, ..., 5 and setup the next boundary value problem for z = ∂u

BVP(3):
where L is the operator in (72) and f (x 1 , x 2 , t) is the given function in (2). We assume that x,t Q T and take and ϕ(x 1 , x 2 ) given in (3), φ(x 1 , x 2 , t) given in (4) are the initial and boundary functions, respectively, u h,τ is the solution of the difference problem in Stage 1 H 2nd .

Lemma 5.
The following inequality holds for r = ωτ h 2 ≤ 3 7 , where u is the solution of the boundary value Problems (2)-(4) and u h,τ is the solution of the difference Problems (10)- (13) in Stage 1 H 2nd and Ω 1 (h, τ) and d are as given in (52) and (55), respectively.

Lemma 6. The following inequality is true
and u h,τ is the solution of the difference problem in Stage 1 H 2nd and Ω 1 (h, τ) and d are as given in (52) and (55), respectively.
Proof. The obtained algebraic linear system of Equations (130)-(134) can be given in matrix form: for k = 0, 1, ..., M − 1, where, A, B are as given in (22) and z k , q k z ∈ R N . On the basis of the assumption that the exact solution z of the BVP(3) belongs to C 6+α,3+ α 2 x,t Q T and using Lemma 1 and induction we get Therefore, the scheme is unconditionally stable.

Theorem 5.
The solution z h,τ of the finite difference problem given in Stage 2 H 2nd ∂u for r = ωτ h 2 ≤ 3 7 , where κ, δ are as given in (146), (147) respectively and z = ∂u ∂x 2 is the exact solution of BVP(3) .

Numerical Results
A test problem is constructed of which the exact solution is known to show the efficiency of the proposed two stage implicit method. The rectangle D is taken as D = (x 1 , , and t ∈ [0, 1]. All the computations are performed using Mathematica in double precision on a personal computer with properties AMD Ryzen 7 1800X Eight Core Processor 3.60GHz. To solve the obtained linear algebraic system of equations, we applied incomplete block-matrix factorization of the block tridiagonal stiffness matrices which are symmetric M-matrices for the all considered pairs of (h, τ) using two-step iterative method for matrix inversion. Then these incomplete block-matrix factorizations are used as preconditioners for the conjugate gradient method as given in [34] (see also [35,36] ). We use the following notations in tables and figures: Analogously, the order of convergence of the approximate solution z 2 −µ ,2 −λ to the exact solution z = ∂u ∂x 2 obtained by using the two stage implicit method H 2nd ∂u ∂x 2 is given by We remark that the values of (170), (171) are ≈ 2 2 showing that the order of convergence of the approximate solution v 2 −µ ,2 −λ to the exact solution v = ∂u ∂x 1 and the order of convergence of the approximate solution z 2 −µ ,2 −λ to the exact solution z = ∂u ∂x 2 are second order both in the spatial variables x 1 , x 2 and in time t, accordingly.
Test Problem:    TCT H 2nd ∂u ∂x 2 , maximum norm of the errors for the same pairs of (h, τ) as in Table 3 and the order of convergence of z h,τ to the exact derivative z = ∂u ∂x 2 with respect to h and τ are used on S h T γ i , i = 1, 3 for the Stage 2 H 2nd ∂u ∂x 1 .   . Table 8 gives CT with second order both in the spatial variables x 1 , x 2 and the time variable t with better error ratios.
is applied on a test problem and the obtained theoretical order of convergence is shown numerically using tables and figures. Remark 1. The methodology given in this research may be used to construct highly accurate implicit splitting schemes (fractional step methods) and alternating direction methods (ADI) (see [40][41][42][43]) for the approximation of the first order partial derivatives of solution of first type boundary value problem of heat equation in three space dimension. Additionally, this approach may be extended to find the spatial derivatives of the solution of the time-fractional structure of the heat equation. Since as well as its solution, the computation of the derivatives of the solution of the fractional model problem are essential. Such as the time-space fractional convection-diffusion equation, see [44], in which for the solution a fast iterative method with a second-order implicit difference scheme was studied. Funding: This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.