Next Article in Journal
On Strongly Continuous Resolving Families of Operators for Fractional Distributed Order Equations
Next Article in Special Issue
Finite-Approximate Controllability of Riemann–Liouville Fractional Evolution Systems via Resolvent-Like Operators
Previous Article in Journal
The Fractional Derivative of the Dirac Delta Function and Additional Results on the Inverse Laplace Transform of Irrational Functions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

Department of Mathematics, Faculty of Arts and Sciences, Eastern Mediterranean University, Famagusta, North Cyprus, Mersin 10, Turkey
*
Author to whom correspondence should be addressed.
Submission received: 30 December 2020 / Revised: 13 February 2021 / Accepted: 20 February 2021 / Published: 26 February 2021

Abstract

:
The first 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 first order derivatives of the solution with respect to the spatial variables. To approximate the solution at the first stage, the unconditionally stable two layer implicit method on hexagonal grids given by Buranay and Arshad in 2020 is used which converges with O h 2 + τ 2 of accuracy on the grids. Here, 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. At the second stage, we propose special difference boundary value problems on hexagonal grids for the approximation of first derivatives with respect to spatial variables of which the boundary conditions are defined by using the obtained solution from the first stage. It is proved that the given schemes in the difference problems are unconditionally stable. Further, for r = ω τ h 2 3 7 , uniform convergence of the solution of the constructed special difference boundary value problems to the corresponding exact derivatives on hexagonal grids with order O h 2 + τ 2 is shown. Finally, the method is applied on a test problem and the numerical results are presented through tables and figures.

1. 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
u t = ω 2 u x 1 2 + 2 u x 2 2 + f x 1 , x 2 , t ,
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 x , t 6 + α , 3 + α 2 , 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 Q ¯ T = D ¯ × 0 , T the solution u x 1 , x 2 , t belongs to the Hölder space C x , t 7 + α , 7 + α 2 Q ¯ T , 0 < α < 1 , where x = x 1 , x 2 D ¯ , t 0 , 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 Section 3 and Section 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.

2. First Type Heat Problem and Second Order Accurate Solution on Hexagonal Grids

Let D = x 1 , x 2 : 0 < x 1 < a 1 , 0 < x 2 < a 2 be a rectangle, where a 2 is multiple of 3 and let γ j , j = 1 , 2 , 3 , 4 , be the sides of D enumerated in anticlockwise direction starting from the side x 1 = 0 . Further, S = j = 1 4 γ j is the boundary of D and denote by D ¯ = D S the closure of D. Let Q T = D × 0 , T , 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):
u t = ω 2 u x 1 2 + 2 u x 2 2 + f x 1 , x 2 , t on Q T ,
u x 1 , x 2 , 0 = φ x 1 , x 2 on D ¯ ,
u x 1 , x 2 , t = ϕ x 1 , x 2 , t on S T ,
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 x , t 7 + α , 7 + α 2 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
D h = x = ( x 1 , x 2 ) D : x 1 = i j 2 h , x 2 = 3 ( i + j ) 2 h , i = 1 , 2 , ; j = 0 ± 1 ± 2 , .
Let γ j h , j = 1 , , 4 be the set of nodes on the interior of γ j and let γ ^ j h = γ j 1 γ j be the jth vertex of D, S h = j = 1 4 ( γ j h γ ^ j h ) , D ¯ h = D h S h . Further, let D * l h , D * r h 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 * l h D * r h and D 0 h = D h \ D * h .
Next, let
γ τ = t k = k τ , τ = T M , k = 1 , , M ,
γ ¯ τ = t k = k τ , τ = T M , k = 0 , , M ,
and the set of internal nodes and lateral surface nodes be defined by
D h γ τ = D h × γ τ = x , t : x = ( x 1 , x 2 ) D h , t γ τ ,
S T h = S h × γ ¯ τ = x , t : x = ( x 1 , x 2 ) S h , t γ ¯ τ ,
respectively. Let D * l h γ τ = D * l h × γ τ D h γ τ and D * r h γ τ = D * r h × γ τ D h γ τ and D * h γ τ = D * l h γ τ D * r h γ τ . In addition, D 0 h γ τ = D h γ τ \ D * h γ τ and D h γ τ ¯ is the closure of D h γ τ .
Let P 0 denote the center of the hexagon and P a t t P 0 denote the pattern of the hexagon consisting the neighboring points P i , i = 1 , , 6 . In addition, u P i k + 1 denotes the exact solution at the point P i and u P A k + 1 denotes the value at the boundary point for the time moment t + τ as follows:
u P 1 k + 1 = u ( x 1 h 2 , x 2 + 3 2 h , t + τ ) , u P 3 k + 1 = u ( x 1 h 2 , x 2 3 2 h , t + τ ) , u P 2 k + 1 = u ( x 1 h , x 2 , t + τ ) , u P 5 k + 1 = u ( x 1 + h , x 2 , t + τ ) , u P 4 k + 1 = u ( x 1 + h 2 , x 2 3 2 h , t + τ ) , u P 6 k + 1 = u ( x 1 + h 2 , x 2 + 3 2 h , t + τ ) , u P 0 k + 1 = u ( x 1 , x 2 , t + τ ) , u P A k + 1 = u ( p ^ , x 2 , t + τ ) , ( p ^ , x 2 , t + τ ) S T h ,
where the value of p ^ = 0 if P 0 D * l h γ τ and p ^ = a 1 if P 0 D * r h γ τ as also given in (21). Analogously, the values u P i k , i = 0 , , 6 and u P A k 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 h , τ , P i k + 1 , i = 0 , , 6 , u h , τ , P A k + 1 , and u h , τ , P i k , i = 0 , , 6 , u h , τ , P A k 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 P 0 k + 1 2 = f ( x 1 , x 2 , t + τ 2 ) , and f P A k + 1 = f ( p ^ , x 2 , t + τ ) . 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 2 n d of the two stage implicit method:
Stage 1 H 2 n d
Θ h , τ u h , τ k + 1 = Λ h , τ u h , τ k + ψ on D 0 h γ τ ,
Θ h , τ * u h , τ k + 1 = Λ h , τ * u h , τ k + Γ h , τ * ϕ + ψ * on D * h γ τ ,
u h , τ = φ x 1 , x 2 , t = 0 on D ¯ h ,
u h , τ = ϕ x 1 , x 2 , t on S T h ,
for k = 0 , , M 1 , where
ψ = f P 0 k + 1 2 ,
ψ * = f P 0 k + 1 2 1 6 f P A k + 1 2 ,
Θ h , τ u k + 1 = 1 τ + 2 ω h 2 u P 0 k + 1 ω 3 h 2 i = 1 6 u P i k + 1 ,
Λ h , τ u k = 1 τ 2 ω h 2 u P 0 k + ω 3 h 2 i = 1 6 u P i k ,
Θ h , τ * u k + 1 = 1 τ + 7 ω 3 h 2 u P 0 k + 1 ω 3 h 2 u ( p + η , x 2 , t + τ ) + u ( p , x 2 + 3 2 h , t + τ ) + u ( p , x 2 3 2 h , t + τ ) ,
Λ h , τ * u k = 1 τ 7 ω 3 h 2 u P 0 k + ω 3 h 2 u ( p , x 2 + 3 2 h , t ) + u ( p , x 2 3 2 h , t ) + u ( p + η , x 2 , t ) ,
Γ h , τ * ϕ = 2 ω 9 h 2 ϕ ( p ^ , x 2 + 3 2 h , t + τ ) + ϕ ( p ^ , x 2 3 2 h , t + τ ) + ϕ ( p ^ , x 2 + 3 2 h , t ) + ϕ ( p ^ , x 2 3 2 h , t ) + 1 6 τ + 8 ω 9 h 2 ϕ ( p ^ , x 2 , t + τ ) + 1 6 τ + 8 ω 9 h 2 ϕ ( p ^ , x 2 , t ) ,
and
p = h , p ^ = 0 , η = h 2 if P 0 D * l h γ τ , p = a 1 h , p ^ = a 1 , η = h 2 if P 0 D * r h γ τ .
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:
A u ˜ k + 1 = B u ˜ k + τ q u k ,
where A , B R N × N are
A = I + ω τ h 2 C , B = I ω τ h 2 C ,
and
C = D 1 1 3 I n c R N × N ,
and u ˜ k , q u k R N . The matrix I n c 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
d 1 , j j = 2 if L j D 0 h γ τ 7 3 if L j D * h γ τ , j = 1 , 2 , , N .
Lemma 1.
(a) The matrix A in (22) is a nonsingular M-matrix and is also a symmetric positive definite matrix.
(b) A 1 2 < 1 and A 1 B 2 < 1 for r = ω τ h 2 > 0 .
Proof. 
Proof is given in ([29]). □
Let
ε h , τ u = u h , τ u on D h γ τ ¯ .
From (10)–(13) and (26), the error function ε h , τ u satisfies the following system as given in [29]
Θ h , τ ε h , τ u , k + 1 = Λ h , τ ε h , τ u , k + Ψ 1 u , k on D 0 h γ τ ,
Θ h , τ * ε h , τ u , k + 1 = Λ h , τ * ε h , τ u , k + Ψ 2 u , k on D * h γ τ ,
ε h , τ u = 0 , t = 0 on D ¯ h ,
ε h , τ u = 0 on S T h ,
where
Ψ 1 u , k = Λ h , τ u k Θ h , τ u k + 1 + ψ ,
Ψ 2 u , k = Λ h , τ * u k Θ h , τ * u k + 1 + Γ h , τ * ϕ + ψ * ,
and ψ , ψ * , and ϕ are the given functions in (10), (11), and (13), respectively.

Pointwise Priory Estimation For the Error Function (27)–(30)

Consider the following systems
Θ h , τ q ^ h , τ k + 1 = Λ h , τ q ^ h , τ k + g ^ 1 k on D 0 h γ τ ,
Θ h , τ * q ^ h , τ k + 1 = Λ h , τ * q ^ h , τ k + Γ h , τ * q ^ ϕ , h , τ + g ^ 2 k on D * h γ τ ,
q ^ h , τ = q ^ φ , h , τ , t = 0 on D ¯ h ,
q ^ h , τ = q ^ ϕ , h , τ on S T h ,
Θ h , τ q ¯ h , τ k + 1 = Λ h , τ q ¯ h , τ k + g ¯ 1 k on D 0 h γ τ ,
Θ h , τ * q ¯ h , τ k + 1 = Λ h , τ * q ¯ h , τ k + Γ h , τ * q ¯ ϕ , h , τ + g ¯ 2 k on D * h γ τ ,
q ¯ h , τ = q ¯ φ , h , τ , t = 0 on D ¯ h ,
q ¯ h , τ = q ¯ ϕ , h , τ on S T h ,
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
A q ^ k + 1 = B q ^ k + τ g ^ k ,
A q ¯ k + 1 = B q ¯ k + τ g ¯ k ,
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 K 1 K 2 0 is nonpositive and K 1 K 2 means that K 1 K 2 0 is nonnegative wherever they present matrices in R N × N or vectors in R N .
Lemma 2.
Let q ^ k + 1 be the solution of the difference Equation (41) and q ¯ k + 1 be the solution of the difference Equation (42). For r = ω τ h 2 3 7 , if
q ¯ 0 0 and g ¯ k 0 ,
q ^ 0 q ¯ 0 ,
g ^ k g ¯ k ,
for k = 0 , , M 1 , then
q ¯ k + 1 0 and q ^ k + 1 q ¯ k + 1 for k = 0 , , M 1 .
Proof. 
On the basis of Lemma 1, A 1 0 and if r = ω τ h 2 3 7 then B 0 and from (43) we have g ¯ k 0 , k = 0 , , M 1 and q ¯ 0 0 . Then, assume that q ¯ k 0 by using induction we have
q ¯ k + 1 = A 1 B q ¯ k + τ A 1 g ¯ k 0 ,
which gives q ¯ k + 1 0 , k = 0 , , M 1 . In addition, q ^ 0 q ¯ 0 from (44). Next assume that q ^ k q ¯ k , by using (45) and induction gives
q ^ k + 1 = A 1 B q ^ k + τ A 1 g ^ k ,
q ^ k + 1 A 1 B q ^ k + τ A 1 g ^ k A 1 B q ¯ k + τ A 1 g ¯ k = q ¯ k + 1 .
Thus, we obtain (46). □
Let
S T γ 1 = γ 1 × 0 , T = 0 , x 2 , t : 0 , x 2 γ 1 , t 0 , T , S T γ 2 = γ 2 × 0 , T = x 1 , 0 , t : x 1 , 0 γ 2 , t 0 , T , S T γ 3 = γ 3 × 0 , T = a 1 , x 2 , t : a 1 , x 2 γ 3 , t 0 , T , S T γ 4 = γ 4 × 0 , T = x 1 , a 2 , t : x 1 , a 2 γ 4 , t 0 , T , S T γ 5 = x 1 , x 2 , 0 : x 1 , x 2 D ¯ , t = 0 ,
and S T h γ i , i = 1 , 2 , , 5 denote the corresponding sets of grid points. In addition, let F = i = 1 5 S T γ i denote the surface of Q T .
Theorem 1.
For the solution of the problem (27)–(30), the following inequality holds true
ε h , τ u d Ω 1 h , τ ρ x 1 , x 2 , t , on D h γ τ ¯ ,
for r = ω τ h 2 3 7 where
Ω 1 h , τ = 1 24 τ 2 ( 1 + 6 ω ) β * + 3 ω 10 h 2 α * ,
α * = max max Q ¯ T 4 u x 1 4 , max Q ¯ T 4 u x 2 4 , max Q ¯ T 4 u x 1 2 x 2 2 ,
β * = max max Q ¯ T 3 u t 3 , max Q ¯ T 4 u x 2 2 t 2 , max Q ¯ T 4 u x 1 2 t 2 ,
d = max a 1 2 ω , a 2 2 ω , 1 ,
and u is the exact solution ofBVP(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
Θ h , τ ε ^ h , τ u , k + 1 = Λ h , τ ε ^ h , τ u , k + Ω 1 h , τ on D 0 h γ τ ,
Θ h , τ * ε ^ h , τ u , k + 1 = Λ h , τ * ε ^ h , τ u , k + 5 6 Ω 1 h , τ on D * h γ τ
ε ^ h , τ u = ε ^ φ , h , τ u = 0 , t = 0 on D ¯ h ,
ε ^ h , τ u = ε ^ ϕ , h , τ u = 0 on S T h ,
and the majorant functions
ε ¯ 1 u x 1 , x 2 , t = 1 2 ω Ω 1 h , τ a 1 x 1 x 1 2 0 on D h γ τ ¯ ,
ε ¯ 2 u x 1 , x 2 , t = 1 2 ω Ω 1 h , τ a 2 x 2 x 2 2 0 on D h γ τ ¯ ,
ε ¯ 3 u x 1 , x 2 , t = Ω 1 h , τ t 0 on D h γ τ ¯ ,
in which ε ¯ i u , i = 1 , 2 , 3 satisfy the difference boundary value problem
Θ h , τ ε ¯ i , h , τ u , k + 1 = Λ h , τ ε ¯ i , h , τ u , k + Ω 1 h , τ on D 0 h γ τ ,
Θ h , τ * ε ¯ i , h , τ u , k + 1 = Λ h , τ * ε ¯ i , h , τ u , k + Γ h , τ * ε ¯ i , ϕ , h , τ u + 5 6 Ω 1 h , τ on D * h γ τ ,
ε ¯ i , h , τ u = ε ¯ i , φ , h , τ u = ε ¯ i u x 1 , x 2 , 0 0 , t = 0 on D ¯ h ,
ε ¯ i , h , τ u = ε ¯ i , ϕ , h , τ u 0 on S T h .
Therefore, we write the difference problems (56)–(59) and (63)–(66) for fixed k 0 in matrix form
A ε ^ u , k + 1 = B ε ^ u , k + τ e ^ u , k ,
A ε ¯ i u , k + 1 = B ε ¯ i u , k + τ e ¯ i u , k , i = 1 , 2 , 3 ,
respectively, where A and B are the matrices given in (22) and e ¯ i u , k , ε ¯ i u , k , i = 1 , 2 , 3 and ε ^ u , k , e ^ u , k R N . Using (52) and (56)–(66) gives ε ¯ i u , 0 0 , ε ^ u , 0 ε ¯ i u , 0 , and e ¯ i u , k 0 , and e ^ u , k e ¯ i u , k , i = 1 , 2 , 3 , for k = 0 , , M 1 . On the basis of Lemma 2, we get ε ^ u , k + 1 ε ¯ i u , k + 1 , k = 0 , , M 1 and using that Ω 1 h , τ Ψ 1 u , k on D 0 h γ τ , and 5 6 Ω 1 h , τ Ψ 2 u , k on D * h γ τ gives
ε h , τ u min i = 1 , 2 , 3 ε ¯ i u x 1 , x 2 , t d Ω 1 h , τ ρ x 1 , x 2 , t on D h γ τ ¯ ,

3. Difference Problem Approximating u x 1 on Hexagonal grids with O ( h 2 + τ 2 ) Order of Accuracy

We use the notations x 1 f P 0 k + 1 2 = f x 1 ( x 1 , x 2 , t + τ 2 ) and x 1 f P A k + 1 2 = f x 1 ( p ^ , x 2 , t + τ 2 ) . Given the boundary value Problems (2)–(4), we denote p i = u x 1 on S T γ i , i = 1 , 2 , , 5 and setup the next boundary value problem for v = u x 1 .
BVP(2):
L v = f x 1 , x 2 , t x 1 on Q T ,
v x 1 , x 2 , t = p i on S T γ i , i = 1 , 2 , , 5 ,
where
L t ω 2 x 1 2 + 2 x 2 2 ,
and f x 1 , x 2 , t is the given function in (2). On the basis of the assumption that u C x , t 7 + α , 7 + α 2 Q ¯ T , we assume that the solution v C x , t 6 + α , 3 + α 2 Q ¯ T .
We take
p 1 h 2 n d = 1 2 h 3 u ( 0 , x 2 , t ) + 4 u h , τ h , x 2 , t u h , τ 2 h , x 2 , t if P 0 D 0 h γ τ 1 3 h 8 u ( 0 , x 2 , t ) + 9 u h , τ h 2 , x 2 , t u h , τ 3 h 2 , x 2 , t if P 0 D * l h γ τ on S T h γ 1 ,
p 3 h 2 n d = 1 2 h 3 u ( a 1 , x 2 , t ) 4 u h , τ a 1 h , x 2 , t + u h , τ a 1 2 h , x 2 , t if P 0 D 0 h γ τ 1 3 h 8 u ( a 1 , x 2 , t ) 9 u h , τ a 1 h 2 , x 2 , t + u h , τ a 1 3 h 2 , x 2 , t if P 0 D * r h γ τ on S T h γ 3 ,
p i h = ϕ x 1 , x 2 , t x 1 on S T h γ i , i = 2 , 4 ,
p 5 h = φ x 1 , x 2 x 1 on S T h γ 5 ,
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 2 n d .
Lemma 3.
The following inequality
p i h 2 n d u h , τ p i h 2 n d u 3 d Ω 1 h , τ , i = 1 , 3 .
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 inStage 1 H 2 n d and Ω 1 h , τ and d are as given in (52) and (55), respectively.
Proof. 
Taking into consideration Theorem 1, and using (51), (73), and (74) when P 0 D 0 h γ τ , we have
p i h 2 n d u h , τ p i h 2 n d u 1 2 h 4 h d Ω 1 h , τ + 2 h d Ω 1 h , τ 3 d Ω 1 h , τ , i = 1 , 3 if P 0 D 0 h γ τ ,
where Ω 1 is as given in (52) and d is the constant given in (55). When P 0 D * h γ τ results
p i h 2 n d u h , τ p i h 2 n d u 1 3 h 9 h 2 d Ω 1 h , τ + 3 h 2 d Ω 1 h , τ 2 d Ω 1 h , τ , i = 1 , 3 if P 0 D * h γ τ .
Thus, we obtain (77). □
Lemma 4.
The following inequality
max S T h γ 1 S T h γ 3 p i h 2 n d u h , τ p i M 1 h 2 + 3 d Ω 1 h , τ , i = 1 , 3 ,
holds true for r = ω τ h 2 3 7 where u h , τ is the solution of the difference problem inStage 1 H 2 n d and M 1 = 1 3 max Q ¯ T 3 u x 1 3 and Ω 1 and d are as given in (52) and (55), respectively.
Proof. 
From the assumption that the exact solution u C x , t 7 + α , 7 + α 2 Q ¯ T , at the end points 0 , η 3 2 h , k τ S T h γ 1 and a 1 , η 3 2 h , k τ S T h γ 3 of each line segment
x 1 , η 3 2 h , k τ : 0 x 1 a 1 , 0 x 2 = η 3 2 h a 2 , 0 t = k τ T ,
difference formulas (73) and (74) give the second order approximation of u x 1 , respectively. From the truncation error formula (see [39]) it follows that
max S T h γ 1 S T h γ 3 p i h 2 n d u p i h 2 3 max Q ¯ T 3 u x 1 3 , i = 1 , 3 if P 0 D 0 h γ τ .
Analogously,
max S T h γ 1 S T h γ 3 p i h 2 n d u p i h 2 8 max Q ¯ T 3 u x 1 3 , i = 1 , 3 if P 0 D * h γ τ .
Using Lemma 3 and the estimations (81) and (82) follows (80). □
We construct the following difference problem for the numerical solution of BVP(2) and denote this stage as
Stage 2 H 2 n d u x 1
Θ h , τ v h , τ k + 1 = Λ h , τ v h , τ k + D x 1 ψ on D 0 h γ τ ,
Θ h , τ * v h , τ k + 1 = Λ h , τ * v h , τ k + Γ h , τ * p 1 h 2 n d + D x 1 ψ * on D * l h γ τ ,
Θ h , τ * v h , τ k + 1 = Λ h , τ * v h , τ k + Γ h , τ * p 3 h 2 n d + D x 1 ψ * on D * r h γ τ ,
v h , τ = p i h 2 n d u h , τ on S T h γ i , i = 1 , 3 ,
v h , τ = p i h on S T h γ i , i = 2 , 4 , 5 ,
where p 1 h 2 n d , p 3 h 2 n d , and p i h , i = 2 , 4 , 5 are defined by (73)–(76) and the operators Θ h , τ , Λ h , τ , Θ h , τ * , Λ h , τ * and Γ h , τ * are the operators given in (16)–(20), respectively. Additionally,
D x 1 ψ = x 1 f P 0 k + 1 2 ,
D x 1 ψ * = x 1 f P 0 k + 1 2 1 6 x 1 f P A k + 1 2 ,
Let
ε h , τ v = v h , τ v on D h γ τ ¯ ,
where v = u x 1 . From (83)–(87) and (90), we have
Θ h , τ ε h , τ v , k + 1 = Λ h , τ ε h , τ v , k + Ψ 1 v , k on D 0 h γ τ ,
Θ h , τ * ε h , τ v , k + 1 = Λ h , τ * ε h , τ v , k + Γ h , τ * ε h , τ * v + Ψ 2 v , k on D * h γ τ ,
ε h , τ v = 0 on S T h γ i , i = 2 , 4 , 5 ,
ε h , τ v = ε h , τ * v = p i h 2 n d u h , τ p i on S T h γ i , i = 1 , 3 ,
where
Ψ 1 v , k = Λ h , τ v k Θ h , τ v k + 1 + D x 1 ψ ,
Ψ 2 v , k = Λ h , τ * v k Θ h , τ * v k + 1 + Γ h , τ * p i + D x 1 ψ * , i = 1 , 3 .
Let
θ 1 = max max Q ¯ T 4 v x 1 4 , max Q ¯ T 4 v x 2 4 , max Q ¯ T 4 v x 1 2 x 2 2 , σ 1 = max max Q ¯ T 3 v t 3 , max Q ¯ T 4 v x 2 2 t 2 , max Q ¯ T 4 v x 1 2 t 2 ,
and
θ = max θ 1 , 40 M 1 3 + 12 d ω α * ,
σ = max σ 1 , 3 d β * ,
where α * , β * are as given in (53), (54), respectively, and M 1 is as given in (80).
Theorem 2.
The implicit scheme given inStage 2 ( u x 1 ) is unconditionally stable.
Proof. 
The obtained algebraic linear system of Equations (83)–(87) can be given in matrix form:
A v ˜ k + 1 = B v ˜ k + τ q v k ,
k = 0 , 1 , , M 1 , where A and B are the matrices given in (22) and v ˜ k , q v k R N . On the basis of the assumption that the exact solution v of the BVP(2) belongs to C x , t 6 + α , 3 + α 2 Q ¯ T and using Lemma 1 and by induction we get
v ˜ k + 1 2 A 1 B 2 v ˜ k 2 + τ A 1 2 q v k 2 v ˜ 0 2 + τ k = 0 k q v k 2 .
Thus, Lax and Richtmyer sufficient condition for stability given in Theorem 1 of [33] is satisfied and the scheme is unconditionally stable. □
Theorem 3.
The solution v h , τ of the finite difference problem given inStage 2 H 2 n d u x 1 satisfies
max D h γ τ ¯ v h , τ v σ 12 1 + 6 ω T + 1 τ 2 + 3 θ 40 h 2 1 + a 1 2 + a 2 2 ,
for r = ω τ h 2 3 7 where θ , σ are as given in (97), (98), respectively, and v = u x 1 is the exact solution ofBVP(2).
Proof. 
Consider the auxiliary system
Θ h , τ ε ^ h , τ v , k + 1 = Λ h , τ ε ^ h , τ v , k + Ω 2 x 1 on D 0 h γ τ ,
Θ h , τ * ε ^ h , τ v , k + 1 = Λ h , τ * ε ^ h , τ v , k + Γ h , τ * ε ^ h , τ v * + Ω 2 x 1 1 6 Ω 2 p ^ on D * h γ τ ,
ε ^ h , τ v = 0 on S T h γ i , i = 2 , 4 , 5 ,
ε ^ h , τ v = ε ^ h , τ v * = p i h 2 n d u h , τ p i on S T h γ i , i = 1 , 3 ,
where
Ω 2 x 1 = σ 24 a 1 1 + 6 ω τ 2 2 a 1 x 1 + 3 θ ω 10 h 2 , σ 24 1 + 6 ω τ 2 + 3 θ ω 10 h 2 Ψ 1 v , k ,
Ω 2 x 1 1 6 Ω 2 p ^ = 1 + 6 ω τ 2 5 72 h 48 a 1 + θ ω 4 h 2 if P 0 D * l h γ τ 1 + 6 ω τ 2 5 144 + h 48 a 1 + θ ω 4 h 2 if P 0 D * r h γ τ Ψ 2 v , k ,
and x 1 = h 2 and p ^ = 0 if P 0 D * l h γ τ and x 1 = a 1 h 2 , p ^ = a 1 if P 0 D * r h γ τ . We take the majorant function
ε ¯ v x 1 , x 2 , t = ε ¯ 1 v x 1 , x 2 , t + ε ¯ 2 v x 1 , x 2 , t ,
where
ε ¯ 1 v x 1 , x 2 , t = σ τ 2 24 a 1 1 + 6 ω t + 1 2 a 1 x 1 0 on D h γ τ ¯ ,
ε ¯ 2 v x 1 , x 2 , t = 3 θ 40 h 2 1 + a 1 2 + a 2 2 x 1 2 x 2 2 0 on D h γ τ ¯ ,
The function in (108) satisfies the difference problem
Θ h , τ ε ¯ h , τ v , k + 1 = Λ h , τ ε ¯ h , τ v , k + Ω 2 x 1 on D 0 h γ τ ,
Θ h , τ * ε ¯ h , τ v , k + 1 = Λ h , τ * ε ¯ h , τ v , k + Γ h , τ * ε ¯ h , τ v * + Ω 2 x 1 1 6 Ω 2 p ^ on D * h γ τ ,
ε ¯ h , τ v = ε ¯ h , τ v * = ε ¯ 1 v 0 , x 2 , t + ε ¯ 2 v 0 , x 2 , t on S T h γ 1 ,
ε ¯ h , τ v = ε ¯ 1 v x 1 , 0 , t + ε ¯ 2 v x 1 , 0 , t on S T h γ 2 ,
ε ¯ h , τ v = ε ¯ h , τ v * = ε ¯ 1 v a 1 , x 2 , t + ε ¯ 2 v a 1 , x 2 , t on S T h γ 3 ,
ε ¯ h , τ v = ε ¯ 1 v x 1 , a 2 , t + ε ¯ 2 v x 1 , a 2 , t on S T h γ 4 ,
ε ¯ h , τ v = ε ¯ 1 v x 1 , x 2 , 0 + ε ¯ 2 v x 1 , x 2 , 0 on S T h γ 5 .
The algebraic system of Equations (102)–(105) and (111)–(117) can be written in matrix form as
A ε ^ v , k + 1 = B ε ^ v , k + τ e ^ v , k ,
A ε ¯ v , k + 1 = B ε ¯ v , k + τ e ¯ v , k ,
respectively, for k = 0 , , M 1 , where A , B are matrices as given in (22) and ε ^ v , k , ε ¯ v , k , e ^ v , k , e ¯ v , k R N . Using (106)–(117), we have ε ¯ v , 0 0 , and e ¯ v , k 0 , and e ^ v , k e ¯ v , k for k = 0 , , M 1 , and ε ^ v , 0 ε ¯ v , 0 . Then, on the basis of Lemma 2, we get ε ^ v , k + 1 ε ¯ v , k + 1 for k = 0 , , M 1 . From
ε ¯ v x 1 , x 2 , t ε ¯ v 0 , 0 , T = σ 12 1 + 6 ω T + 1 τ 2 + 3 θ 40 h 2 1 + a 1 2 + a 2 2 ,
and using (106) and (107) follows (101). □

4. Second Order Approximation of u x 2 on Hexagonal grids

Let the BVP(1) be given, then we use the notations x 2 f P 0 k + 1 2 = f x 2 ( x 1 , x 2 , t + τ 2 ) and x 2 f P A k + 1 2 = f x 2 ( p ^ , x 2 , t + τ 2 ) 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 x 2 ,
BVP(3):
L z = f x 1 , x 2 , t x 2 on Q T ,
z x 1 , x 2 , t = q i on S T γ i , i = 1 , 2 , , 5 ,
where L is the operator in (72) and f x 1 , x 2 , t is the given function in (2). We assume that the solution z C x , t 6 + α , 3 + α 2 Q ¯ T and take
q 2 h 2 n d = 1 2 3 h 3 u ( x 1 , 0 , t ) + 4 u h , τ x 1 , 3 h , t
u h , τ x 1 , 2 3 h , t on S T h γ 2 , q 4 h 2 n d = 1 2 3 h 3 u ( x 1 , a 2 , t ) 4 u h , τ x 1 , a 2 3 h , t
+ u h , τ x 1 , a 2 2 3 h , t on S T h γ 4 ,
q i h = ϕ x 1 , x 2 , t x 2 on S T h γ i , i = 1 , 3 ,
q 5 h = φ x 1 , x 2 x 2 on S T h γ 5 ,
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 2 n d .
Lemma 5.
The following inequality holds
q i h 2 n d u h , τ q i h 2 n d u 3 d Ω 1 h , τ , i = 2 , 4 ,
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) inStage 1 H 2 n d and Ω 1 h , τ and d are as given in (52) and (55), respectively.
Proof. 
Taking into consideration Theorem 1, and using (51), (122), and (123), we have
q i h 2 n d u h , τ q i h 2 n d u 1 2 3 h 4 3 h d Ω 1 h , τ + 2 3 h d Ω 1 h , τ 3 d Ω 1 h , τ , i = 2 , 4 ,
thus, we obtain (126). □
Lemma 6.
The following inequality is true
max S T h γ 2 S T h γ 4 q i h 2 n d u h , τ q i M 2 h 2 + 3 d Ω 1 h , τ , i = 2 , 4 ,
for r = ω τ h 2 3 7 , where M 2 = max Q ¯ T 3 u x 2 3 and u h , τ is the solution of the difference problem inStage 1 H 2 n d and Ω 1 h , τ and d are as given in (52) and (55), respectively.
Proof. 
From the assumption that the exact solution u C x , t 7 + α , 7 + α 2 Q ¯ T , at the end points ϑ h , 0 , k τ S T h γ 2 and ϑ h , a 2 , k τ S T h γ 4 of each line segment
ϑ h , x 2 , k τ : 0 x 1 = ϑ h a 1 , 0 x 2 a 2 , 0 t = k τ T ,
difference formulas (122) and (123) give the second order approximation of u x 2 , respectively. From the truncation error formula (see [39]), it follows that
max S T h γ 2 S T h γ 4 q i h 2 n d u q i h 2 max Q ¯ T 3 u x 2 3 , i = 2 , 4 .
Taking M 2 = max Q ¯ T 3 u x 2 3 and using Lemma 5 and the estimation (127) and (129) follows (128). □
We construct the following difference problem for the numerical solution of BVP(3) and denote this stage as
Stage 2 H 2 n d u x 2
Θ h , τ z h , τ k + 1 = Λ h , τ z h , τ k + D x 2 ψ on D 0 h γ τ ,
Θ h , τ * z h , τ k + 1 = Λ h , τ * z h , τ k + Γ h , τ * q 1 h + D x 2 ψ * on D * l h γ τ ,
Θ h , τ * z h , τ k + 1 = Λ h , τ * z h , τ k + Γ h , τ * q 3 h + D x 2 ψ * on D * r h γ τ ,
z h , τ = q i h 2 n d u h , τ on S T h γ i , i = 2 , 4 ,
z h , τ = q i h on S T h γ i , i = 1 , 3 , 5 ,
where q 2 h 2 n d , q 4 h 2 n d , and q i h , i = 1 , 3 , 5 are defined by (122)–(125) and the operators Θ h , τ , Λ h , τ , Θ h , τ * , Λ h , τ * , and Γ h , τ * are the operators given in (16)–(20), respectively. In addition,
D x 2 ψ = x 2 f P 0 k + 1 2 ,
D x 2 ψ * = x 2 f P 0 k + 1 2 1 6 x 2 f P A k + 1 2 .
Let
ε h , τ z = z h , τ z on D h γ τ ¯ .
From (130)–(134) and (137), we have
Θ h , τ ε h , τ z , k + 1 = Λ h , τ ε h , τ z , k + Ψ 1 z , k on D 0 h γ τ ,
Θ h , τ * ε h , τ z , k + 1 = Λ h , τ * ε h , τ z , k + Ψ 2 z , k on D * h γ τ ,
ε h , τ z = 0 on S T h γ i , i = 1 , 3 , 5 ,
ε h , τ z = q i h 2 n d u h , τ q i on S T h γ i , i = 2 , 4 ,
where q i h are defined by (122)–(125) and
Ψ 1 z , k = Λ h , τ z k Θ h , τ z k + 1 + D x 2 ψ ,
Ψ 2 z , k = Λ h , τ * z k Θ h , τ * z k + 1 + Γ h , τ * q i + D x 2 ψ * , i = 1 , 3 .
Let
κ 1 = max max Q ¯ T 4 z x 1 4 , max Q ¯ T 4 z x 2 4 , max Q ¯ T 4 z x 1 2 x 2 2 ,
δ 1 = max max Q ¯ T 3 z t 3 , max Q ¯ T 4 z x 2 2 t 2 , max Q ¯ T 4 z x 1 2 t 2 ,
and
κ = max κ 1 , 40 M 2 3 + 12 d ω α * ,
δ = max δ 1 , 3 d β * ,
α * , β * are as given in (53), (54), respectively, and M 2 is the constant given in Lemma 6 and z is the solution of BVP(3).
Theorem 4.
The implicit scheme given inStage 2 H 2 n d u x 2 is unconditionally stable.
Proof. 
The obtained algebraic linear system of Equations (130)–(134) can be given in matrix form:
A z ˜ k + 1 = B z ˜ k + τ q z k ,
for k = 0 , 1 , , M 1 , where, A , B are as given in (22) and z ˜ k , q z k R N . On the basis of the assumption that the exact solution z of the BVP(3) belongs to C x , t 6 + α , 3 + α 2 Q ¯ T and using Lemma 1 and induction we get
z ˜ k + 1 2 A 1 B 2 z ˜ k 2 + τ A 1 2 q z k 2 z ˜ 0 2 + τ k = 0 k q z k 2 .
Therefore, the scheme is unconditionally stable. □
Theorem 5.
The solution z h , τ of the finite difference problem given inStage 2 H 2 n d u x 2 satisfies
max D h γ τ ¯ z h , τ z δ 12 1 + 6 ω T + 1 τ 2 + 3 κ 40 1 + a 1 2 + a 2 2 h 2 ,
for r = ω τ h 2 3 7 , where κ , δ are as given in (146), (147) respectively and z = u x 2 is the exact solution ofBVP(3).
Proof. 
Consider the auxiliary system
Θ h , τ ε ^ h , τ z , k + 1 = Λ h , τ ε ^ h , τ z , k + Ω 3 x 2 on D 0 h γ τ ,
Θ h , τ * ε ^ h , τ z , k + 1 = Λ h , τ * ε ^ h , τ z , k + 5 6 Ω 3 x 2 on D * h γ τ
ε ^ h , τ z = 0 on S T h γ i , i = 1 , 3 , 5 ,
ε ^ h , τ z = q i h 2 n d u h , τ q i on S T h γ i , i = 2 , 4 ,
where q 2 h 2 n d , q 4 h 2 n d , q i h , i = 1 , 3 , 5 , are defined by (122)–(125) and
Ω 3 x 2 = δ 24 a 2 1 + 6 ω τ 2 2 a 2 x 2 + 3 κ ω 10 h 2
δ 24 1 + 6 ω τ 2 + 3 κ ω 10 h 2 Ψ 1 z , k , 5 6 Ω 3 x 2 = 5 δ 144 a 2 1 + 6 ω τ 2 2 a 2 x 2 + κ ω 4 h 2
5 δ 144 1 + 6 ω τ 2 + κ ω 4 h 2 Ψ 2 z , k .
We take the majorant function
ε ¯ z x 1 , x 2 , t = ε ¯ 1 z x 1 , x 2 , t + ε ¯ 2 z x 1 , x 2 , t ,
where
ε ¯ 1 z x 1 , x 2 , t = δ 24 a 2 τ 2 1 + 6 ω ( t + 1 ) 2 a 2 x 2 0 on D h γ τ ¯ ,
ε ¯ 2 z x 1 , x 2 , t = 3 κ 40 h 2 1 + a 1 2 + a 2 2 x 1 2 x 2 2 0 on D h γ τ ¯ .
The majorant function in (157) satisfies the difference problem
Θ h , τ ε ¯ h , τ z , k + 1 = Λ h , τ ε ¯ h , τ z , k + Ω 3 x 2 on D 0 h γ τ ,
Θ h , τ * ε ¯ h , τ z , k + 1 = Λ h , τ * ε ¯ h , τ z , k + Γ h , τ * ε ¯ h , τ z * + 5 6 Ω 3 x 2 on D * h γ τ ,
ε ¯ h , τ z = ε ¯ h , τ z * = ε ¯ 1 z 0 , x 2 , t + ε ¯ 2 z 0 , x 2 , t on S T h γ 1 ,
ε ¯ h , τ z = ε ¯ 1 z x 1 , 0 , t + ε ¯ 2 z x 1 , 0 , t on S T h γ 2 ,
ε ¯ h , τ z = ε ¯ h , τ z * = ε ¯ 1 z a 1 , x 2 , t + ε ¯ 2 z a 1 , x 2 , t on S T h γ 3 ,
ε ¯ h , τ z = ε ¯ 1 z x 1 , a 2 , t + ε ¯ 2 z x 1 , a 2 , t on S T h γ 4 ,
ε ¯ h , τ z = ε ¯ 1 z x 1 , x 2 , 0 + ε ¯ 2 z x 1 , x 2 , 0 on S T h γ 5 .
We write the algebraic system of Equations (151)–(154) and (160)–(166) for fixed k 0 in matrix form
A ε ^ z , k + 1 = B ε ^ z , k + τ e ^ z , k ,
A ε ¯ z , k + 1 = B ε ¯ z , k + τ e ¯ z , k ,
respectively, where A , B are as given in (22) and ε ^ z , k , ε ¯ z , k , e ^ z , k , e ¯ z , k R N . Using (155)–(156), we get e ¯ z , k 0 and e ^ z , k e ¯ z , k for k = 0 , 1 , , M 1 and ε ¯ z , 0 0 , ε ^ z , 0 ε ¯ z , 0 . Then, on the basis of Lemma 2 follows ε ^ z , k + 1 ε ¯ z , k + 1 , k = 0 , 1 , , M 1 . From
ε ¯ z x 1 , x 2 , t ε ¯ z 0 , 0 , T = δ 12 1 + 6 ω T + 1 τ 2 + 3 κ 40 1 + a 1 2 + a 2 2 h 2 ,
and using (155), (156) follows (150). □

5. 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 , x 2 : 0 < x 1 < 1 , 0 < x 2 < 3 2 , 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:
H 2 n d u x 1
denotes the proposed two stage implicit method on hexagonal grids for the approximation of the derivative u x 1 .
H 2 n d u x 2
denotes the proposed two stage implicit method on hexagonal grids for the approximation of the derivative u x 2 .
C T H 2 n d u x 1
presents the Central Processing Unit time in seconds ( C P U s ) per time level for the method H 2 n d u x 1 .
C T H 2 n d u x 2
presents the Central Processing Unit time in seconds ( C P U s ) per time level for the method H 2 n d u x 2 .
T C T H 2 n d u x 1
shows the total Central Processing Unit time in seconds required for the solution at t = 1 , by the method H 2 n d u x 1 .
T C T H 2 n d u x 2
shows the total Central Processing Unit time in seconds required for the solution at t = 1 , by the method H 2 n d u x 2 .
The proposed two stage implicit method for the approximation of the derivatives u x 1 , u x 2 are denoted as the methods H 2 n d u x 1 , and H 2 n d u x 2 . Additionally, the corresponding solutions are denoted by v 2 μ , 2 λ , and z 2 μ , 2 λ , respectively, for h = 2 μ and τ = 2 λ where μ , λ are positive integers. On the grid points D h γ τ ¯ , which is the closure of D h γ τ we present the error function ε h , τ obtained by H 2 n d u x 1 , and H 2 n d u x 2 by ε H 2 n d u x 1 and by ε H 2 n d u x 2 , respectively. In addition, maximum norm of the errors max D h γ τ ¯ ε H 2 n d u x 1 and max D h γ τ ¯ ε H 2 n d u x 2 are presented by ε H 2 n d u x 1 and ε H 2 n d u x 2 , accordingly. Further, we denote the order of convergence of the approximate solution v 2 μ , 2 λ to the exact solution v = u x 1 obtained by using the two stage implicit method H 2 n d u x 1 by
H 2 n d u x 1 = ε H 2 n d u x 1 ( 2 μ , 2 λ ) ε H 2 n d u x 1 ( 2 μ + 1 , 2 λ + 1 ) .
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 2 n d u x 2 is given by
H 2 n d u x 2 = ε H 2 n d u x 2 ( 2 μ , 2 λ ) ε H 2 n d u x 2 ( 2 μ + 1 , 2 λ + 1 ) .
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:
u t = 0.5 2 u x 1 2 + 2 u x 2 2 + f x 1 , x 2 , t on Q T , u x 1 , x 2 , 0 = 0.0001 x 1 57 8 1 x 1 + cos x 2 57 8 3 2 x 2 on D ¯ , u x 1 , x 2 , t = u ^ x 1 , x 2 , t on S T ,
where
f x 1 , x 2 , t = 0.00035625 t 41 16 6.125 x 1 41 8 + 8.125 x 1 49 8 + 3 3249 912 7.125 x 2 x 2 49 4 cos x 2 57 8 + 3 2793 912 8.125 x 2 x 2 41 8 sin x 2 57 8 u ^ x 1 , x 2 , t = 0.0001 t 57 16 + x 1 57 8 1 x 1 + cos x 2 57 8 3 2 x 2 ,
are the heat source and exact solution. Table 1 demonstrates C T H 2 n d u x 1 , T C T H 2 n d u x 1 , maximum norm of the errors for h = 2 μ , μ = 4 , 5 , 6 , 7 when τ = 2 λ , λ = 13 , 14 , 15 , 16 , that is r = 0.5 τ h 2 3 7 and the order of convergence of v h , τ to the exact derivatives v = u x 1 with respect to h and τ obtained by using the constructed two stage implicit method H 2 n d u x 1 . Table 2 shows C T H 2 n d u x 2 , T C T H 2 n d u x 2 , maximum norm of the errors for the same pairs of h , τ as in Table 1 and the order of convergence of z h , τ to the exact derivative z = u x 2 with respect to h and τ obtained by using the constructed two stage implicit method H 2 n d u x 2 . Table 1 and Table 2 justify the theoretical results given such that the approximate solutions v h , τ and z h , τ of the proposed method converge to the corresponding exact derivatives v = u x 1 and z = u x 2 with second order both in the spatial variables x 1 , x 2 and the time variable t for r 3 7 .
Table 3 presents the C T H 2 n d u x 1 , T C T H 2 n d u x 1 , maximum norm of the errors for h = 2 μ , μ = 4 , 5 , 6 , 7 , 8 when τ = 2 λ , λ = 8 , 9 , 10 , 11 , 12 , that is r = 0.5 τ h 2 > 3 7 and the order of convergence of v h , τ to the exact derivative v = u x 1 with respect to h and τ obtained by using the constructed two stage implicit method H 2 n d u x 1 . Table 4 shows the C T H 2 n d u x 2 , T C T H 2 n d 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 τ obtained by using the constructed two stage implicit method H 2 n d u x 2 . Numerical results given in Table 3 and Table 4 demonstrate that when r > 3 7 , the approximate solutions v h , τ and z h , τ of the proposed method also converge with second order both in the spatial variables x 1 , x 2 and the time variable t to their corresponding exact derivatives v = u x 1 and z = u x 2 .
Figure 3 illustrates the absolute error functions ε H 2 n d u x 1 ( 2 4 , 2 13 ) , ε H 2 n d u x 1 ( 2 5 , 2 14 ) , ε H 2 n d u x 1 ( 2 6 , 2 15 ) , and ε H 2 n d u x 1 ( 2 7 , 2 16 ) at time moment t = 0.2 obtained by using H 2 n d u x 1 . Figure 4 demonstrates the absolute error functions ε H 2 n d u x 2 ( 2 4 , 2 13 ) , ε H 2 n d u x 2 ( 2 5 , 2 14 ) , ε H 2 n d u x 2 ( 2 6 , 2 15 ) , and ε H 2 n d u x 2 ( 2 7 , 2 16 ) at time moment t = 0.2 obtained by using H 2 n d u x 2 . The exact derivative v = u x 1 and the grid function v 2 6 , 2 15 for h = 2 6 , τ = 2 15 at time moment t = 0.2 obtained by using H 2 n d u x 1 are presented in Figure 5. Further, Figure 6 shows the exact derivative z = u x 2 and grid function z 2 6 , 2 15 at time moment t = 0.2 obtained by using H 2 n d u x 2 .
Table 5 shows the C T H 2 n d u x 1 , T C T H 2 n d u x 1 , maximum norm of the errors for r 3 7 , and the order of convergence of v h , τ to the exact derivative v = u x 1 with respect to h and τ obtained when third order approximations for v = u x 1
p 1 h 3 r d = 1 6 h 11 u ( 0 , x 2 , t ) + 18 u h , τ h , x 2 , t 9 u h , τ 2 h , x 2 , t + 2 u h , τ 3 h , x 2 , t if P 0 D 0 h γ τ 1 60 h 184 u ( 0 , x 2 , t ) + 225 u h , τ h 2 , x 2 , t 50 u h , τ 3 h 2 , x 2 , t + 9 u h , τ 5 h 2 , x 2 , t if P 0 D * l h γ τ on S T h γ 1 ,
p 3 h 3 r d = 1 6 h 11 u ( a 1 , x 2 , t ) 18 u h , τ a 1 h , x 2 , t + 9 u h , τ a 1 2 h , x 2 , t 2 u h , τ a 1 3 h , x 2 , t if P 0 D 0 h γ τ 1 60 h 184 u ( a 1 , x 2 , t ) 225 u h , τ a 1 h 2 , x 2 , t + 50 u h , τ a 1 3 h 2 , x 2 , t 9 u h , τ a 1 5 h 2 , x 2 , t if P 0 D * r h γ τ on S T h γ 3 ,
are used on S T h γ i , i = 1 , 3 for the Stage 2 H 2 n d u x 1 . Table 6 shows C T H 2 n d u x 2 , T C T H 2 n d u x 2 , maximum norm of the errors for r 3 7 and the order of convergence of z h , τ to the exact derivative z = u x 2 with respect to h and τ obtained when third order approximations for z = u x 2 .
q 2 h 3 r d = 1 6 3 h 11 u ( x 1 , 0 , t ) + 18 u h , τ x 1 , 3 h , t 9 u h , τ x 1 , 2 3 h , t + 2 u h , τ x 1 , 3 3 h , t on S T h γ 2 ,
q 4 h 3 r d = 1 6 3 h 11 u ( x 1 , a 2 , t ) 18 u h , τ x 1 , a 2 3 h , t + 9 u h , τ x 1 , a 2 2 3 h , t 2 u h , τ x 1 , a 2 3 3 h , t on S T h γ 4 ,
are used on S T h γ i , i = 2 , 4 for the Stage 2 H 2 n d u x 2 . Table 7 presents C T H 2 n d u x 1 , T C T H 2 n d u x 1 , maximum norm of the errors for r > 3 7 , and the order of convergence of v h , τ to the exact derivatives v = u x 1 with respect to h and τ obtained by using the difference formulae (172), (173) on S T h γ i , i = 1 , 3 for the Stage 2 H 2 n d u x 1 . Table 8 gives C T H 2 n d u x 2 , T C T H 2 n d u x 2 , maximum norm of the errors for r > 3 7 , and the order of convergence of z h , τ to the exact derivative z = u x 2 with respect to h and τ obtained by using the difference formulae (174), (175) on S T h γ i , i = 2 , 4 for the Stage 2 H 2 n d u x 2 . Numerical results given in Table 5, Table 6, Table 7 and Table 8 demonstrate that the approximate solution v h , τ and z h , τ of the proposed method converges to the corresponding exact derivatives v = u x 1 and z = u x 2 with second order both in the spatial variables x 1 , x 2 and the time variable t with better error ratios.

6. Discussions and Conclusion

A second order accurate two stage implicit method on hexagonal grids is proposed for the approximation of the first order derivatives of the solution to first type boundary value problem of the heat equation with respect to spatial variables on rectangle. At the first stage, for the error function, we provided a pointwise prior estimation depending on ρ x 1 , x 2 , t , which is the distance from the current grid point in the domain to the surface of Q T . At the second stage, we constructed special difference problems for the approximation of u x 1 and u x 2 . Uniform convergence with order O h 2 + τ 2 of the solution of the constructed difference problems to the corresponding exact derivatives u x 1 and u x 2 on the hexagonal grids when r = ω τ h 2 3 7 is proved. Furthermore, the given two stage implicit method 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.

Author Contributions

Conceptualization, S.C.B.; methodology, S.C.B.; software, N.A. and A.H.M.; validation, S.C.B., A.H.M., and N.A.; formal analysis, S.C.B., A.H.M., and N.A; investigation, A.H.M. and N.A.; resources, S.C.B., A.H.M., and N.A.; writing—review and editing, S.C.B., A.H.M., and N.A.; visualization, S.C.B., A.H.M., and N.A.; supervision, S.C.B.; project administration, S.C.B. All authors have read and approved the final manuscript.

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Doungmo Goufo, E.F. Evolution equations with a parameter and application to transport-convection differential equations. Turk. J. Math. 2017, 41, 636–654. [Google Scholar] [CrossRef]
  2. Doungmo Goufo, E.F. Solvability of chaotic fractional systems with 3D four-scroll attractors. Chaos Solitons Fractals 2017, 104, 443–451. [Google Scholar] [CrossRef]
  3. Doungmo Goufo, E.F. The Proto-Lorenz system in its chaotic fractional and fractal structure. Int. J. Bifurcation Chaos 2020, 30, 2050180. [Google Scholar] [CrossRef]
  4. Doungmo Goufo, E.F.; Khan, Y. A new auto-replication in systems of attractors with two and three merged basins of attraction via control. Commun. Nonlinear Sci. Numer. Simulat. 2021, 96, 105709. [Google Scholar] [CrossRef]
  5. Owolabi, K.M.; Hammouch, Z. Mathematical modeling and analysis of two-variable system with noninteger-order derivative. Chaos 2019, 29, 013145. [Google Scholar] [CrossRef] [PubMed]
  6. Doungmo Goufo, E.F.; Khan, Y.; Mugisha, S. Control parameter & solutions to generalized evolution equations of stationarity, relaxation and diffusion. Results Phys. 2018, 9, 1502–1507. [Google Scholar]
  7. Doungmo Goufo, E.F.; Kumar, S.; Mugisha, S.B. Similarities in a fifth-order evolution equation with and with no singular kernel. Chaos Solitons Fractals 2020, 130, 109467. [Google Scholar] [CrossRef]
  8. Heydari, M.H.; Hooshmandasl, M.R.; Mohammadi, F.; Cattani, C. Wavelets method for solving systems of nonlinear singular fractional Volterra integro-differential equations. Commun. Nonlinear Sci. Numer. Simulat. 2014, 19, 37–48. [Google Scholar] [CrossRef]
  9. Čermák, J.; Kisela, T. Delay-dependent stability switches in fractional differential equations. Commun. Nonlinear Sci. Numer. Simulat. 2019, 79, 104888. [Google Scholar] [CrossRef]
  10. Tuladhar, R.; Santamaria, F.; Stamova, I. Fractional Lotka-Volterra-Type Cooperation Models: Impulsive control on their stability behavior. Entropy 2020, 22, 970. [Google Scholar] [CrossRef]
  11. Zhang, H.; Sandu, A. A second-order diagonally-implicit-explicit multi-stage integration method. Procedia Computer Sci. 2012, 9, 1039–1046. [Google Scholar] [CrossRef] [Green Version]
  12. Hasan, M.K.; Ahamed, M.S.; Alam, M.S.; Hossain, M.B. An implicit method for numerical solution of singular and stiff initial value problems. J. Comput. Eng. 2013, 2013, 1–5. [Google Scholar] [CrossRef] [Green Version]
  13. Ashyralyev, A.; Hincal, E.; Kaymakamzade, B. A Crank Nicolson difference scheme for system of nonlinear observing epidemic models. In AIP Conference Proceedings; AIP Publishing LLC: Melville, NY, USA, 2019; Volume 2183. [Google Scholar]
  14. Ahmed, N.; Ali, M.; Baleanu, D.; Rafiq, M.; ur Rehman, M.A. Numerical analysis of diffusive susceptible-infected-recovered epidemic model in three space dimension. Chaos Solitons Fractals 2020, 132, 1–10. [Google Scholar] [CrossRef]
  15. Sadourney, R.; Arakawa, A.; Mintz, Y. Integration of the nondivergent barotropic vorticity equation with an icosahedral-hexagonal grid for the sphere. Mon. Wea. Rev. 1968, 96, 351–356. [Google Scholar] [CrossRef] [Green Version]
  16. Williamson, D. Integration of the barotropic vorticity equation on a spherical geodesic grid. Tellus 1968, 20, 642–653. [Google Scholar] [CrossRef]
  17. Sadourney, R. Numerical integration of the primitive equation on a spherical grid with hexagonal cells. In Proceedings of the WMO/IUGG Symposium on Numerical Whether Prediction, Tokyo, Japan, 26 November–4 December 1969; pp. 45–77. [Google Scholar]
  18. Sadourney, R.; Morel, P. A finite-difference approximation of the primitive equations for a hexagonal grid on a plane. Mon. Wea. Rev. 1969, 97, 439–445. [Google Scholar] [CrossRef]
  19. Masuda, Y. A finite difference scheme by making use of hexagonal mesh-points. In Proceedings of the WMO/IUGG Symposium on Numerical Whether Prediction, Tokyo, Japan, 26 November–4 December 1969; pp. VII35–VII44. [Google Scholar]
  20. Masuda, Y.; Ohnishi, H. An integration scheme of the primitive equation model with an icosahedral-hexagonal grid system and its application to the shallow water equation. J. Meteorol. Soc. Jpn. 1986, 64, 317–326. [Google Scholar] [CrossRef] [Green Version]
  21. Thacker, W.C. Irregular grid finite difference techniques: Simulations of oscillations in shallow circular basins. J. Phys. Oceanogr. 1977, 7, 284–292. [Google Scholar] [CrossRef] [Green Version]
  22. Thacker, W.C. Comparison of finite element and finite difference schemes, Part II: Two dimensional gravity wave motion. J. Phys. Oceanogr. 1978, 8, 680–689. [Google Scholar] [CrossRef] [Green Version]
  23. Salmon, R.; Talley, D.L. Generalization of Arakawa’s Jacobian. J. Comput. Physics 1989, 83, 247–259. [Google Scholar] [CrossRef]
  24. Ničkovič, S. On the use of hexagonal grids for simulation of atmospheric processes. Contrib. Atmos. Phys. 1994, 67, 103–107. [Google Scholar]
  25. Ničkovič, S.; Gavrilov, M.B.; Tosič, I.A. Geostrophic adjuctment on hexagonal grids. Mon. Wea. Rev. 2001, 130, 668–683. [Google Scholar] [CrossRef]
  26. Pruess, K.; Bodvarsson, G.S. A Seven-Point Finite Difference method For Improved Grid Orientation Performance in Pattern Steam Floods; Lawrence Berkerly National Laboratory: Berkeley, CA, USA, 1983; pp. 1–32. [Google Scholar]
  27. Lee, D.; Tien, H.-C.; Luo, C.P.; Luk, H.-N. Hexagonal grid methods with applications to partial differential equations. Int. J. Comput. Math. 2014, 91, 1986–2009. [Google Scholar] [CrossRef]
  28. Richtmyer, R.D.; Morton, K.W. Difference Methods for Initial-Value Problems, 2nd ed.; Interscience Publishers a Division of Jhon Wiley and Sons: Hoboken, NJ, USA, 1967. [Google Scholar]
  29. Buranay, S.C.; Arshad, N. Hexagonal grid approximation of the solution of heat equation on special polygons. Adv. Differ. 2020, 2020, 309. [Google Scholar] [CrossRef]
  30. Arshad, N. Hexagonal Grid Approximation of the Solution of Two Dimensional Heat Equation. Ph.D. Thesis, Eastern Mediterranean University, Famagusta, Cyprus, 2020. [Google Scholar]
  31. Karaa, S. High-order approximation of 2D convection-diffusion equation on hexagonal grids. Numer. Methods Partial. Differ. Equations 2006, 22, 1238–1246. [Google Scholar] [CrossRef]
  32. Dosiyev, A.A.; Celiker, E. Approximation on the hexagonal grid of the Dirichlet problem for Laplace’s equation. Boundary Value Probl. 2014, 73, 1–19. [Google Scholar] [CrossRef] [Green Version]
  33. Lax, P.D.; Richtmyer, R.D. Survey of the stability of linear finite difference equations. Commun. Pure Applied Mathematics 1956, 9, 267–293. [Google Scholar] [CrossRef]
  34. Buranay, S.C.; Iyikal, O.C. Incomplete block-matrix factorization of M-matrices using two step iterative method for matrix inversion and preconditioning. Math. Methods Appl. Sci. 2020, in press. [Google Scholar] [CrossRef]
  35. Concus, P.; Golub, G.H.; Meurant, G. Block preconditioning for the conjugate gradient method. SIAM J. 1985, 6, 220–252. [Google Scholar] [CrossRef]
  36. Axelsson, O. A general incomplete block matrix factorization method. Linear Algebra Its Appl. 1986, 74, 179–190. [Google Scholar] [CrossRef] [Green Version]
  37. Azzam, A.; Kreyszig, E. On solutions of parabolic equations in regions with edges. Bull. Austral. Math. Soc. 1980, 22, 219–230. [Google Scholar] [CrossRef] [Green Version]
  38. Azzam, A.; Kreyszig, E. Smoothness of solutions of parabolic equations in regions with edges. Nagoya Math. J. 1981, 84, 159–168. [Google Scholar] [CrossRef] [Green Version]
  39. Burden, R.L.; Faires, J.D. Numerical Analysis, 9th ed.; Cengage Learning: Boston, MA, USA, 2011. [Google Scholar]
  40. Peaceman, D.W.; Rachford, H.H., Jr. The numerical solution of parabolic and elliptic differential equations. J. Soc. Industrial Appl. Math. 1955, 3, 28–41. [Google Scholar] [CrossRef]
  41. Douglas, J. On the numerical integration of 2 u x 2 + 2 u y 2 = u t by implicit methods. J. Soc. Industrial Appl. Math. 1955, 3, 42–65. [Google Scholar]
  42. Bagrinovskii, K.A.; Godunov, S.K. Difference schemes for multidimensional problems. Dokl. Akad. Nauk USSR 1957, 115, 431–433. [Google Scholar]
  43. Marchuk, G.I. Splitting and Alternating Methods. In Handbook of Numerical Analysis; Elsevier Science Publishers: Amsterdam, The Netherlands, 1990. [Google Scholar]
  44. Gu, X.M.; Huang, T.Z.; Ji, C.C.; Carpentieri, B.; Alikhanov, A.A. Fast iterative method with a second-order implicit difference scheme for time-space fractional convection–diffusion equation. J. Sci. Comput. 2017, 72, 957–985. [Google Scholar] [CrossRef]
Figure 1. The solution u P 2 k and u P 2 k + 1 on the left ghost points at time moments t = k τ and k + 1 τ , respectively.
Figure 1. The solution u P 2 k and u P 2 k + 1 on the left ghost points at time moments t = k τ and k + 1 τ , respectively.
Fractalfract 05 00019 g001
Figure 2. The solution u P 5 k and u P 5 k + 1 on the right ghost points at time moments t = k τ and k + 1 τ , respectively.
Figure 2. The solution u P 5 k and u P 5 k + 1 on the right ghost points at time moments t = k τ and k + 1 τ , respectively.
Fractalfract 05 00019 g002
Figure 3. The absolute error functions at time moment t = 0.2 obtained by using H 2 n d u x 1 .
Figure 3. The absolute error functions at time moment t = 0.2 obtained by using H 2 n d u x 1 .
Fractalfract 05 00019 g003
Figure 4. The absolute error functions at time moment t = 0.2 obtained by using H 2 n d u x 2 .
Figure 4. The absolute error functions at time moment t = 0.2 obtained by using H 2 n d u x 2 .
Fractalfract 05 00019 g004
Figure 5. The exact solution v = u x 1 and the approximate solution v 2 6 , 2 15 at t = 0.2 .
Figure 5. The exact solution v = u x 1 and the approximate solution v 2 6 , 2 15 at t = 0.2 .
Fractalfract 05 00019 g005
Figure 6. The exact solution z = u x 2 and the approximate solution z 2 6 , 2 15 at t = 0.2 .
Figure 6. The exact solution z = u x 2 and the approximate solution z 2 6 , 2 15 at t = 0.2 .
Fractalfract 05 00019 g006
Table 1. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when r = 0.5 τ h 2 3 7 .
Table 1. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when r = 0.5 τ h 2 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 1 TCT H 2 nd u x 1 ε H 2 nd u x 1 H 2 nd u x 1
2 6 2 4 , 2 13 0.03 197.34 9.3475 × 10 6 3.1457
2 5 2 5 , 2 14 0.09 1187.55 2.97147 × 10 6 3.5508
2 4 2 6 , 2 15 0.59 18 , 501.80 8.36840 × 10 7 3.7737
2 3 2 7 , 2 16 3.69 144 , 505.21 2.21757 × 10 7
Table 2. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when r = 0.5 τ h 2 3 7 .
Table 2. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when r = 0.5 τ h 2 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 2 TCT H 2 nd u x 2 ε H 2 nd u x 2 H 2 nd u x 2
2 6 2 4 , 2 13 0.02 181.88 3.72134 × 10 6 1.7362
2 5 2 5 , 2 14 0.13 1187.55 2.14336 × 10 6 2.6720
2 4 2 6 , 2 15 0.70 21 , 557.80 8.02154 × 10 7 3.2757
2 3 2 7 , 2 16 4.09 169 , 305.04 2.44880 × 10 7
Table 3. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when r = 0.5 τ h 2 > 3 7 .
Table 3. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when r = 0.5 τ h 2 > 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 1 TCT H 2 nd u x 1 ε H 2 nd u x 1 H 2 nd u x 1
2 1 2 4 , 2 8 0.02 4.75 9.34796 × 10 6 3.1458
1 2 5 , 2 9 0.08 37.30 2.97159 × 10 6 3.5508
2 2 6 , 2 10 0.42 347.70 8.36871 × 10 7 3.7737
2 2 2 7 , 2 11 3.47 3988.83 2.21765 × 10 7 3.8889
2 3 2 8 , 2 12 41.25 68 , 313.10 5.70258 × 10 8
Table 4. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when r = 0.5 τ h 2 > 3 7 .
Table 4. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when r = 0.5 τ h 2 > 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 2 TCT H 2 nd u x 2 ε H 2 nd u x 2 H 2 nd u x 2
2 1 2 4 , 2 8 0.03 7.52 3.72102 × 10 6 1.7361
1 2 5 , 2 9 0.13 64.38 2.14327 × 10 6 2.6720
2 2 6 , 2 10 0.59 533.53 8.02135 × 10 7 3.2757
2 2 2 7 , 2 11 3.83 5122.09 2.44877 × 10 7 3.6202
2 3 2 8 , 2 12 42.91 73 , 957.51 6.76426 × 10 8
Table 5. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when (172) and (173) are used and r = 0.5 τ h 2 3 7 .
Table 5. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when (172) and (173) are used and r = 0.5 τ h 2 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 1 TCT H 2 nd u x 1 ε H 2 nd u x 1 H 2 nd u x 1
2 6 2 4 , 2 13 0.03 216.25 3.93819 × 10 6 4.5863
2 5 2 5 , 2 14 0.09 1695.39 8.58690 × 10 7 4.6031
2 4 2 6 , 2 15 0.63 18945.40 1.86547 × 10 7 4.6131
2 3 2 7 , 2 16 3.67 218517.01 4.04385 × 10 8
Table 6. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when (174) and (175) are used and r = 0.5 τ h 2 3 7 .
Table 6. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when (174) and (175) are used and r = 0.5 τ h 2 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 2 TCT H 2 nd u x 2 ε H 2 nd u x 2 H 2 nd u x 2
2 6 2 4 , 2 13 0.03 251.27 3.37221 × 10 6 2.5722
2 5 2 5 , 2 14 0.13 2088.16 1.31103 × 10 6 4.3277
2 4 2 6 , 2 15 0.63 18945.40 3.02939 × 10 7 4.4163
2 3 2 7 , 2 16 3.85 234313.60 6.85956 × 10 8
Table 7. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when (172) and (173) are used and r = 0.5 τ h 2 > 3 7 .
Table 7. Computational time, maximum norm of the errors, and the order of convergence of v h , τ to the exact derivative v = u x 1 when (172) and (173) are used and r = 0.5 τ h 2 > 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 1 TCT H 2 nd u x 1 ε H 2 nd u x 1 H 2 nd u x 1
2 1 2 4 , 2 8 0.02 5.08 3.93866 × 10 6 4.5862
1 2 5 , 2 9 0.08 38.19 8.58815 × 10 7 4.6030
2 2 6 , 2 10 0.44 352.03 1.86579 × 10 7 4.4176
2 2 2 7 , 2 11 3.52 3994.16 4.22355 × 10 8
Table 8. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when (174) and (175) are used and r = 0.5 τ h 2 > 3 7 .
Table 8. Computational time, maximum norm of the errors, and the order of convergence of z h , τ to the exact derivative z = u x 2 when (174) and (175) are used and r = 0.5 τ h 2 > 3 7 .
r = 0.5 τ h 2 h , τ CT H 2 nd u x 2 TCT H 2 nd u x 2 ε H 2 nd u x 2 H 2 nd u x 2
2 1 2 4 , 2 8 0.02 5.89 3.67669 × 10 6 2.8278
1 2 5 , 2 9 0.11 45.67 1.30019 × 10 6 4.4268
2 2 6 , 2 10 0.50 414.27 2.93712 × 10 7 4.5165
2 2 2 7 , 2 11 3.72 4475.91 6.50300 × 10 8
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Buranay, S.C.; Matan, A.H.; Arshad, N. Two Stage Implicit Method on Hexagonal Grids for Approximating the First Derivatives of the Solution to the Heat Equation. Fractal Fract. 2021, 5, 19. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract5010019

AMA Style

Buranay SC, Matan AH, Arshad N. Two Stage Implicit Method on Hexagonal Grids for Approximating the First Derivatives of the Solution to the Heat Equation. Fractal and Fractional. 2021; 5(1):19. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract5010019

Chicago/Turabian Style

Buranay, Suzan Cival, Ahmed Hersi Matan, and Nouman Arshad. 2021. "Two Stage Implicit Method on Hexagonal Grids for Approximating the First Derivatives of the Solution to the Heat Equation" Fractal and Fractional 5, no. 1: 19. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract5010019

Article Metrics

Back to TopTop