NUMERICAL ANALYSIS OF A TIME RELAXATION FINITE DIFFERENCE METHOD FOR THE HEAT EQUATION

. In this study, we first consider the time-relaxation model, which consists of adding the term κ ( u − u ) to the heat equation. Then, an explicit discretization scheme for the model is introduced to find the finite difference solutions. We first obtain the solutions by using the scheme and then investigate the method’s consistency, stability, and convergence properties. We prove that the method is consistent and unconditionally stable for any given value of r and appropriate values of κ and δ . As a result, the method obtained by adding the time relaxation term to the first-order finite-difference explicit method behaves like the second-order implicit method. Finally, we apply the method to some test examples.


Introduction
The heat equation has a fundamental importance in various scientific fields.Heat is a form of energy that exists in any material.For example, the temperature in an object changes with time and the position within the object.Furthermore, this equation can be applied to solve the heat flow related to science and engineering.The numerical and analytical methods can be used to solve the heat equation problem in science or engineering fields.The finite difference method is one of several techniques for numerical solutions to boundary value problems.Morton and Mayers [16], and Cooper [5] provide a modern introduction to partial differential equations theory by considering the development of finite difference methods and 1078 Ö. İLHAN, O. R. ISIK, S. BOZKURT numerical methods in more detail.Fletcher [7] described the method to implement finite differences to solve boundary value problems.
The time relaxation (based on spatial averaging) regularizations considered by Rosenau [19], Schochet and Tadmor [20], Stolz et.al. [22] are often referred to as "secondary regularization" or "time relaxation" to regularize the flow.Stolz and Adams [1] studied the approximate deconvolution model (ADM) for the large-eddy simulation of incompressible flows and applied it to turbulent channel flow.Layton and Neda [13] developed a time relaxation regularization of flow problems suggested by Stolz and Adams.Ervin et al. [6] studied the numerical errors in finite elements discretization of a time relaxation model of fluid motion.They showed that the point of the relaxation term is to approach the unresolved fluctuations in a computational simulation to zero exponentially fast by an appropriate and depending on the problem choice of its coefficient; thus, they concluded that this relaxation term is intermediate between a tunable numerical stabilization and a continuum modeling term.Neda [17] investigated a high-order family of time relaxation models based on approximate deconvolution.This problem was also considered in [14] to regularize the flow.They studied the analogous approach based on time scales, time filtering, and damping of under-resolved temporal features and investigated theoretical and practical aspects of temporally damped fluid-flow simulations.Isik et al. [12] analysed the NSE time relaxation model which is obtained adding term "κ (u − u)" where u denotes the time filter of u introduced by Pruet et al [18].An advantage of adding the term time-relaxation is that it provides a faster approach to steady-state solutions [12], [10] [11].In some problems obtained by adding this term, the solutions can be achieved in some cases where the numerical solutions can not be obtained [12].Cibik et al. [4] introduced the backward Euler method obtained by adding the time relaxation term to MHD flow.They proved that the method improves the accuracy of the solution without a significant change in the complexity of the system.
We know that leapfrog discretization time filter methods are applied to geophysical fluid dynamics.While these methods decrease spurious oscillations to improve estimations, they reduce numerical accuracy and excessively dampen the physical model [2].Williams proposed a successfully tuned model which reduces undesired numerical damping of [2] with higher-order accuracy [23].This model is also studied by others see; [15], [24].Otherwise, it is often preferable to use the backward Euler method for a steady-state problem in practice.This method is stable, but the solutions are inefficient and time accurate transient [8].Guzel and Layton stabilized the backward Euler discretization using time filter for the classical numerical ODE theory by improving solutions [9].
In this paper, we investigate the finite difference solution of the model, which is obtained by adding the time relaxation term "κδu t " to the heat equation, where, u is the differential filter defined as [18], Here, κ and δ are positive real constants.It is well known that the finite difference solution of Eqn (1) obtained by the explicit method is stable for r ≤ 1 2 whereas the implicit method is unconditionally stable for κ = 0. Also, the explicit method is consistent and the convergence order for this method is O(k + h 2 ) for κ = 0. To obtain the finite difference solutions of eqn.(1) we will develop an algorithm by using the explicit method.We analyze the consistency and stability of the finite difference solution obtained by this algorithm.Thus, the proposed method is consistent on an unbounded domain subset of κδ-plane.For any given value of r, the method will be stable for some particular values of κ and δ.Although this algorithm seems structurally explicit, the behavior of stability results is similar to the implicit method.We applied the method to some numerical examples.We found that the numerical results are consistent with the theoretical results.For each example, we show that the convergence order for the method is 1 and 2, and the method is stable for any given r.All computations are done using MATLAB R2020b.
The paper is organized as follows.An algorithm that depends on the finite difference method is introduced for Eqn (1) in Section 2. Then consistency and stability of the finite difference method obtained by the algorithm are investigated, and some properties are given in the same section.Some numerical examples are given to verify the theoretical results in Section 3. It is seen that the numerical results obtained by numerical examples are consistent with the theoretical results in this section.While the solutions blow-up for κ = 0 and r = 5, convergence results are obtained for appropriate values of κ and δ.Finally, the conclusions of this study are given in Section 4.

Finite Difference Approximation
The finite difference method is a numerical method in which the solution of a differential problem given on an interval at a point is approached with finite differences.The methods generally generate the solutions that are either as accurate as the data warrant or as accurate as is necessary for the technical purpose for which the solutions are required.One of the finite difference approximations to the "Heat equation", is where U is the exact solution of the heat equation and u i,j+1 is the unknown temperature at the (i, j + 1)th mesh point in terms of known "temperatures" along the jth time-row.This method is referred as explicit method.The consistency and stability results for explicit method can be found in any book involves finite difference methods, i.e., [21].
Theorem 1. (Lax's Equivalence Theorem) Given a properly posed linear initialvalue problem and a linear finite-difference approximation to it that satisfies the consistency condition, then stability is the necessary and sufficient condition for convergence [21].
We consider an explicit discretization of model ( 1): where κ ≥ 0, δ > 0 and u i,j := u(x i , t j ) is the known value of temperature at the (i, j)th mesh point.Note that u i,j+1 and u i,j+1 can be written as and respectively.Thus, by substituting Eqn (7) into Eqn (4), we obtain the following explicit scheme Now, we will show the von-Neumann stability of Eqn (8).
Lemma 1.Let u i,j be the finite difference solution of model (1) obtained by Eqn (8).Then, the method is unconditionally stable on an unbounded region in the κδ−plane.
Proof.Substituting Eqn (6) into Eqn (8) yields Then, putting u i,j = e iβph ξ q into the finite difference scheme yields If we simplify "e iβph ξ q " term from both sides, we get Multiplying both sides with "ξ" gives Then, we can find two solutions of Eqn (9): Let us find the stability condition |ξ| ≤ 1 for ξ 1 .
On the other hands, for the second root, we can use the same process: Let us find κ and δ as k + 2δ − 2rδ − kr + kr cos θ + kκδ + 2rδ cos θ ≥ 0. Since 2rδ + kr − kr cos θ − 2rδ cos θ ≤ 4rδ + 2kr and yields the stability of the scheme.It can be said that the increased δ will give stability for kκ − 4r > 0. □ 2.1.Consistency.In this subsection, we deal with the consistency of the method for various κ and δ.We will prove that the convergence order of the method is . Specifically, we proved that for a constant δ and under the condition κδ ≤ 2 k , the method is consistent and the convergence order is 1.Proposition 1.Let k, h, κ and δ be fixed positive real numbers.If u = u(x, t) is smooth enough, then the method is consistent and the convergence order of the method is O(δ k 2 + k 2 + κδ k 2 − 1).Proof.Firstly, when u(x, t) is expanded to Taylor polynomials at (x, t) = (x i , t j ) point and then these polynomials are substituted into Eqn (4), we can get following, If necessary simplifications are done, we get Hence, T (x i , t j ) → 0 is satisfied for h, k → 0, so, the method is consistent with the convergence order O(δ k 2 + k 2 + κδ k 2 − 1) in time.As a special case, we will prove the consistency for a constant δ and under the condition κδ ≤ 2 k .To show the consistency, let us show the following inequality, Simplifying Eqn (11), we can get the following inequality Finally, we can obtain the following inequality which satisfies Eqn (11).It can be easily seen that lim (h,k)→(0,0) T (x j , t j ) = 0 is satisfied for a constant δ and under the condition κδ ≤ 2 k .Therefore, for a constant δ and under the condition κδ ≤ 2 k , the method is consistent and the convergence order for the method is 1. □ Note that it can come to the question to mind: Is it possible to get the convergence order from 1 to 2 in time?A way of doing this can order 2 be found in time for If necessary operations are performed here, we obtain By applying the stability condition ( 10) and ( 14), we get following inequality: It means that the method is order of 2 for values of κ in this interval.In addition, for any value of r and values of κ and δ satisfying the (15) condition under the stability condition, the method is stable and the convergence order is 2.
Corollary 1.If the stability condition and the (15) condition are satisfied for a given r value, a method that is convergent to the second order is obtained.

Numerical examples
In this section, we give some numerical examples to illustrate the method.We use MATLAB R2020b program for all computation.For the term u i,j+1 , the following relation can be used for each time step, so there is only one unknown term in time.
The results are given in Tables 1, 2, 3 and 4 and Figure 1.
Otherwise, the solutions are convergent with the order of 1. Figure 1 shows the κδ-stability domain which is unbounded subset of R 2 .The solutions are stable in this domain.In other words, for r = 6 and any chosen κδ-pairs which are in this region, the solutions do not blow-up.
It is seen in Table 3 that the solutions are convergent with the order of 2 in time for the δ = 0.0001 and selected κ values.Although r increases, the solutions are convergent with the order of 2.
It is shown in Table 4 that for different κ values the solutions are convergent with the order of 1 in time.These κ values are satisfied to Eqn. 13.
Example 2. Consider the heat equation with boundary condition for 0 ≤ x ≤ 2π and in time interval 0 ≤ t ≤ T = 0.5.
The exact solution is u(x, t) = e −t cos(x).
The results are given in Tables 5-10 and Figures 2 and 3.    2 and 3 show the κδ-stability domains which are unbounded subsets of R 2 .The solutions are stable in these domains.Namely, for r = 6.4038 and r = 50.6606and any chosen κδ-pairs which are in these domains, the solutions do not blow-up.
It is seen in Table 8 that the solutions are convergent with the order of 2 in time for the δ = 10000 and selected κ values.Although r increases, the solutions are convergent with the order of 2. Table 9 shows that for a constant value of h the solutions are almost convergent with the order of 2 in time.It is shown in Table 10 that for different κ values the solutions are convergent with the order of 1 in time.These κ values are satisfied to Eqn. 13.

Conclusion
In this study, we present a model which is obtained by adding the time-relaxation term κ (u − u) to the heat equation.We develop an algorithm using the explicit method to find the solutions of this model.We analyze the consistency, stability and convergence properties of the solution.We find that the method is consistent on an unbounded domain subset of κδ-plane.Moreover, we see that for any given value of r, the method is stable for some particular values of κ and δ which are selected from an unbounded consistency region which is the subset of R 2 .The algorithm seems structurally explicit.Besides, for any value of r, the behavior of the stability results shows that this algorithm is similar to the implicit method, but this method is still explicit.On the other hand, we obtain the convergence order of convergence as 1 when stability condition holds.If the condition 15 is satisfied in addition to the stability condition, the convergence order of the method is increased from 1 to 2, which means that the convergence is accelerated.Consequently, the presented model is an efficient model that stabilizes an unstable method by specifying unbounded κ − δ region.Moreover adding the term time-relaxation to the heat equation in this presented model expands the stability range for the explicit method.We give two examples to validate theoretical results and how the method works.We observe that the blow-up solutions which is obtained by κ = 0 are stabilized by selecting suitable pair of κ and δ as given in examples.This shows that for any value of r, the convergence order for the method is 1 and 2 for the appropriate values of κ and δ.As a result, the numerical results obtained by using the algorithm are consistent with the theoretical results.As further works, the time-relaxation term will be added to elliptic or hyperbolic equations.

Table 4 .
Numerical results of Example 1, for 5850 ≤ κ ≤ 3881.2 in the domain satisfying the stability condition, δ = 0.01 It can be seen from Table 1 that, for

Table 5 .
Numerical results of Example 2, for 11599 ≤ κ ≤ 664800 in the domain satisfying the stability condition, δ = 0.01 and r =

Table 6 .
Numerical results of Example 2, for 10439 ≤ κ ≤ 167040 in the domain satisfying the stability condition, 0.1 ≤ δ and r =

Table 10 .
Numerical results of Example 2, for κ = 4012.4 in the domain satisfying the stability condition, δ = 10000