A Set of New Stable, Explicit, Second Order Schemes for the Non-Stationary Heat Conduction Equation

: This paper introduces a set of new fully explicit numerical algorithms to solve the spatially discretized heat or diffusion equation. After discretizing the space and the time variables according to conventional ﬁnite difference methods, these new methods do not approximate the time derivatives by ﬁnite differences, but use a combined two-stage constant-neighbour approximation to decouple the ordinary differential equations and solve them analytically. In the ﬁnal expression for the new values of the variable, the time step size appears not in polynomial or rational, but in exponential form with negative coefﬁcients, which can guarantee stability. The two-stage scheme contains a free parameter p and we analytically prove that the convergence is second order in the time step size for all values of p and the algorithm is unconditionally stable if p is at least 0.5, not only for the linear heat equation, but for the nonlinear Fisher’s equation as well. We compare the performance of the new methods with analytical and numerical solutions. The results suggest that the new algorithms can be signiﬁcantly faster than the widely used explicit or implicit methods, particularly in the case of extremely large stiff systems.


Introduction
It is generally known that conductive heat transfer is described by a second-order linear parabolic partial differential equation (PDE), the so-called heat equation: or more generally where u = u(x, t) or, in the case of Equation (2) , u = u( r, t) is the temperature, α = k/(cρ) > 0 is the thermal diffusivity, q, k, c, and ρ are the intensity of heat sources (electromagnetic radiation, chemical reactions, radioactive decay, etc.), heat conductivity, specific heat and (mass) density, respectively. These equations and their generalizations such as the convection (or advection)-diffusion equation and reaction-diffusion equation are used to model the diffusion of various particles in chemical and biological systems [1] as well as in electronic devices [2]. Mathematically very similar equations are used to treat the flow of fluids through porous media [3], like moisture [4], ground water or crude oil in reservoirs. When diffusion or fluid flow is modelled, the u variable denotes the concentration or the pressure. We will examine one of the nonlinear reaction-diffusion equations, namely Fisher's equation [5], which originally was introduced to model how advantageous gene-variants spread in space and time.

The Description of the Methods
To solve Equation (1) numerically, we begin with the most typical starting step, as in the standard method of lines. We discretize the second space derivatives by the most widely used second order central difference formula ( [38], p. 48) First, we examine the one dimensional heat equation in a homogeneous medium. After defining an equidistant grid of N grid nodes we get the following ordinary differential equation (ODE) system: where Q i is the value of the q(x) source function at node i. This system of equations can be written in a briefer matrix-form: The diagonal elements of the matrix M are the sums of the non-diagonal elements of the appropriate rows with a negative sign: We introduce the characteristic time or time-constant of the cell i:

The One-Stage Constant-Neighbour (CN) Method
In this point we will recall the already published [36] method to solve the ODE system (4), which can be derived through the following steps:

1.
We make an approximation: When we calculate the new value of a variable u n+1 i , we neglect that other variables are also changing during the time step h = ∆t . This means that we consider u j a constant if j = i , so we have to solve uncoupled, linear ODEs: where the u i unknowns are considered as a function of time with initial values u n i , which are the temperatures at the beginning of the n-th time step. We introduced 2.
The ODE system (5) is a highly simplified approximation of the original PDE, and it is straightforward to solve it analytically: We use this solution at the end of the time step to obtain the new values of the u variable: For a standard one dimensional homogeneous system with an equidistant grid, we have the following one-stage first order method, which was called "constantneighbour method" or CN scheme.
The CN method algorithm while the actual values can be denoted as u CN i . Formula (7) can be written as: where the simple identity e ϕ = 1 + (e ϕ − 1) has been applied. Using the power series of the exponential function, we get: Inserting this back to the previous (8) formula and keeping only the first term yields: which is exactly the explicit Euler method, also called Forward Time Central Space (FTCS) scheme in the case of PDEs. One can see that regarding one time step, the one-stage CN method is identical to the explicit Euler method, but only up to first order. However, in spite of this apparent similarity, there is a fundamental difference between the methods: During step (1) we do not fix the actual variable u i to its initial value (we fix only the neighbours), but keep it as a variable. This means that at the end of this step we still have differential equations and not difference equations to be solved analytically, i.e., here time derivatives are not approximated by finite difference formulas as in Finite Difference Methods (FDM). In our CN method, the time step size h in the Euler/FTCS scheme is replaced by 1 − exp − 2αh ∆x 2 ∆x 2 2α and because of these exponential terms, Formulas (6)-(8) contain h up to infinite order, which is crucial for the unconditional stability.

The Two-Stage Constant-Neighbour Method
Now we show how to create a more accurate, second order two-stage "combined CN" algorithm. In Table 1, one can see the Butcher tableau of the well-known explicit generic second order Runge-Kutta (RK) method. In three cases these methods obtained distinct names: for p = 1, one has the explicit trapezoidal rule, which is also called Heun's method or improved Euler's method. • for for p = 2/3, it is called the Ralston-method, which provides a minimum bound on the truncation error for the second-order RK algorithms [39]. The first attempt could be to replace the time step h in these second order RK schemes 2α to stabilize it. We obtained that although this replacement indeed yields stable algorithms, they are still not second order. That is why we must return to the original logic by which we derived the first order CN method and combine it with the logic of the two-stage RK method above. At the first stage, we can calculate new "predictor" values of the variables, using (6) but with h 1 = ph time step: At the second stage, we can calculate the linear combination: and take a full time step using these combinations as the more accurate constant values u comb j of the neighbours towards which the u i unknown functions are tending with increasing time This means we have to replace a i in (6) to: For a standard one dimensional homogeneous equidistant grid, we suggest the following Algorithm 1: Stage 2. Take a full h time step with the CN method: The C 1 2 C , C 2 3 C , C1C versions are analogous to the (explicit) midpoint, Ralston and trapezoidal methods, respectively. In Section 3 we will analyse the general CpC algorithm and prove that it is second order in the time step size h. One might think that this whole procedure of constructing new methods using the already known RK schemes is trivial. To refute this opinion it is enough to note that following the logic explained above and using the Butcher-tableau of third order RK methods, the obtained methods are not third order, sometimes not even second order.

Extension to a General Grid
For a realistic system, the discretization must reflect the geometrical and material properties of the system; thus, we start from the more general Equation (2), where α, k, c and ρ are not considered uniform in space. Using Formula (3) for a one dimensional, equidistant grid, we have Now we change to cell variables, where the indices refer to whole cells: Here u i is the temperature of the cell i, is the heat capacity of the cell in units[J/K], ∆x is the length, V = A∆ is the volume while m is the mass of the cell. We introduce two other quantities, the heat source term Q of the cells: In the case of an irregular grid, the distances between the cell-centres are d ij = ∆x i + ∆x j /2 and the resistances can be Using these notations we have: This can be generalised easily to obtain the ODE system for a general (possibly unstructured) grid, which gives the time derivative of each temperature independently of any coordinate-system: Formally we can write this equation system into the same matrix-form as in (4). Now an off-diagonal m ij = 1/ R ij C i element of the M matrix can be nonzero only if the cells i and j are neighbours, i.e., there is direct heat conduction between them. Using this fact we can write the diagonal elements of the matrix as: After this generalization of the system matrix and time constant τ we can immediately use not only the CN algorithm (6) but the Formulas (9)-(12) as well for the CpC method.

The Properties of the Methods for the Heat Equation
The new CN and CpC methods can be used to solve general ODE systems as well as the special ODE systems obtained by the semidiscretization of PDEs. In the following pages, we will consider both aspects, but we certainly do not state that they can be recommended for any kind of ODE systems and PDEs. First, we collect those properties of the methods which are obvious and do not require to be proved:

•
They are explicit; we use matrices only to store data in a structured way, but do not perform operations with them and in the loop for the time steps, we use only vectors, not matrices. It also implies that the methods are easily parallelizable. • They are one-step methods. When the u n+1 i values are calculated, only the values u n j at the beginning of the current timestep are used. This implies that the methods are self-starting, the step size can be changed without any difficulty and the memory requirements are quite low. • They can be easily applied for any space dimensions, unstructured grid or inhomogeneous heat conduction medium provided that the heat capacity of the cells and the thermal resistance between the cells can be calculated.
At this point, we have to underline again that these new methods may not be considered as Finite Difference Methods. Here time derivatives are not approximated by finite difference formulas as in FDM. Because of the exponential terms, our methods might also seem to be similar to the so-called generalized Runge-Kutta Formula [40] or exponential integrators [41]. However, those methods use matrix exponentials, while we do not even need to use matrices during the calculation; thus, our methods are fundamentally different from those.

Convergence
In this section, we analyse the convergence properties of the methods as solvers for the spatially discretized systems. In Section 3.3 these will be examined regarding the methods as PDE solvers. Theorem 1. The new CpC method, defined by Equations (9)-(12), is a second order numerical method for arbitrary nonzero p ∈ R and for the general: linear ODE initial value problem, where M is an arbitrary nonsingular constant matrix, while Q and u(0) are arbitrary vectors.
Proof of Theorem 1. To help the reader follow the argument, we present the matrix form of the ODE system: The exact solution of the initial value problem (13) can be written as follows: Without the loss of generality we examine only u 1 at the first time step. It implies that we use the initial values u j (0) at the beginning of the time step. The 0 th and first order terms in the exact solution at t = h are the following The second order terms: When we apply the CpC method for the general initial value problem (13), at the first stage we obtain: u The result of the second stage in the CpC method is the following: The first term, up to second order is simple: To examine the second term we have to do the substitutions using (11) and (10): which can be written as: As the last term contains terms at least first order, we can omit the second order terms in the large bracket and obtain: Notice that the parameter p has been completely cancelled out. Now we can easily see that the zeroth and first order terms in (16) and (17) together are the same as in (14). Let us collect the second order terms: which is equivalent to (15), i.e., the CpC method is (at least) 2nd order.

Stability
We recall that when the stability of the ODE solvers are examined, typically the following scalar equation is used Absolute stability means that arbitrarily large time steps yield a bounded solution for this equation [42] and [43] (p. 13). It is well-known that the stability function of any explicit Runge-Kutta method is a polynomial, so they can never be A-stable, which limits their use for solving PDEs ( [43], p. 24). Now, for this scalar equation, i.e., in the case of one cell, our methods can be defined only in one meaningful way: To provide the analytical solution. This means that in this sense, they are trivially A-stable. We therefore examine the stability when u is a vector, as it is the case for spatially discretized PDEs. Unlike Runge-Kutta methods, Formulas (6), (9) and (12) use an exponential expression of the time step size h to calculate the new values of the variable u. At the first stage, Equation (9) or Stage 1 in Algorithm 1 the temperature u i of each cell tends to a i τ i with increasing h, which is the weighted average of the temperatures of the neighbouring cells at the beginning of the time step. At the second stage, u i tends to the combination of the beginning temperatures of the neighbours and the neighbours of the neighbours with increasing h. We will formulate this statement in Theorem 2, but first we evoke the following simple lemma, the associativity of convex combinations ( [44], p. 28) on which we build during the proof of the theorem. Proof of Theorem 2. We will use the facts that τ i = −1/m ii and m ij,j =i are non-negative quantities; thus, 0 < e − h/τ i = e m ii h ≤ 1 holds because of physical reasons, namely the second law of thermodynamics. Again, without the loss of generality, we examine only u 1 at the first time step.
(A) Let us examine the first stage: The term a 1 τ 1 is the following: Let us calculate the coefficients of the u j (0) initial values: One can see that they are nonnegative and their sum is one, which implies that a 1 τ 1 is a convex combination u j (0). Now it is clear that the coefficients e The square bracket is a convex combination if 0 ≤ 1 2p ≤ 1 and 0 ≤ 1 − is a convex combination of the initial values.
Theorem 2 implies that any error produced during the calculation will be distributed quickly and cannot be amplified; thus, the methods are stable for the heat equation (and as we will show later, for Fisher's equation). Moreover, in the homogeneous (Q = 0) case, Theorem 2 guarantees that the solution follows the maximum and minimum principles ( [43], p. 87), i.e., the extreme values of u occur among the initial or the prescribed boundary values. This means that starting from positive initial values and assuming no heat sinks are present, the variables always remain positive. This property is especially advantageous when Equations (1) and (2) are used to model physical processes, where the temperature measured in Kelvin (or analogous variables such as concentrations of chemical species) cannot be negative.

Analysis of the Errors
We analytically calculated the first two terms of the error of the CpC scheme for a system of 2 identical cells. At the end of the first time step, the error for the first cell is the following: We note that we obtained similar expressions for larger systems and for the higher order error terms as well, but of course with more terms and with larger exponents. We can draw two conclusions.
First, for not too large time steps, the error is smaller for smaller parameter values p, so p should be chosen as small as possible. Unfortunately, this advice contradicts the requirement of stability, which is not guaranteed for p < 0.5. We will return to this question during the discussion of the numerical experiments.
For the second consequence we recall that τ ∼ ∆x 2 , which means that for small ∆x , the main error terms contain powers of h ∆x 2 . We have to conclude that the methods are only conditionally consistent in the sense that when the size of the spatial cells ∆x tends to zero, the time step size h also has to go to zero faster than the second power of ∆x . If ∆x → 0 while h/(∆x) 2 is kept constant, then the global errors tend to some constant values, and if one decreases h faster, the errors tend to zero. However, if h is decreased slower, the errors are growing, but not without limits: They approach their maximum possible values which are restricted by the above mentioned maximum and minimum principles. It implies that considering as solvers for the original PDE, the methods are conditionally convergent. However, as we have proved above, they are unconditionally convergent for any fixed spatial grid, i.e., the solution converges to the exact solution of the ODE system obtained by the discretization of space variables when the time step size h tends to zero. It is worth mentioning that this second consequence holds for the CN method as well.

Comparison with Exact Results
We note that all simulations are performed using the MATLAB R2019b software on a desktop computer with an Intel Core i3-8100 as CPU.

The Heat Equation
We solve PDE (1) on the interval x ∈ [0, 1] with the following heat source term: q(x, t) = q 0 sin(πx), where q 0 = 2 and α = 1 . The initial condition is also a sine function but with smaller wavelength: u(x, t = 0) = sin(3πx) and we use zero Dirichlet boundary conditions: It is easy to check that the analytical solution of this problem is: We discretize x by setting: x j = j∆x , j = 0, ..., 100 , ∆x = 0.01.
The solutions are examined and compared at final time t fin = 0.2 . We define the (global) error as the maximum of the absolute value of the difference between the exact temperature u exact i and the temperature u num i obtained by our numerical methods at the end of the examined time interval: In Figure 1, one can see these global errors as a function of time step size h. One can see that the error tends to a small nonzero value in all cases. This non-vanishing error is caused by the discretization of the space variable. More precisely, it is the truncation errors ( [43], p. 7) of the central difference Formula (3) accumulated during the subsequent time steps. In Figure 2, we illustrate the behaviour of the maximum errors for different space steps and for p = 0.5. The minimum of these functions appears in the form of a cusp in each case. The reason for this is that the largest error terms (19) and (21), coming from the time-integration method and the spatial discretization, respectively, are of opposite signs and cancel each other at these special, optimal time step size and grid spacing combinations. Unfortunately, we observed this behaviour only when the source term Q strongly dominates the effect of the (non-uniform) initial conditions. Moreover, the location of the cusp, i.e., the optimal h depends on so many factors that this favourable property can hardly be exploited in practical problems.

The Nonlinear Fisher's Equation
We consider the following equation with α = 1 and β > 0 , subject to the following initial condition: The analytical solution of this problem is the following travelling wave function [45,46] The appropriate Dirichlet boundary conditions are prescribed at the two ends of the interval: The values of the exact solution obviously lie in the unit interval, u(x, t) ∈ [0, 1] , x, t ∈ R. We take the nonlinear term into account in a "semi-implicit" way: where u n+1,pred i is the result of the first and the second stage using the same formulas as in the previous section, i.e., with taking into account the effect of the diffusion term only. Note that now the new value u n+1 i appeared on the right hand side of the equation as well. However, this equation can be rearranged into a fully explicit form: This operation is performed in a separate third stage after Stage 2 (a separate second stage after the first one in case of the CN method), where a loop is going through all the nodes. We have done several tests with this semi-implicit treatment and obtained that the methods behave similarly well as in the linear case. The errors as a function of the time step size h are presented in Figure 3 for β = 2.5, x 0 = 0, x fin = 4 , t fin = 2 , and ∆x = 0.01 . (23), it is obvious that the new value u n+1 i cannot be negative, and also cannot exceed 1; thus, 0 ≤ u n+1 i ≤ 1.

Corollary 1.
In the case of Fisher's equation (22) if the initial values are in the interval u 0 i ∈ [0, 1], then the values of u remain in this interval for arbitrary β and time step size h, for the whole calculation. This is a special case of the maximum and minimum principles, which implies that the CN and the CpC methods supplemented with (23) is unconditionally stable for Fisher's equation.
Proof. Theorem 2 and Remark 1 immediately imply the statement.
However, we do not state that we have found the optimal algorithm to solve Fisher's equation. The purpose of this subsection is to demonstrate that these methods can be successfully applied for nonlinear equations as well. The systematic investigation of the application possibilities of the methods for nonlinear equations is clearly out of the scope of this paper and planned to be published in a series of papers in the future.

Comparison with Numerical Results
In this section rectangle-shaped lattices are examined with N = N x × N z cells, where the edge of the system is thermally isolated regarding conductive type heat transfer, which means closed (zero Neumann) boundary conditions. To help the reader to imagine the system, we present the arrangement of the variables in Figure 4. We emphasize that the size and the shape of the cells are not necessarily identical. To obtain a reference solution, we used the implicit ode15s solver of MATLAB, which is a variable-step, variable-order solver based on the numerical differentiation formulas (NDFs) of orders 1 to 5, where the letter s indicates that the codes were suggested to use in the case of stiff systems. During this calculation, we applied stringent error tolerance ('RelTol' and 'AbsTol' were both 10 −16 ). We compare the performance of the new CN and CpC methods with professionally coded and widely used MATLAB routines. In addition to ode15s, the following solvers are employed, with loose error tolerance to increase the speed: • ode45, the explicit Runge-Kutta-Dormand-Prince formula of order 4(5). • ode23, the explicit Runge-Kutta-Bogacki-Shampine method of order 2(3).
• ode23t, an implementation of the trapezoidal rule using a "free" interpolant . • ode23s, a modified Rosenbrock formula of order 2. Besides the above listed MATLAB solvers, we have coded the previously mentioned unconditionally positive finite-difference (UPFD) method of Chen-Charpentier et al. [27] for comparison purposes because of its simplicity. Applying that scheme to (1) we obtain: This formula can be rewritten in the following explicit form: which can be generalized in the manner we explained in Section 2.3, so the UPFD formula we use: Here we define the (global) error in the same way as in (20) but now we have a reference solution instead of the exact analytical solution:

First Case: A Very Stiff System
We set N x = 100, N z = 40; thus, the number of cells is N = 4000. Different random values are given to the heat capacities and to the thermal resistances: where rand is a random number generated by the MATLAB uniformly in the (0,1) interval for each quantity. This means that the capacities (the resistances) follow a log-uniform distribution between 0.01 and 100 (between 0.001 and 1000). The initial temperatures have a uniform distribution between 0 and 1, i.e., u i (0) = rand, while the source-terms have a uniform distribution between −0.5 and 0.5 . Starting from t 0 = 0s we solve this system for the temperatures at final time t f in = 1s.
As the system is thermally isolated, the system matrix has a zero eigenvalue, all other eigenvalues are negative. If we denote the eigenvalues with the largest (smallest) absolute value with λ MAX (λ MIN ), then the stiffness ratio can be calculated as For the explicit Euler method, the maximum possible time step size is: above this threshold instability necessarily appears. In Figure 5, we present the errors for the one-stage CN and the two-stage CpC methods in the case of different p parameters as a function of the step size h. From the figure one can see that for small h, the global error decreases with the first power of the step size in the case of the Constant-neighbour methods and with the second power for the CpC methods, which confirms that the methods are first and second order, respectively. It can be seen that for p = 1/3 the method is unstable in the region of large step sizes, but in the region of small h-s it is slightly more accurate than for larger p-s. For p = 1/2 the error-function has a relative elevation around h = 10 −2 ...10 −3 , because this parameter value is at the border of the stability region. We also examine the global errors as a function of the total running times. The results are presented in Figure 6. For smaller step sizes the running times of several runs have been averaged to avoid random fluctuations due to very short running times. We note that the explicit solvers ode23 and ode45 can hardly provide any meaningful results, because for large error tolerance they diverge, while tightening the tolerance yields an error message "array exceeds maximum array size".
In Table 2 we summarize some results obtained by MATLAB routines ode15s, ode23s, ode23, ode23t, the UPFD scheme, and our methods C 1/2 C and C 2/3 C. One can see that, if we set the same moderate precision requirements, our methods are much faster than the conventional explicit or implicit methods, even without adaptive stepsize control, parallelization, or vectorization. The UPFD method is also very fast, but its accuracy is far from being optimal.  Table 2. Performance of different solvers for the system consisting of N = 4000 cells described in Section 5.1. CpC means the new two-stage method. The absolute and the relative tolerance 'RelTol' and 'AbsTol' have been set to the same values denoted by 'tol'. The error is defined by (25).

Method
Running

Second Case: A Highly Anisotropic System
Now N x = 100, N z = 100; thus, N = 10,000. The capacities and the thermal resistances followed a log-uniform distribution again: Note that, on average, there is four orders of magnitude difference between the resistances in the x and in the z direction. The initial temperatures are u i (0) = 2rand − 1 while the source-terms are Q i = rand − 0.5 . The task is to solve this initial value problem for the temperatures at t f in = 1s starting from t 0 = 0 s. The stiffness ratio is rather high, 2.15 × 10 9 . The maximum possible time step size for the explicit Euler method is h EE MAX = 9.6 × 10 −5 s, a few times larger than in the first case.
In Figure 7, we present the error for the CN and the CpC methods as a function of the step size h. We also present the global errors as a function of the total running times in Figure 8. In this case, due to the larger h EE MAX , the explicit solvers ode23 and ode45 are able to produce some results, albeit still in a relatively narrow region of tolerance values 'RelTol' and 'AbsTol'. One can see that our methods are much faster than the conventional explicit or implicit methods. One might think that the two-stage C½C method produces too large errors for medium step sizes, e.g., for h = 5 × 10 −4 . However, in this case the average absolute error is only 0.002 and from Figure 9 one can see that the graph of the produced solutionfunction is almost indistinguishable from the reference curve. This means that our methods can produce acceptable results three orders of magnitude faster than the well-established MATLAB routines. In Table 3, we summarize some results obtained by MATLAB solvers ode23, ode45, ode15s, ode23s, ode23t, the UPFD scheme, and our two-stage CpC methods.  We emphasize that the N = 10,000 number of cells is still much smaller than in many real-world applications. Since for larger systems, implicit methods have a more serious disadvantage, we think that our methods are even more competitive in those cases. We note that the reason we have not used a larger number of cells for testing purposes is that our computer cannot solve them using the MATLAB routines because of the too large memory requirements.

Conclusions
The purpose of this paper is to introduce and examine a set of new explicit one-step numerical schemes. We propose these methods to solve the ODE system which is obtained after spatially discretizing diffusion term. We prove that one-stage "constant-neighbour" (CN) method is first order, while the two-stage "CpC" methods are second order for the heat equation in the time step size h for arbitrary parameter values p. For p ≥ 0.5 the variables at the end of the time-step are the convex combinations of the values at the beginning of the time step for the heat equation in the absence of external sources; thus, the methods are stable for arbitrarily large values of the time step size. This holds for the nonlinear Fisher's equation as well, if one treats the nonlinear term as we suggest. These statements have been proved analytically and confirmed by several numerical experiments.
The schemes have been verified by comparisons to analytical solutions of the inhomogeneous heat equation and Fisher's equation. The performance of the methods has been tested in the case of two discretized two-dimensional systems with highly inhomogeneous random parameters and initial conditions. The obtained data show that if quick results are required, the proposed schemes have a significant advantage in speed, even without adaptive stepsize control, parallelization or vectorization. We must add that they are easy to implement and parallelize and they can be applied regardless of space dimensions and grid irregularity. On the one hand, increasing number of cells makes our methods more competitive against the widely used implicit schemes. On the other hand, if the range of the eigenvalues of the system matrix is wide (which is in connection with stiffness), the new methods have a significant advantage over the traditional explicit schemes, which are only conditionally stable. However, when highly accurate results are required, first and second order methods are not the best choice; thus, we have already started to search for higher-order versions of the new methods.