Hexagonal Grid Computation of the Derivatives of the Solution to the Heat Equation by Using Fourth-Order Accurate Two-Stage Implicit Methods

: We give fourth-order accurate implicit methods for the computation of the ﬁrst-order spatial derivatives and second-order mixed derivatives involving the time derivative of the solution of ﬁrst type boundary value problem of two dimensional heat equation. The methods are constructed based on two stages: At the ﬁrst stage of the methods, the solution and its derivative with respect to time variable are approximated by using the implicit scheme in Buranay and Arshad in 2020. Therefore, O (cid:0) h 4 + τ (cid:1) of convergence on constructed hexagonal grids is obtained that the step sizes in the space variables x 1 , x 2 and in time variable are indicated by h , √ 32 h and τ , respectively. Special difference boundary value problems on hexagonal grids are constructed at the second stages to approximate the ﬁrst order spatial derivatives and the second order mixed derivatives of the solution. Further, O (cid:0) h 4 + τ (cid:1) order of uniform convergence of these schemes are shown for r = ωτ h 2 ≥ 116 , ω > 0. Additionally, the methods are applied on two sample problems.


Introduction
The modeling of numerous phenomena in diverse scientific fields leads us to consider conventional or fractional boundary value problems of time dependent differential equations on a modeling domain such as the first and second type boundary value problems to heat equation or diffusion equation. For example, the Brownian motion problem in statistics is modeled by heat equation via the Fokker-Planck equation (Adriaan Fokker [1] and Max Planck [2]). It is also named as the Kolmogorov forward equation, who discovered the concept in 1931, see in [3] independently. The stock market fluctuations represent one of the several important real-world applications of the mathematical model of Brownian motion. It was first given in the PhD thesis titled as "The theory of speculation", by Louis Bachelier (see Mandelbrot and Hudson [4]) in 1900.
Another representative sample of problems that mathematical modeling brings about the heat equation is the image processing problems appearing through many applied sciences from archaeology to zoology. Examples of archaeological investigations include a camcorder for 3D underwater reconstruction of archeological objects in the study of Meline et al. [5]. Furthermore, a recent investigation by Woźniak and Polap [6] gave soft trees with neural components as image processing technique for archeological excavations. In zoology, a study of image reconstruction problem by the application of magnetic resonance imaging was given by Ziegler et al. [7] and in medical sciences as medical image reconstruction was studied in Zeng [8]. Furthermore, tomography, and medical and industrial applications are archetypal examples where substantial mathematical manipulation is required. In some cases, the aim is humble denoising or de-blurring. Witkin [9] and Koenderink [10] gave the modeling of blurring of an image by the heat equation. Later, a problem of solving the reverse heat equation known as de-blurring is studied in Rudin et al. [11] and Guichard and Morel [12].
Additionally, in mathematical biology, Wolpert [13,14] gave a phenomenological concept of pattern formation and differentiation known as positional information. The pre-programming of the cells for reacting to a chemical concentration and differentiate accordingly, into different kinds of cells such as cartilage cells was proposed. Afterwards, the animal coat patterns, pattern formation on growing domains as alligators, snakes and bacterial patterns were modeled by reaction diffusion equations in Murray [15]. Furthermore, therein, gliomas or glioblastomas, which are highly diffusive brain tumors, are analyzed and a mathematical model for the spatiotemporal dynamics of tumor growth was developed. Therefore, the basic model in dimensional form was given by the diffusion equation where c(x, t) is the number of cells at a position x and time t, ρ represents the net rate of growth of cells including proliferation and death (or loss), and J diffusional flux of cells taken J = D∇c, where D(x) (distance 2 /time) is the diffusion coefficient of cells in brain tissue and ∇ is the gradient operator. In general, finding analytical solutions of these modeled problems is a difficult task or even not possible. Approximations are needed when a mathematical model is switched to a numerical model. Finite difference methods (FDM) are a class of numerical techniques for solving differential equations that each derivative appearing in the partial differential equation has to be replaced by a suitable divided difference of function values at the chosen grid points, see Grossman et al. [16]. In the last decade, the use of advanced computers has led to the widespread use of FDM in modern numerical analysis. For example, recently, a study on fractional diffusion equation-based image denoising model using Crank-Nicholson and Grünwald Letnikov difference schemes (CN-GL) have been given in Abirami et al. [17]. Another example is the most recent investigation by Buranay and Nouman [18] in which computation of the solution to heat equation on special polygons, where ω > 0 and f is the heat source by using implicit schemes defined on hexagonal grids was given. Therein, under some smoothness assumptions of the solution, two implicit methods were developed both on two layers with 14-point that has convergence orders of O h 2 + τ 2 and O h 4 + τ accordingly to the solution on the grids. Besides the solution of the modeled problem, the high accurate computation of the derivatives of the solution are fundamental to determine some important phenomena of the considered model problem. Such as for the diffusion problem (1) the functions ∂c ∂t and J gives the rate of change of the cells and diffusional flux of cells, respectively.
In the literature, exhaustive studies exist for the approximation of the derivatives of the solution to Laplace's equation under some smoothness conditions of the boundary functions and compatibility conditions. For the 2D Laplace equation, research was conducted by Volkov [19] and Dosiyev and Sadeghi [20]. For the 3D Laplace equation on a rectangular parallelepiped, studies were given by Volkov [21] and Dosiyev and Sadeghi [22], and recently by Dosiyev and Abdussalam [23], and Dosiyev and Sarikaya [24].
For the heat equation, the derivative of the solution of one-dimensional heat equation with respect to the space variable was given in Buranay and Farinola [25]. Within this paper, two implicit schemes were developed that converge to the corresponding exact spatial derivative with O h 2 + τ and O h 2 + τ 2 accordingly. Most recently, in Buranay et al. [26] numerical methods using implicit schemes with hexagonal grids approximating the derivatives of the solution of (2) on a rectangle has been given. The smoothness condition x,t , 0 < α < 1 in the Hölder space was required and uniform convergence on the grids to the respective spatial derivatives of O h 2 + τ 2 of accuracy for r = ωτ h 2 ≤ 3 7 was proved.
In regard to the equilateral triangulation with a regular hexagonal support, we remark the research by Barrera et al. [27] where a new class of quasi-interpolant was constructed which has remarkable properties such as high order of regularity and polynomial reproduction. Furthermore, on the Delaunay triangulation, we mention the study by Guessab [28] that approximations of differentiable convex functions on arbitrary convex polytopes were given. Further, optimal approximations were computed by using efficient algorithms accessed by the set of barycentric coordinates generated by the Delaunay triangulation.
The motivation of the contributions of this research is the need of highly accurate and time-efficient numerical algorithms that compute the derivatives of the solution u(x 1 , x 2 , t) to the heat Equation (2). The achievements of this study are summarized as follows. 1.
The first type boundary value problem (Dirichlet problem) for the heat Equation (2) on a rectangle D is considered. The smoothness of the solution u is taken from the and Q T = D × (0, T) also D, Q T denote the closure of D, Q T , respectively. At the first stage, an implicit scheme on hexagonal grids given in Buranay and Nouman [18] with O h 4 + τ order of accuracy is used to approximate the solution u(x 1 , x 2 , t).
The step sizes h and √ 3 2 h are taken for the spatial variables x 1 and x 2 , respectively, while τ is taken for the time variable. An analogous implicit method is also given to approximate the derivative of the solution with respect to time.

2.
At the second stages, computation of the first-order spatial derivatives and secondorder mixed derivatives involving time derivatives of the solution u(x 1 , x 2 , t) of (2) are developed. When r = ωτ h 2 ≥ 1 16 uniform convergence of the approximate derivative to the exact derivatives ∂u ∂x i , ∂u ∂t , and ∂ 2 u ∂x i ∂t , i = 1, 2 with order O h 4 + τ of accuracy on the hexagonal grids are proved.

3.
Numerical examples are given and for the solution of the obtained algebraic linear systems preconditioned conjugate gradient method is used. The incomplete block matrix factorization of the M-matrices given in Buranay and Iyikal [29] (see also Concus et al. [30], Axelsson [31]) is applied for the preconditioning.

Hexagonal Grid Approximation of the Heat Equation and the Rate of Change by Using Fourth Order Accurate Difference Schemes
where we require a 2 to be multiple of √ 3. Next, let γ j , j = 1, 2, 3, 4, be the sides of D that starting from the side x 1 = 0 are labeled in anticlockwise direction. Furthermore, the boundary of D is shown by S = 4 j=1 γ j . Further, we indicate the closure of D by D = D ∪ S. Let x = (x 1 , x 2 ) and is the closure of Q T . Let s be a non-integer positive number, C s, s 2 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 ∂ ξ+s 1 +s 2 u ∂t ξ ∂x with bounded norm where u (j) u (s) further, u α x , u β t for α, β ∈ (0, 1) are defined as Volkov gave the differentiability properties of solutions of boundary value problems for the Laplace and Poisson equations on rectangle in the study [32]. On cylindrical domains with smooth boundary, the differentiability properties of solutions of the parabolic equations were given in Ladyženskaja et al. [33] and Friedman [34]. On regions with edges, Azzam and Kreyszig studied the smoothness of solutions of parabolic equations for the Dirichlet problem in [35] and for the mixed boundary value problem in [36].

Dirichlet Problem of Heat Equation and Difference Problem: Stage 1 H 4th (u)
Our interest is the following problem for the heat equation: where ω is positive constant. This problem is framed as first type (Dirichlet) boundary value problem.
We assume that the initial and boundary functions ϕ(x 1 , x 2 ), φ(x 1 , x 2 , t), respectively, also the heat source function f (x 1 , x 2 , t) possess the necessary smoothness and satisfy the conditions that the problem (11) has unique solution u ∈ C 9+α, 9+α 2 x,t Q T . We define hexagonal grids on D with the step size h, such that h = a 1 /N 1 , and N 1 is positive integer and present this set by D h as Let γ h j , j = 1, ..., 4 be the set of nodes on the interior of γ j , and let γ h Further, we denote by D * lh and D * rh the set of interior nodes whose distance from the boundary is h 2 , thus the hexagon is irregular hexagon with a ghost point that emerges through the left (x 1 = 0) or right (x 1 = a 1 ) side of the rectangle, respectively. The illustration of the exact solution at the irregular hexagons with a ghost point at time levels t − τ, t and t + τ is given in Figure 1.
Next, we give the set of interior hexagonal points and the lateral surface points by respectively. Let D * lh γ τ = D * lh × γ τ ⊂ D h γ τ and D * rh γ τ = D * rh × γ τ ⊂ D h γ τ and D * h γ τ = D * lh γ τ ∪ D * rh γ τ . Furthermore, D 0h γ τ = D h γ τ \D * h γ τ and D h γ τ is the closure of D h γ τ . We denote the center of the hexagon by P 0 and Patt(P 0 ) is the pattern of the hexagon consisting the neighboring points P i , i = 1, ..., 6. Furthermore, the exact solution at the neighboring points P i , i = 0, 1, ..., 6 for the time moment t + τ is presented by u k+1 is the value on the boundary point as given: where ( p, x 2 , t + τ) ∈ S h T and the value of p = 0 if P 0 ∈ D * lh γ τ and p = a 1 if P 0 ∈ D * rh γ τ . Moreover, u k+1 h,τ,P i , i = 0, ..., 6, u k+1 h,τ,P A , present the numerical solution at the same space coordinates of P i , i = 0, ..., 6 and P A accordingly for time moments t + τ. We also use the following notations in Table 1 to denote the values and partial derivatives of the heat source function f and f t = ∂ f ∂t with respect to the space variables.
For computing numerically the solution of the BVP(u) we use the following difference problem given in Buranay and Arshad [18] and call this Stage 1 H 4th (u) .
where ϕ, φ are the initial and boundary functions in (11), respectively, also and where and ϕ, φ are the initial and boundary functions in (11).
x,t Q T , fourth-order accurate implicit schemes for the solution of the BVP ∂u ∂t is proposed with the following difference problem. This stage is called (24), respectively, and

M−Matrices and Convergence of Finite Difference Schemes in Stage 1 H 4th (u)
and . Analogous notation is also used for the vectors. Further, let w be a vector with coordinates w j , j = 1, 2, ..., N, the vector with coordinates w j is denoted by |w|. For a fixed time level k ≥ 0 we present the Equations (17) and (28) in matrix form with N unknown interior grid points L j , j = 1, 2, ..., N, labeled using standard ordering as respectively, where A, B ∈ R N×N and u k , q k u , u k t , q k u t ∈ R N and and Inc is the neighboring topology matrix,Ȇ 1 ,Ȇ 2 are diagonal matrices with entries respectively (see Buranay and Arshad [18]). (31) are symmetric positive definite (spd) matrices

Lemma 1. (Buranay and Arshad [18]) (a) The matrices A and B in
Proof. Taking into consideration Lemma 1, the matrix A is a spd matrix. Further, using the Equations (32)-(35), A is strictly diagonally dominant matrix with positive diagonal entries. Furthermore, off-diagonal entries are non-positive for r = ωτ From (17) and (36) the error function (36) satisfies the following system as given in Buranay and Arshad: [18] Θ h,τ ξ u,k+1 where and ψ, ψ * and φ are presented in (17). Analogously, using (28) and (37) the error function (37) satisfies the following system: where and φ t , ψ t , and ψ * t are the given functions in (27), (29) and (30) respectively. Further, the following systems are considered: for k = 0, ..., M − 1, where κ k 1 , κ k 2 and κ k 1 , κ k 2 are given functions. The algebraic systems (44) and (45) at a fixed time level k ≥ 0 may be given in matrix representation as accordingly. In these equations, w k , w k , κ k , κ k ∈ R N and the matrices A and B are given in (32).
Additionally, let Theorem 1. For the solution of the system (38) and (41) when r = ωτ h 2 ≥ 1 16 , the following pointwise error estimations hold true: respectively, where and and u is the solution of BVP(u) and ρ(x 1 , x 2 , t) is the function giving the distance from the considered hexagonal grid point (x 1 , x 2 , t) ∈ D h γ τ to the surface of Q T .
Proof. We give the proof of (57) by considering the auxiliary system and the majorant functions which ξ u l (x 1 , x 2 , t), satisfy the following difference problem for l = 1, 2, 3, respectively.
Therefore, difference problems (62) and (66) in matrix form are accordingly, and A and B are as given in (32) and on D 0h γ τ , and 5 6 Ω 1 (h, τ) ≥ Ψ u,k 2 on D * h γ τ and on the basis of Lemma 3 we obtain The proof of (58) is analogous and follows from Lemma 3 by taking the majorant functions where Ω t,1 (h, τ) is as given in (60).

Second Stages of the Implicit Methods Approximating
and the corresponding sets of grid points is shown by S h T γ i , i = 1, 2, ..., 5. on S T γ i , i = 1, 2, ..., 5 and use the next problem given in Buranay et al. [26].

Boundary Value Problem for
where, f (x 1 , x 2 , t) is the given heat source function in (11) and x,t Q T . Further, we take where ϕ(x 1 , x 2 ), φ(x 1 , x 2 , t) are as in (11), and u h,τ is obtained by using Stage 1 H 4th (u) .

Lemma 4.
Let u be the solution of BVP(u) in (11) and u h,τ be the solution of (17) in Stage 1 H 4th (u) . Then, it holds that where Ω 1 (h, τ) in (59) and d in (61) was defined.

Lemma 5.
Let u h,τ be the solution of the problem (17) in Stage 1 H 4th (u) . Then, it holds that and Ω 1 (h, τ) in (59) and d in (61) was defined.

Theorem 2.
The solution v h,τ of the finite difference problem given in Stage 2 H 4th ∂u for r = ωτ h 2 ≥ 1 16 where θ, σ are as given in (93)
Proof. The proof basically is analogous with the proof of Theorem 2 and follows from the Let the BVP(u) be given. First, we apply Stage 1 H 4th (u) and obtain the approximate solution u h,τ on the hexagonal grids. Then, by denoting q i = ∂u ∂x 2 on S T γ i , i = 1, 2, ..., 5 we use the next problem for z = ∂u ∂x 2 , proposed in Buranay et al. [26] Boundary Value Problem for ∂u

Second Stages of the Implicit Methods Approximating
We take and ϕ(x 1 , x 2 ), φ(x 1 , x 2 , t) given in (11) are the initial and boundary functions, respectively, u h,τ is the solution taken by using Stage 1 H 4th (u) .

Lemma 6.
Let u be the solution of (11) and u h,τ be the approximation achieved by using Stage 1 H 4th (u) . Then, the following inequality holds true for r ≥ 1 16 where, Ω 1 (h, τ) is given in (59) and d is defined in (61).
and u h,τ be the approximation taken by using Stage 1 H 4th (u) .
Then, the following inequality is true: where Ω 1 (h, τ) is given in (59) and d is defined in (61).
Proof. The proof is analogous to the proof of Theorem 4, and follows from the requirement x,t Q T .

Experimental Investigations
The proposed fourth order two stage implicit methods are applied on two test problems such that for the first example the exact solution is known. However, for the second example the exact solution is not given. We take D = (x 1 , , and t ∈ [0, 1]. Further, Mathematica is used for the realization of the algorithms in machine precision. Also we used preconditioned conjugate gradient method with the preconditioning approach given in Buranay and Iyikal [29] (see also Concus et al. [30] and Axelsson [31]). We define the following:

while the error function resulting by the methods
Further, we denote the order of convergence of the approximate solution v 2 −µ ,2 −λ and z 2 −µ ,2 −λ to the functions v = ∂u ∂x 1 and z = ∂u ∂x 2 obtained by using the fourth-order implicit Furthermore, the order of convergence of the approximate solutions v t,2 −µ ,2 −λ and z t,2 −µ ,2 −λ to their corresponding exact solutions v t = ∂ 2 u ∂x 1 ∂t and z t = ∂ 2 u ∂x 2 ∂t obtained by H 4th ∂ 2 u ∂x i ∂t , i = 1, 2 are given by We remark that the computed values of (161) and (162) are ≈ 2 4 showing the fourth order convergence of the given methods in x 1 , x 2 and linear convergence in t.

Conclusions
Numerical methods using implicit schemes defined on hexagonal grids are proposed for computing the derivatives of the solution to Dirichlet problem of heat equation on rectangle. For the required smoothness conditions of the solution and when r = ωτ h 2 ≥ 1 16 the uniform convergence of the constructed difference schemes on the grids to the respective exact derivatives ∂u ∂x i and ∂ 2 u ∂x i ∂t , i = 1, 2 is shown to be O h 4 + τ . Novelty Statement: In Buranay et al. [26], we gave a second-order hexagonal grid approximation of the first-order spatial derivatives of the solution to BVP(u) in (11) with the smoothness condition u ∈ C 7+α, 7+α 2 x,t , 0 < α < 1 in the Hölder space. The method was established in two stages. In this study, we require that u ∈ C 9+α, 9+α 2 x,t , and give hexagonal grid computation of all the first order derivatives and the mixed order second order derivatives involving the time derivative by developing two stage implicit methods of fourth order accurate in space variables.