New Stable, Explicit, Shifted-Hopscotch Algorithms for the Heat Equation

Abstract

Our goal was to find more effective numerical algorithms to solve the heat or diffusion equation. We created new five-stage algorithms by shifting the time of the odd cells in the well-known odd-even hopscotch algorithm by a half time step and applied different formulas in different stages. First, we tested 10⁵ = 100,000 different algorithm combinations in case of small systems with random parameters, and then examined the competitiveness of the best algorithms by testing them in case of large systems against popular solvers. These tests helped us find the top five combinations, and showed that these new methods are, indeed, effective since quite accurate and reliable results were obtained in a very short time. After this, we verified these five methods by reproducing a recently found non-conventional analytical solution of the heat equation, then we demonstrated that the methods worked for nonlinear problems by solving Fisher’s equation. We analytically proved that the methods had second-order accuracy, and also showed that one of the five methods was positivity preserving and the others also had good stability properties.

1. Introduction

Heat and mass transport has a central importance in many scientific and engineering applications. One of the most fundamental ways for these phenomena is the diffusion of heat or particles, described, in the simplest case, by the following linear parabolic partial differential equation (PDE), the so-called heat or diffusion equation. If the medium in which the diffusion takes place is not homogeneous, one can use a more general form.

(1)$\frac{\partial u}{\partial t}=\alpha {\nabla }^{2}u$

(2)$c\rho \frac{\partial u}{\partial t}= \nabla \left(k\nabla u\right)$

where, in case of diffusion of particles, $u=u\left(\overrightarrow{r},t\right)$ denotes the concentration and α is the diffusion coefficient. In the case of conductive heat transfer, u is the temperature and $\alpha =k/(c\rho )$ is the thermal diffusivity, while $k=k\left(\overrightarrow{r},t\right)$, $c=c\left(\overrightarrow{r},t\right)$, and $\rho =\rho \left(\overrightarrow{r},t\right)$ refer to the following nonnegative quantities: heat conductivity, specific heat, and mass density, respectively. These equations and their generalizations, such as the diffusion-convection-reaction equation [1], can describe the diffusive transport of particles in very different physical, chemical, and biological systems, such as in semiconductor devices [2] and the human brain [3]. Furthermore, mathematically similar PDEs can be applied to simulate damped material flow through porous media, such as moisture [4], ground water, crude oil, or natural gas in reservoirs below the Earth’s surface [5]. In this work, we examined only one nonlinear reaction-diffusion equation, the so-called Fisher’s equation [6], apart from the original, linear diffusion equation.

The simple diffusion Equation (1) or more sophisticated equations containing a diffusion term have more and more analytical solutions [7,8,9,10,11,12]. Unfortunately, most of these solutions consider the parameters like the diffusion coefficient or the heat conductivity as constants, which do not depend on space, time, or the dependent variable u. There are counter-examples, such as the work of Zoppou and Knight, in which “analytical solutions have been derived for the two- and three-dimensional advection-diffusion equation with a particular form of the spatially variable coefficients” [13]. However, for space-dependent coefficients in general, one needs to use numerical calculations; therefore, finding effective numerical methods are still important. This is especially true in the nontrivial case of systems where the physical properties are extremely different at different points even in the vicinity of each other [5,14] (p. 15). This often implies that the eigenvalues of the problem have a range of several orders of magnitude; therefore, the problem can be severely stiff. When the PDE is spatially discretized, one obtains a system of ordinary differential equations (ODEs). If the number of variables is very large (which is the case almost always in three space dimensions), the numerical solution of these systems still forms a challenge either by explicit or by implicit methods. It is well known that the traditional explicit methods (such as the Runge–Kutta types) are only conditionally stable, and very small (sometimes unacceptably small) time step sizes have to be applied. Implicit methods are typically unconditionally stable (there are exceptions, e.g., the higher order backward difference methods), but a system of algebraic equations must be solved at each time-step. These calculations can be still very time consuming (note that the matrix is generally not tridiagonal) and, therefore, immense efforts have been made to develop several sophisticated tricks and modifications of the implicit methods (see, e.g., [15] and the references therein). Currently, the implicit methods with these extensions are usually proposed to solve these kinds of problems [16,17,18]. The main obstacle to further progress in the feasibility of large-scale simulations is that the parallelization of these implicit methods are nontrivial, albeit, some progress has been achieved [19,20]. However, there is a tendency towards increasing parallelism in high-performance computing [21,22] due to the halt of the formerly rapid increase of the CPU clock frequencies in the last decades.

Keeping in mind these facts, we are working on the development of novel, easily parallelizable, explicit, and unconditionally stable methods. One of the most important example of these is the two-stage odd-even hopscotch (OEH) algorithm, which was introduced by Gordon [23], then reformulated and analyzed by Gourlay [24,25,26] (see also [27]). This method has been modified to increase its reliability and accuracy, but, as far as we know, always in the direction of a larger extent of implicitness. More concretely, a hierarchy of algorithms was constructed, from the fully explicit OEH through the line hopscotch to the alternating direction implicit (ADI) hopscotch, with increasing accuracy at the expense of programming and running time [25]. These methods were soon applied for the calculation of the heat flow in a thermal print head by Morris and Nicoll. They found that, in case of isotropic media, the OEH method is indeed faster than its more implicit versions with the same accuracy, but in anisotropic cases, the OEH and the line hopscotch give very inaccurate results and the ADI hopscotch is necessary to give a meaningful solution [28]. Later, the OEH method was applied to many problems, e.g., the incompressible Navier–Stokes Equations [29], the Frank–Kamenetskii [30] and the Gray–Scott reaction-diffusion equations [31], and even to the nonlinear Dirac equation [32]. Goede and Boonkkamp applied a vectorized implementation of the OEH scheme to the two-dimensional system of Burgers’ equations and found that the vectorization seriously increased the speed and the obtained solver was very powerful [33]. Recently, Maritim et al. constructed hybrid algorithms containing the hopscotch, the Crank–Nicolson, the Du Fort Frankel, etc. schemes for the two-dimensional system of Burgers’ equations, and concluded that their implicit algorithms are stable and accurate [34,35]. Surprisingly, all modifications of the OEH algorithm increased the level of implicitness and we found no traces of trying to construct similar schemes while keeping the fully explicit property. In our previous paper series [36,37,38], we constructed three new hopscotch combinations by choosing other formulas instead of the original explicit and implicit Euler schemes. The tests showed [36] that, for stiff systems, the original OEH method can produce a so huge inaccuracy for large time step sizes that the relative errors can reach 10⁴, which, in fact, can be more dangerous than instability if not noticed by an inexperienced user. We found that two of the three new combinations behaved much better. In this paper, we extended our research by modifying the underlying space and time structure as well.

The organization of the paper is as follows. In Section 2 we review the OEH structure very briefly and introduce the new algorithms, first for the simplest, one dimensional, equidistant mesh, then, for a general, arbitrary mesh as well. In Section 3, first, we report very concisely the results of the numerical tests for the first assessment of the 10⁵ methods. Then, in Section 3.3 and Section 3.4, we present two numerical experiments in two space-dimensional stiff systems consisting of a large number of cells. Based on the results, we selected the top five combinations with the most valuable properties. In Section 3.5 we verify these five methods by comparing the numerical results produced by them to analytical solutions of the heat equation in one space dimension. The behavior of the algorithms is examined in the case of Fisher’s equation in Section 3.6. Finally, we analyze the convergence and stability properties of these methods. Section 4 will be about the conclusions and our future research goals.

2. The New Methods

First, we recall the original OEH method for the simple one dimensional equation form of Equation (1). The space and the time variable is discretized by setting nodes according to the most standard rules: $x_i=i\Delta x , i=0,\dots ,N-1$, $t_n=nh , n=0,\dots ,T-1$. To introduce the OEH algorithm, we defined a so-called bipartite grid, in which the set of nodes is divided into two similar subsets A and B (with odd and even space index) such that all nearest neighbors of nodes in subset A belong to subset B and vice versa, like on a checkerboard. This space structure is in an intimate connection with the time discretization in the following way. At the first stage of the first time step, the new values of u were calculated only at the points of A, using the latest available values of the neighbors, i.e., the values at the beginning of the time step (see green arrows in Figure 1a). At the second stage, the remaining node values, which belong to subset B, were calculated, using again the latest values at subset A, which were now valid at the end of the time step and already calculated at stage one (red arrows). At the next time step, the roles of subsets A and B were interchanged, etc., as one can see in Figure 1a.

The original OEH method applied the well-known explicit Euler formula at the first stage and the implicit Euler formula at the second stage. However, when the second stage calculations began, the new values of the neighbors $u_{i-1}^{n+1}$ and $u_{i+1}^{n+1}$ were known, which made the implicit formula effectively explicit. We note, again, that earlier we created new OEH methods by applying different formulas instead of the abovementioned. Now, we modified not only the formulas, but the underlying space and time structure as well, as follows.

The calculation started with taking a half-sized time step for the odd nodes (subset A) using the already calculated $u_i^n$ values. Then, a full time step was made for the even nodes (subset B), then for the odd cells and the even nodes again. Finally, a half-size time step closed the calculation of the values, as one can see in Figure 1b. In each stage, we used the latest available u values of the neighbors, which meant that the constructed methods were fully explicit and the previous values did not need to be stored at all. Thus, we had a structure consisting of five stages, which corresponded to five partial time steps, which altogether spanned two time steps for odd and even cells, too.

By applying the well-known central difference formula

$\frac{{\partial }^2}{\partial x^2}f(x_i,t_j)\approx \frac{\frac{f(x_{i+1},t_j)-f(x_i,t_j)}{\Delta x}+\frac{f(x_{i-1},t_j)-f(x_i,t_j)}{\Delta x}}{\Delta x}$

to Equation (1) in one dimension, a system of ordinary differential equations (ODEs) were obtained for nodes $i=1,\dots ,N-2$:

(3)$\frac{d{u}_i}{dt}=\alpha \frac{u_{i-1}-2u_i+u_{i+1}}{\Delta x^2}.$

The form of this equation for the first and last node depended on the concrete boundary conditions, which will be discussed later. We defined a matrix M with the following elements:

(4)$m_{ii}=-\frac{2\alpha }{\Delta x^2} (1<i<N), m_{i,i+1}=\frac{\alpha }{\Delta x^2} (1\le i<N),\hspace{1em} m_{i,i-1}=\frac{\alpha }{\Delta x^2} (1<i\le N),$

which is tridiagonal in the currently discussed one-dimensional case. Now, equation-system Equation (3) can be written into a condensed matrix-form:

(5)$\frac{d\overrightarrow{u}}{dt}=M\overrightarrow{u}.$

We recall that the following general time discretization

$\frac{u_i^{n+1}-u_i^n}{\Delta t}=\frac{\alpha }{\Delta x^2}\left[\theta \left(u_{i-1}^n-2u_i^n+u_{i+1}^n\right)+\left(1-\theta \right)\left(u_{i-1}^{n+1}-2u_i^{n+1}+u_{i+1}^{n+1}\right)\right],$

leads to the so-called theta method:

(6)$u_i^{\mathrm{n}+1}=u_i^{\mathrm{n}}+r\left[\theta \left(u_{i-1}^{\mathrm{n}}-2u_i^{\mathrm{n}}+u_{i+1}^{\mathrm{n}}\right)+\left(1-\theta \right)\left(u_{i-1}^{\mathrm{n}+1}-2u_i^{\mathrm{n}+1}+u_{i+1}^{\mathrm{n}+1}\right)\right],$

where $r=\frac{\alpha h}{\Delta x^2}=-\frac{m_{ii}h}{2}>0 , 0<i<N-1$ is the usual mesh ratio and $\theta \in \left[0,1\right]$. For $\theta =0, \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right., \mathrm{and} 1$ one obtains the implicit Euler, the Crank–Nicolson, and the explicit Euler (or, more concretely, the forward-time central-space, FTCS) schemes, respectively [39]. If $\theta <1$, the theta method is implicit. Now, in our shifted-hopscotch scheme, the neighbors were always taken into account at the same, latest time level. Thus, we inserted $u_{i\pm 1}^m$ into Equation (6) instead of $u_{i\pm 1}^n$ and $u_{i\pm 1}^{n+1}$, where $m=n, n+\frac{1}{2} , \mathrm{or} n+1$ at the first, middle, and last stages, respectively (see the colored rows in Figure 1b as well). Now, instead of Equation (6) we can write

(7)$u_i^{n+1}=u_i^n-2r\theta u_i^n-2r\left(1-\theta \right)u_i^{n+1}+r\left(u_{i-1}^m+u_{i+1}^m\right),$

i.e., our final formula reads as follows:

(8)$u_i^{n+1}=\frac{\left(1-2r\theta \right)u_i^n+r\left(u_{i-1}^m+u_{i+1}^m\right)}{1+2r\left(1-\theta \right)}.$

In the case of $\theta =0$, this formula gives back the UPFD method [40,41] with m = n, which takes the form for a half and a full time step, respectively:

(9)$u_i^{n+1}=\frac{u_i^n+\raisebox{1ex}{$r$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(u_{i-1}^m+u_{i+1}^m\right)}{1+r} ,u_i^{n+1}=\frac{u_i^n+r\left(u_{i-1}^m+u_{i+1}^m\right)}{1+2r}.$

The other formula we used was the constant neighbor (CNe) method, which was introduced in our papers [42,43] and now briefly restated here. The starting point is Equation (3), where an approximation is made: When the new value of a variable $u_i^{n+1}$ was calculated, we neglected the fact that the neighbors $u_{i-1}^n$ and $u_{i+1}^n$ were also changing during the time step. It means that the values of u_j $\left(j\ne i\right)$ were considered as constants (that is why we call it constant-neighbor method), so a set of uncoupled ODEs remained:

(10)$\frac{d{u}_i(t)}{dt}=a_i+m_{ii}{u}_i(t),$

where the unknowns u_i are still functions of the (continuous) time, and a quantity $a_i$ is introduced to collect information about the neighbors of cell i:

$a_i={\displaystyle \sum _{j\in \mathrm{neighbours}(i)}{m}_{ij}}{u}_j^n=\alpha \frac{u_{i-1}^n+u_{i+1}^n}{\Delta x^2}=\frac{r}{h}\left(u_{i-1}^n+u_{i+1}^n\right).$

These $a_i$ quantities were considered as constant during one time step, but must be recalculated after each time step. The CNe method takes the analytical solution of these simple Equation (10) at the end of the time step or stage to obtain the new values of the u variable:

(11)$u_i^{n+1}=u_i^n\cdot e^{m_{ii}h}-\frac{a_i}{m_{ii}}\left(1-e^{m_{ii}h}\right).$

At the one-dimensional case, this can be written as follows:

(12)$u_i^{n+1}=u_i^n\cdot e^{- 2r}+\frac{u_{i-1}^n+u_{i+1}^n}{2}\left(1-e^{- 2r}\right).$

At this point, we extended the used methods to general cases where, as happens frequently in real-life problems, the quantities α, k, c, and ρ, describing material properties, are not constants but functions of the space (and time) variables. The interested reader can find more details about this treatment for the case of petroleum reservoirs in Chapter 5 of the book [44]. Let us start with a one-dimensional, still equidistant grid to solve Equation (2). We again discretized the second-order space derivatives by the central difference formula, but now k is a function of x and cannot be simply merged with c and ρ. Thus, we have

${c(x_i)\rho (x_i)\frac{\partial u}{\partial t}|}_{x_i}=\frac{1}{\Delta x} \left[k\left(x_i+\frac{\Delta x}{2}\right)\frac{u(x_i+\Delta x)-u(x_i)}{\Delta x}+k\left(x_i-\frac{\Delta x}{2}\right)\frac{u(x_i-\Delta x)-u(x_i)}{\Delta x}\right] .$

Now, we change from node to cell variables, and $u_i$, $c_i$, and ${\rho }_i$ will be the (average) temperature, specific heat, and density of cell i, respectively, while $k_{i,i+1}$ will be the heat conductivity between cell i and its (right) neighbor. The previous formula will have the form

$\frac{d{u}_i}{dt}=\frac{1}{c_i{\rho }_i\Delta x} \left(k_{i,i+1}\frac{u_{i+1}-u_i}{\Delta x}+k_{i-1,i}\frac{u_{i-1}-u_i}{\Delta x}\right)$

We go on with a still one-dimensional, but non-equidistant, grid with non-uniform cross section. Let us denote by $\Delta x_i$ and $A_i$ the length and the (average) cross section of the cell. One can write the distance between the cell-center of the cell and its arbitrary neighbor j as $d_{ij}=\left(\Delta x_i+\Delta x_j\right)/2$, and let us approximate the area of the interface between the two neighboring cells as $A_{i,i\pm 1}\approx A_i\approx A_{i\pm 1}$. Using these we can write, more generally than above, that

$\frac{d{u}_i}{dt}=\frac{1}{c_i{\rho }_i\Delta x_i{A}_i} \left(A_{i,i+1}{k}_{i,i+1}\frac{u_{i+1}-u_i}{d_{i,i+1}}+A_{i,i-1}{k}_{i,i-1}\frac{u_{i-1}-u_i}{d_{i,i-1}}\right)$

Now, the volume and the heat capacity of the cell is $V_i=A_i \Delta x_i$, and $C_i=c_i{\rho }_i{V}_i$, respectively, while the inverse thermal resistance between these cells can be approximated as $1/R_{ij}\approx k_{ij}{A}_{ij}/d_{ij}$. It is obvious that this quantity is zero between non-neighboring cells. Using these new quantities, we obtain the expression for the time derivative of each cell-variable:

$\frac{d{u}_i}{dt}=\frac{u_{i-1}-u_i}{R_i{{}_{-1,}}_i{C}_i}+\frac{u_{i+1}-u_i}{R_i{{}_{+1,}}_i{C}_i},$

which can be easily generalized further to the case with arbitrary number of neighbors:

(13)$\frac{d{u}_i}{dt}={\displaystyle \sum _{j\ne i}\frac{u_j-u_i}{R_i{}_{,j}{C}_i}}.$

Let us introduce the characteristic time or time constant τ_i of cell i as follows:

(14)${\tau }_i\doteq {\left({\displaystyle \sum _{j\ne i}\frac{1}{R_i{}_{,j}{C}_i}}\right)}^{-1}=\frac{-1}{m_{ii}}.$

Here, m_ii is a diagonal element of the matrix M if we write Equation (13) into the same matrix form as in Equation (5). To restate the used formulas for the general case Equation (13), we introduce the following notations

(15)$r_i=\frac{h}{{\tau }_i} \mathrm{and} A_i=h{\displaystyle \sum _{j\ne i}{m}_{ij}{u}_j^m}=h{\displaystyle \sum _{j\ne i}\frac{u_j^m}{C_i{R}_{ij}}}.$

Now, the general form of Equation (7) is

$u_i^{n+1}=u_i^n-r_i\theta u_i^n-r_i\left(1-\theta \right)u_i^{n+1}+h{\displaystyle \sum _{j\ne i}\frac{u_j^m}{C_i{R}_{ij}}};$

thus, the generalized theta method for integer time steps reads as follows:

(16)$u_i^{n+1}=\frac{\left(1-r_i\theta \right)u_i^n+A_i}{1+r_i\left(1-\theta \right)}.$

Similarly, the generalized CNe formula is

(17)$u_i^{n+1}=u_i^n\cdot e^{- r_i}+\frac{A_i}{r_i}\left(1-e^{- r_i}\right),$

and, of course, for halved time steps, r_i and A_i must be divided by 2.

For the sake of brevity, we will use a compact notation of the individual combinations, where five data are given in a bracket and the numbers are the values of the parameter θ, while the letter ‘C’ is for the CNe constant-neighbor method. For example, (¼, ½, C, ½, ¾) means the following five-stage algorithm will be selected into the top five algorithm in Section 3.3 and named as S2.

Example 1.

Algorithm S2(¼, ½, C, ½, ¾), general from.

Stage 1. Take a half time step with the Equation (16) formula with θ = ¼ for odd cells:

$u_i^{n+1}=\frac{\left(1-\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\right)u_i^n+A_{i,\mathrm{half}}}{1+\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(1-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\right)}, A_{i,\mathrm{half}}=\frac{h}{2}{\displaystyle \sum _{j\ne i}\frac{u_j^m}{C_i{R}_{ij}}}.$

Stage 2. Take a full time step with the Equation (16) formula with θ = ½ for even cells:

$u_i^{n+1}=\frac{\left(1-\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)u_i^n+A_i}{1+r_i\left(1-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}, A_i=h{\displaystyle \sum _{j\ne i}\frac{u_j^m}{C_i{R}_{ij}}}.$

Stage 3. Take a full time step with the Equation (17) formula for odd cells:

$u_i^{n+1}=u_i^n\cdot e^{- r_i}+\frac{A_i}{r_i}\left(1-e^{- r_i}\right), A_i=h{\displaystyle \sum _{j\ne i}\frac{u_j^m}{C_i{R}_{ij}}}.$

Stage 4. The same as Stage 2.

Stage 1. Take a half time step with the Equation (16) formula with θ = ¾ for odd cells:

$u_i^{n+1}=\frac{\left(1-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$8$}\right.r_i\right)u_i^n+A_{i,\mathrm{half}}}{1+\raisebox{1ex}{$r_i$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(1-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.\right)}, A_{i,\mathrm{half}}=\frac{h}{2}{\displaystyle \sum _{j\ne i}\frac{u_j^m}{C_i{R}_{ij}}}.$

All other combinations can be constructed in this manner straightforwardly.

3. Results

3.1. General Definitions and Circumstances

If not stated otherwise, we examined 2-dimensional, rectangle-structured lattices with $N=N_x\times N_z$ cells, similar to what can be seen in Figure 2. We solved Equation (13) subjected to randomly generated initial conditions $u_i(0)=rand$, where rand is a (pseudo)random number with a uniform distribution in the interval (0, 1), generated by the MATLAB for each cell. We also generated different random values for the heat capacities and for the thermal resistances, but with a log-uniform distribution:

(18)$C_i={10}^{({\alpha }_C-{\beta }_C\times rand)} , R_{x,i}={10}^{({\alpha }_{Rx}-{\beta }_{Rx}\times rand)} , R_{z,i}={10}^{({\alpha }_{Rz}-{\beta }_{Rz}\times rand)},$

where the coefficients ${\alpha }_{C }, \dots , {\beta }_{Rz}$ in the exponents were concretized later.

We used zero Neumann boundary conditions, i.e., the system was thermally isolated. This was implemented naturally at the level of Equation (13) since it was enough to omit those terms of the sum, which had infinite resistivity in the denominator due to the isolated border. This implied that the system matrix M had one zero eigenvalue, belonging to the uniform distribution of temperatures, and all other eigenvalues must be negative. Let us denote by ${\lambda }_{\mathrm{MIN}} \left({\lambda }_{\mathrm{MAX}}\right)$ the (nonzero) smallest (largest) absolute value eigenvalues of matrix M. The stiffness ratio of the system was defined as ${\lambda }_{\mathrm{MAX}}/{\lambda }_{\mathrm{MIN}}$. The maximum possible time step size for the FTCS (explicit Euler) scheme (from the point of view of stability) was exactly calculated as $h_{\mathrm{MAX}}^{\mathrm{FTCS}}=\left|2/{\lambda }_{\mathrm{MAX}}\right|$, above which the solutions were expected to blow up. This threshold time step $h_{\mathrm{MAX}}^{\mathrm{FTCS}}$ (often called CFL limit) was essentially valid for the second-order explicit Runge–Kutta (RK) methods as well, since, in the case of them, the border of the stability domain in the negative real axis was the same as for the first-order RK method [45]. We used these two numbers to characterize the “difficulty level” of the problem.

We calculated the numerical error by comparing our numerical solutions $u_{\mathrm{j}}^{\mathrm{num}}$ with the reference solution $u_{\mathrm{j}}^{\mathrm{ref}}$ at final time $t_{\mathrm{fin}}$. In Section 3.5, the reference solution will be an analytical solution; otherwise, it is a very accurate numerical solution, which has been calculated by the ode15s built-in solver of MATLAB with very strict error tolerance. We used the following three types of (global) error. The first one was the maximum of the absolute differences:

(19)$\mathrm{Error}(L_{\infty })=\underset{0\le \mathrm{j}\le N}{\max }\left|u_{\mathrm{j}}^{\mathrm{ref}}(t_{\mathrm{fin}})-u_{\mathrm{j}}^{\mathrm{num}}(t_{\mathrm{fin}})\right|.$

The second one was the average absolute error:

(20)$\mathrm{Error}(L_1)=\frac{1}{N}{\displaystyle \sum _{0\le \mathrm{j}\le N}\left|u_{\mathrm{j}}^{\mathrm{ref}}(t_{\mathrm{fin}})-u_{\mathrm{j}}^{\mathrm{num}}(t_{\mathrm{fin}})\right|}.$

The third one gave the error in terms of energy in case of the heat equation. It took into account that an error of the solution in a cell with a large volume or heat capacity had more significance in practice than in a very small cell

(21)$\mathrm{Error}(Energy)=\frac{1}{N}{\displaystyle \sum _{0\le \mathrm{j}\le N}{C}_{\mathrm{j}}\left|u_{\mathrm{j}}^{\mathrm{ref}}(t_{\mathrm{fin}})-u_{\mathrm{j}}^{\mathrm{num}}(t_{\mathrm{fin}})\right|}.$

It is well known that the true solution always follows the maximum and minimum principles [39] (p. 87). We say a method is positivity preserving if it never violates this principle, i.e., in our case no value of u was outside of the $\left[0,1\right]$ interval. We were interested in how these errors depend on the time step size in different concrete situations. As one can see in Figure 1, there were five time steps (five stages) altogether, instead of four in the shifted hopscotch structure. Therefore, for the sake of honesty, we had to calculate the effective time step size as $h_{\mathrm{EFF}}=\raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$5$}\right.h$, and the errors were plotted as a function of this quantity.

We used a desktop computer with an Intel Core i5-9400 as CPU, 16.0 GB RAM, for the simulations with the MATLAB R2020b software, in which there was a built-in tic-toc function to measure the running time.

3.2. Preliminary Tests

We applied the following nine different values for parameter theta: $\theta \in \left\{0, \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$5$}\right., \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$4$}\right., \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$3$}\right., \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right., \raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right., \raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right., \raisebox{1ex}{$4$}\!\left/ \!\raisebox{-1ex}{$5$}\right., 1 \right\}$ in Equation (16). It means that, together with the CNe formula, we had 10 different formulas and we inserted all of these into the shifted-hopscotch structure in all possible combinations. As there were five stages in the structure, we had 10⁵ = 100,000 different algorithm-combinations. We wrote a code, which systematically constructed and tested all these combinations. After some tests, we chose the few best combinations, and continued the work only with them. For this, we needed an automatic assessment of the performance of the combinations. The difficulty was in the fact that methods that are very inaccurate or even unstable for large time step sizes can be the most accurate for small time step sizes. Therefore, we chose two different final times, $t_{\mathrm{fin}}=0.1 , 10$, and calculated the solution with a large time step size (typically $t_{\mathrm{fin}} /4$), then repeated the calculation for subsequently halved time step sizes R times until h reached a small value (typically around $2\times {10}^{-6}$). We introduced aggregated relative error (ARE) quantities for each type of errors defined above, which were calculated for the $L_{\infty }$ error as follows:

(22)$\mathrm{ARE}(L_{\infty })=\frac{1}{R}{\displaystyle \sum _{i=1}^R\left(\log {\left(\mathrm{Error}(L_{\infty })\right)}_{\mathrm{OEH}}-\log {\left(\mathrm{Error}(L_{\infty })\right)}_{\mathrm{shifted}}\right)},$

which means that $\mathrm{ARE}(L_{\infty })$ was the average of the difference between the error of the original OEH method and the actual shifted combination in terms of orders of magnitude. Then the code calculated the simple average of these errors:

(23)$\mathrm{ARE}=\frac{1}{3}\left(\mathrm{ARE}(L_{\infty })+\mathrm{ARE}(L_1)+\mathrm{ARE}(Energy)\right),$

and finally sorted the 100,000 combinations in descending order, according to this quantity. In the obtained list, usually positive ARE values were assigned the first few thousands of combinations, with the largest ones typically around 2, which means that some combinations were roughly two orders of magnitude more accurate than the original OEH method. We performed this procedure in the case of four different small systems with $N_x\times N_z=2\times 2 , 2\times 6 , 4\times 4 , \mathrm{and} 3\times 5$. The parameters ${\alpha }_{C }, {\beta }_C , {\alpha }_{Rx} , {\beta }_{Rx} , {\alpha }_{Rz} , {\beta }_{Rz}$ of the distribution of the mesh-cells data were chosen to construct test problems with various stiffness ratios and $h_{\mathrm{MAX}}^{\mathrm{FTCS}}$, for example, ${\alpha }_{C }=1, 2, \mathrm{or} 3 , {\beta }_C=2, 4, \mathrm{or} 6$. We give the best 12 combinations in their short form:

(24) (0, ½, ½, ½, 1), (½, ½, ½, ½, ½), (0, C, ½, C, 1), (0, C, C, C, 1),(¾, ⅔, ½, ⅓, ¼), (¼, ½, C, ½, ¾), (⅓, ⅔, C, ⅓, ⅔), (C, ½, C, ½, C), (⅕, ½, ½, ½, ⅘), (¼, ½, ½, ½, ¾), (⅓, ½, ½, ½, ⅔), (0, ½, ½, C, 1).

Later, we proved that formulas θ = 1 and CNe preserved positivity of the solution and, therefore, if only these two formulas are used in a combination, the whole algorithm will preserve positivity. Since this property is considered valuable [40,46], we repeated the above experiments for these 2⁵ = 32 combinations (instead of the 100,000 above). We concluded that the (C, C, C, C, C) combination was the most accurate among these; therefore, we further investigate d13 combinations altogether. We emphasize that these are the results of only preliminary (one might say tentative) tests, with the sole purpose of reducing the huge number of combinations into a manageable number, and we have not stated anything exactly until this point.

3.3. Case Study I and Comparison with Other Solvers

We examined a grid similar to the one in Figure 2 with isolated boundary, but the sizes were fixed to $N_x=100$ and $N_z=100$; thus, the total cell number was 10,000, while the final time was $t_{\mathrm{fin}}=0.1$. The exponents introduced above were set to the following values

(25)${\alpha }_C=2, {\beta }_{\mathrm{C}}=4, {\alpha }_{Rx}={\alpha }_{Rz}=1, {\beta }_{Rx}={\beta }_{Rz}=2,$

which means that log-uniformly distributed values between 0.01 and 100 were given to the capacities. The generated system can be characterized by its stiffness ratio and $h_{\mathrm{MAX}}^{\mathrm{FTCS}}$ values, which are $3.1\times {10}^7$ and $7.3\times {10}^{-4}$, respectively. The performance of new algorithms was compared with all available MATLAB solvers:

• ode15s, a first- to fifth-order (implicit) numerical differentiation formula with variable step and variable order (VSVO), developed for solving stiff problems;

• ode23s, a second-order modified (implicit) Rosenbrock formula;

• ode23t, applies (implicit) trapezoidal rule with using free interpolant;

• ode23tb, combines backward differentiation formula and trapezoidal rule;

• ode45, a fourth/fifth-order explicit Runge–Kutta–Dormand–Prince formula;

• ode23, second/third-order explicit Runge–Kutta–Bogacki–Shampine method;

• ode113, first- to 13th- order VSVO Adams–Bashforth–Moulton numerical solver.

For all used MATLAB solvers, tolerances were changed over many orders of magnitude, from the maximum value $‘\mathrm{AbsTol}’=‘\mathrm{RelTol}’ \doteq ‘\mathrm{Tol}’={10}^3$ to the minimum value $‘\mathrm{AbsTol}’=‘\mathrm{RelTol}’ \doteq ‘\mathrm{Tol}’={10}^{-5}$. We also used the following methods for comparison purposes. The original UPFD, the CNe method, and the Dufort–Frankel (DF) algorithm

(26)$u_i^{n+1}=\frac{\left(1-r_i\right)u_i^{n-1}+2A_i}{1+r_i}$

are explicit unconditionally stable schemes. The DF scheme is a two-step method, which is not self-starter; thus, we calculated the first time step by two subsequent UPFD steps with halved time step sizes. Finally, the widely used Crank–Nicolson (CN) is a widely used implicit scheme:

$\underset{A}{\underbrace{\left(I-\frac{h}{2}M\right)}}{\overrightarrow{u}}^{n+1}=\underset{B}{\underbrace{\left(I+\frac{h}{2}M\right)}}{\overrightarrow{u}}^n ,$

where I is the unit matrix. The A and B matrixes did not depend on time here, so they were calculated only once, before the first time step. This scheme was implemented two different ways: by calculating $Y=A^{-1}B$ before the first time step and then each time step was just a simple matrix multiplication ${\overrightarrow{u}}^{n+1}=Y {\overrightarrow{u}}^n$, which is denoted by ‘CrN invert’, and by the ${\overrightarrow{u}}^{n+1}=\mathrm{linsolve}\left(A,B{\overrightarrow{u}}^n \right)$ command, which is denoted by ‘CrN lins’. One will see that this implicit method was still very slow for the actual system size.

We plotted the $L_{\infty }$, the $L_1$, and the energy errors as a function of the effective time step size $h_{\mathrm{EFF}}$, and, based on this (and on similar data belonging to the more stiff system in the next subsection), we selected the following top five combinations from those listed in Equation (24) and after that:

S1 (C, C, C, C, C),S2 (¼, ½, C, ½, ¾),S3 (¼, ½, ½, ½, ¾),S4 (0, ½, ½, ½, 1),S5 (0, ½, ½, C, 1)

In Figure 3, we present the error functions only for these top five combinations and some other, known methods, while in Figure 4 and Figure 5 one can see the $L_{\infty }$ and the energy errors vs. the total running times. Furthermore, Table 1 presents some results that were obtained by our numerical schemes, the two implementations of the Crank–Nicolson method and the “ode” routines of MATLAB.

3.4. Case Study II and Comparison with Other Solvers

We tested our new algorithms and the conventional solvers for a harder problem as well. Thus, new values were set for the α and β exponents:

(27)${\alpha }_C=3, {\beta }_{\mathrm{C}}=6, {\alpha }_{Rx}=3, {\alpha }_{Rz}=1, {\beta }_{Rx}={\beta }_{Rz}=4.$

This means that the width of the distribution of the capacities and thermal resistances were increased. Now, the largest capacity was six orders or magnitude larger than the smallest one. The system acquired some anisotropy as well, since the resistances in the x direction were two orders of magnitude larger than in the z direction, on average, $R_{x,i}\in \left[{10}^{-1}, {10}^3\right] , R_{z,i}\in \left[{10}^{-3}, {10}^1\right]$. With this modification we gained a system with a much higher stiffness ratio, $2.5\times {10}^{11}$, while the maximum allowed time step size for the standard FTCS was $h_{\mathrm{MAX}}^{\mathrm{EE}}=1.6\times {10}^{-6}$. All other parameters and circumstances remained the same as in Section 3.3. In Figure 6 and Figure 7, energy errors and the $L_{\infty }$ errors are presented as a function of the time step size and the total running time, respectively, while in Table 2 some results are collected for the same system.

One can see that the implicit MATLAB solvers performed much better than the explicit ones, unlike in the previous case. As all of our cases had considerable stiffness, the implicit Crank–Nicolson method was more accurate than our explicit methods, but the difference was not always large. On the other hand, all implicit methods required much longer time for running than the unconditionally stable explicit algorithms. However, we remind the reader that this system of 10,000 cells was still rather small, and, for a growing number of cells, the implicit methods have an increasing disadvantage. Moreover, the performance of the standard explicit solvers declined with increasing stiffness, so solving very large and stiff systems is a very tough problem. In Figure 7, the otherwise very powerful ode45 solver is represented only by one point (right side of the figure) because, for large tolerances, it diverged, while for low tolerances, the running of the program was stopped by an error message. In fact, we did not test our methods in the case of larger systems (with 20,000 or more cells) because no solvers of MATLAB could produce a reference solution.

3.5. Verification by Comparison to Analytical Results

We considered very recent nontrivial analytical solutions of Equation (1) found by Barna and Mátyás [10] by a similarity transformation technique. Both of them are given, on the whole, a real number line for positive values of t, as follows

(28)$u_1^{exact}(x,t)=\frac{x}{t^{3/2}}{e}^{-\frac{x^2}{4\alpha t}},$

and

(29)$u_2^{exact}=\frac{x}{t^{5/2}}\left(1-\frac{x^2}{6\alpha t}\right)e^{-\frac{x^2}{4\alpha t}}.$

We reproduced these solutions only in finite space and time intervals $x\in \left[x_1,x_2\right]$ and $t\in \left[t_0, t_{\mathrm{fin}}\right]$, where $x_1=-5 , x_2=5 , t_0=0.5 , t_{\mathrm{fin}}=1$. The space interval was discretized by creating nodes as follows: $x_j=x_1+j\Delta x , j=0,\dots ,1000 , \Delta x=0.01$. We prescribed the appropriate Dirichlet boundary conditions at the two ends of the interval:

(30)$u_1(x=x_b,t)=\frac{x_b}{t^{3/2}}{e}^{-\frac{x_b^2}{4\alpha t}},$

and

(31)$u_2(x=x_b,t)=\frac{x_b}{t^{5/2}}\left(1-\frac{x_b^2}{6\alpha t}\right)e^{-\frac{x_b^2}{4\alpha t}},$

where $x_b\in \left\{x_1 , x_2\right\}$. We obtained that the new methods were convergent and the order of convergence was 2. In Figure 8, the $L_{\infty }$ errors as a function of the effective time step size h_EFF are presented for the case of the u₂ solution for the top five algorithms and a first-order “reference curve” for the original CNe method. We note that very similar curves were obtained for the u₁ solution, as well as for other space and time intervals.

We also present the graphs of the initial function $u_2(x, t_0=0.5)$, the analytical solution $u_2^{exact}\left(x,t_{\mathrm{fin}}\right)$, and two corresponding numerical solutions for $h=2\times {10}^{-4}$ in Figure 9. We emphasize that, at this time step size range, the algorithm S4 was the most accurate, but we did not present its graph, as it would be completely indistinguishable from the exact solution.

3.6. Verification for the Nonlinear Fisher’s Equation

We considered

$\frac{\partial u}{\partial t}=\alpha {\nabla }^{2}u+\beta u\left(1-u\right),$

for α = 1 and solved it, subject to the following initial condition:

$u(x,t=0)={\left(1+e^{\sqrt{\frac{\beta }{6}} x}\right)}^{-2}.$

The analytical solution of this equation was known [11,12] and can be given as

$u^{exact}(x,t)={\left(1+e^{\sqrt{\frac{\beta }{6}} x- \frac{5}{6}\beta t}\right)}^{-2},$

The Dirichlet boundary conditions appropriate to the exact solution were prescribed at both ends of the interval:

$u(x=x_0,t)={\left(1+e^{\sqrt{\frac{\beta }{6}} x_0- \frac{5}{6}\beta t}\right)}^{-2} , \mathrm{and}\hspace{1em}u(x=x_{\mathrm{fin}},t)={\left(1+e^{\sqrt{\frac{\beta }{6}} x_{\mathrm{fin}}- \frac{5}{6}\beta t}\right)}^{-2}.$

We reproduced this analytical solution in the case of several values of the parameters $\beta , x_0, x_{\mathrm{fin}}, t_{\mathrm{fin}}$, and $\Delta x$, but here we present results only for $\beta =2, x_0=0, x_{\mathrm{fin}}= 5 , t_{\mathrm{fin}} = 1$, and $\Delta x=0.01$. In order to have good stability properties, we took account the effect of the nonlinear term via the following semi-implicit treatment

$u_{\mathrm{i}}^{\mathrm{new}}=u_{\mathrm{i}}^{\mathrm{pred}}+\beta u_{\mathrm{i}}^{\mathrm{pred}}\left(1-u_{\mathrm{i}}^{\mathrm{new}}\right)h,$

which can be rearranged into a fully explicit form:

$u_{\mathrm{i}}^{\mathrm{new}}=\frac{1+\beta h}{1+\beta h{u}_{\mathrm{i}}^{\mathrm{pred}}}{u}_{\mathrm{i}}^{\mathrm{pred}},$

where $u_{\mathrm{i}}^{\mathrm{pred}}$ is the actual value of $u_i$ after performing the previous stage calculations, which deal with the diffusion term. This operation was performed in two extra, separate stages, where a loop was going through all the nodes (both odd and even). The first extra stage was after (the original) Stage 2, while the second was after Stage 4, at the same moments in which the boundary conditions were refreshed.

With this procedure we obtained that the new methods behaved very similarly as in the linear case. In Figure 10, the $L_{\infty }$ errors as a function of the effective time step size h_EFF are presented for the top five algorithms and the original CNe method. We note that very similar curves were obtained for other space and time intervals and values of parameter β.

We emphasize that this subsection is only to demonstrate that these shifted-hopscotch methods can be successfully applied for other, possibly nonlinear, equations as well, and there is some chance that the relative strength of the S1–S5 algorithms is preserved. The systematic investigation of the behavior of the new methods in case of nonlinear equations was out of the scope of this paper and plans to be published in future works.

3.7. Analytical Investigations

We started with the analysis of the convergence properties of the five most advantageous methods in the one-dimensional case for constant values of mesh ratio r: In other words, when the equidistant spatial discretization was fixed.

Theorem 1.

The order of convergence of the S1…S5 shifted hopscotch algorithms is 2 for the Equation (5) system of linear ODEs:

(32)$\frac{d\overrightarrow{u}}{dt}=M\overrightarrow{u} ,\hspace{1em}\overrightarrow{u}(\mathrm{t}=0)={\overrightarrow{u}}^0,$

where M is defined in Equation (4) and ${\overrightarrow{u}}^0$ is an arbitrary vector.

The proof is presented in the Appendix A. Now, we turn our attention to the question of stability. First, we examined the two types of formulas (the UPFD and the CNe), about which we already stated that they are positivity preserving.

Lemma 1.

The new$u_i^{n+1}$values are the convex combinations of the old$u_i^n$ values in Equations (9), (11) and (17) formulas, as well as in the Equation (16) formula for$\theta =0$.

The proofs of these statements are very easy, and they are given in our previous papers (see Theorem 2 in [37] for the UPFD and also Theorem 2 in [43] for the CNe). We also recalled a lemma on the associativity of convex combinations [47] (p. 28).

Lemma 2.

A convex combination$x={\displaystyle \sum a_i{x}_i}$of convex combinations$x_i={\displaystyle \sum b_{ij}{y}_{ij}}$ is again a convex combination:

$x={\displaystyle \sum {\displaystyle \sum \left(a_i{b}_{ij}\right)y_{ij}}}$

for any $y_{ij}\in {ℝ}^n$.

Corollary 1.

All shifted hopscotch algorithms containing only the CNe and the UPFD formulas mentioned in Lemma 1, the new $u_i^{n+2}$values are the convex combinations of the old$u_i^n$ values.

This immediately implies that all algorithms mentioned in Corollary 1, especially S1 (C, C, C, C, C), followed the minimum and maximum principle, and, therefore, preserved positivity, which is a much stronger property than unconditional stability.

We examined the stability of the S2–S5 algorithms through their eigenvalues. We constructed the five stages of these algorithms in a matrix form. If we denote these matrices by M₁,…, M₅, respectively, the matrix of a two-step, five-stage algorithm can be written as $H=M_5\left(M_4\left(M_3\left(M_2{M}_1\right)\right)\right)$, and ${\overrightarrow{u}}^{n+2}=H {\overrightarrow{u}}^n$. We calculated and plotted the absolute values of the eigenvalues of this matrix H as a function of the mesh ratio r for several values of N, and found that the largest one was always exactly 1, which implies that no perturbation vector $\overrightarrow{\epsilon }$ could grow unboundedly. In Figure 11, we present here only the eigenvalue functions of S4 (0, ½, ½, ½, 1), which produced the largest errors for large time step sizes (see, e.g., the top right corner of Figure 8) and, therefore, one might question its stability first. These results verify that the S2…S4 methods were unconditionally stable, even if they could not be considered as exact analytical proofs.

4. Discussion and Summary

In this article, we constructed and tested novel numerical algorithms to solve the non-stationary diffusion (or heat) equation. The new algorithms were fully explicit time-integrators obtained by applying half and full time steps in the odd-even hopscotch structure. All of the algorithms consisted of five stages, but they were one-step methods in the sense that when the new values of the unknown function u were calculated, only the most recently calculated u values were used. Thus, the methods could be implemented such that only one array of storage was required for the u variable, which means that the memory requirement was very low. We applied the conventional theta method with nine different values of $\theta$ and the non-conventional CNe method to construct 10⁵ combinations and chose the top five of them via numerical experiments. We showed concrete results of these five algorithms in the case of two 2-dimensional stiff systems containing 10,000 cells with highly inhomogeneous, randomly generated parameters and discontinuous initial conditions. These experiments suggest that the proposed methods are, indeed, competitive, as they can give fairly accurate results orders of magnitude faster than the well-optimized MATLAB routines or the Crank–Nicolson method, and they are also significantly more accurate for stiff systems than the UPFD, the Dufort–Frankel, or the original odd-even hopscotch method. We think that, if high accuracy is required, the S4 (0, ½, ½, ½, 1) combination can be proposed; however, when preserving positivity is crucial, the S1 (C, C, C, C, C) algorithm should be used.

We verified these new algorithms by reproducing a recently found nontrivial analytical solution of the heat equation, and then also showed that the nonlinear Fisher’s equation can be solved by them. We proved analytically that the new methods are second order for fixed spatial discretization. It must be noted that the original explicit Euler, UPFD, and CNe formulas are only first order and the shifted-hopscotch structure made their combinations second order. On the other hand, the original theta method for $\theta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ (which is the trapezoidal rule) is second order, but it is either an implicit method (the CN scheme) or, if made explicit by a predictor step, it requires an extra (predictor) stage and, at the same time, becomes only conditionally stable. At the end of the paper, we conveyed evidence about the good stability properties of the new methods.

In our next papers we would like to extend our investigations in order to find even more effective algorithms and perform more comprehensive tests to clarify which are the best under specific circumstances. Then, we also plan to apply the methods for more complicated, mostly nonlinear equations, as well as in case of real-life engineering problems such as heat transfer by conduction and radiation in buildings or machine tools. We already started the application of the methods to equations where additional terms are present besides the diffusion term, such as the advection-diffusion and the Kardar–Parisi–Zhang equations, but these investigations are still at an embryonic stage. In principle, the shifted hopscotch structure and the concrete proposed methods can be adapted to hyperbolic equations like the Euler equations as well, but we have doubts about their stability properties under those circumstances.

Author Contributions

Conceptualization, E.K. and H.K.; methodology, E.K. and M.S.; software, Á.N. and I.O.; validation, E.K., H.K. and I.O.; formal analysis, M.S.; investigation, Á.N., I.O. and H.K.; resources, E.K.; data curation, Á.N.; writing—original draft preparation, E.K. and M.S.; writing—review and editing, E.K. and Á.N.; visualization, Á.N. and I.O.; supervision, E.K.; project administration, E.K. and Á.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data are available at the following link: http://dx.doi.org/10.17632/sdxrxf68jg.1#folder-649cdd07-097d-43de-b61b-87a74f9b2eec (accessed on 24 August 2021).

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Proof of Theorem 1.

First, we give the general form of the new vector elements $u_i^{n+2}$ expressed by the old elements $u_j^n$ for odd and even cells far from the boundaries for the top 5 algorithms.

S1 (C, C, C, C, C), i is odd:

(A1)

$u_i^{n+2}=\left[\begin{array}{l}\frac{1}{8}\left(3e^{-8r}-3e^{-7r}+3e^{-6r}+e^{-5r}+e^{-4r}-e^{-3r}+e^{-2r}+3e^{-r}\right)u_i^n\\ +\left(-\frac{5}{16}{e}^{-8r}+\frac{1}{4}{e}^{-7r}-\frac{1}{8}{e}^{-6r}-\frac{1}{4}{e}^{-5r}+\frac{1}{4}{e}^{-3r}+\frac{1}{8}{e}^{-2r}-\frac{1}{4}{e}^{-r}+\frac{5}{16}\right)\left(u_{i-1}^n+u_{i+1}^n\right)\\ +\left(\frac{1}{4}{e}^{-8r}-\frac{1}{4}{e}^{-7r}+\frac{1}{4}{e}^{-5r}-\frac{1}{4}{e}^{-4r}-\frac{1}{4}{e}^{-3r}+\frac{1}{4}{e}^{-r}\right)\left(u_{i-2}^n+u_{i+2}^n\right)\\ +\left(-\frac{5}{32}{e}^{-8r}+\frac{3}{16}{e}^{-7r}+\frac{1}{16}{e}^{-6r}-\frac{1}{16}{e}^{-5r}+\frac{1}{16}{e}^{-3r}-\frac{1}{16}{e}^{-2r}-\frac{3}{16}{e}^{-r}+\frac{5}{32}\right)\left(u_{i-3}^n+u_{i+3}^n\right)\\ +\frac{1}{16}{\left(e^{-2r}-1\right)}^3\left(e^{-r}-1\right)e^{-r}\left(u_{i-4}^n+u_{i+4}^n\right)-\frac{1}{32}{\left(e^{-2r}-1\right)}^3{\left(e^{-r}-1\right)}^2\left(u_{i-5}^n+u_{i+5}^n\right)\end{array}\right]$

S1 (C, C, C, C, C), i is even:

(A2)

$u_i^{n+2}=\frac{1}{8}\left[\begin{array}{l}\left(3e^{-7r}+e^{-6r}-e^{-5r}+e^{-4r}+e^{-3r}+3e^{-2r}-3e^{-r}+3\right)u_i^n-\left(3e^{-7r}-e^{-5r}+e^{-3r}-3e^{-r}\right)\left(u_{i-1}^n+u_{i+1}^n\right)\\ +2\left(e^{-7r}-e^{-5r}-e^{-4r}+e^{-3r}-e^{-r}+1\right)\left(u_{i-2}^n+u_{i+2}^n\right)-e^{-r}{\left(e^{-2r}-1\right)}^3\left(u_{i-3}^n+u_{i+3}^n\right)\\ +\frac{1}{2}{\left(e^{-2r}-1\right)}^3\left(e^{-r}-1\right)\left(u_{i-4}^n+u_{i+4}^n\right)\end{array}\right]$

S2 (¼, ½, C, ½, ¾), i is odd:

(A3)

$u_i^{n+2}=\frac{1}{(4+r){(1+r)}^2(4+3r)}\left[\begin{array}{l}(r-4)\left(2e^{-2r}{r}^3-e^{-2r}{r}^2+r^3-e^{-2r}r-5r^2-4e^{-2r}-4r\right)u_i^n\\ +\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left((-r^4 +16r^3-3r^2 +16r-16)e^{-2r} +r^4 +16r^3+3r^2 +16r+16\right)\left(u_{i-1}^n+u_{i+1}^n\right)\\ +\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left((-4+r)r((r^2-3r+4)e^{-2r}-r^2-5r-4)\right)\left(u_{i-2}^n+u_{i+2}^n\right)\\ +\left(4r^2(e^{-2r}r-2e^{-2r}+r+2)\right)\left(u_{i-3}^n+u_{i+3}^n\right)\\ +r^3(e^{-2r}-1)(-4+r)\left(u_{i-4}^n+u_{i+4}^n\right)-2r^4(e^{-2r}-1)\left(u_{i-5}^n+u_{i+5}^n\right)\end{array}\right]$

S2 (¼, ½, C, ½, ¾), i is even:

(A4)

$u_i^{n+2}=\frac{1}{{(1+r)}^2(4+3r)}\left[\begin{array}{l}\left((r^3 +5r^2 -4r)e^{-2r} +2r^3+r^2 -r+4\right)u_i^n+\frac{r}{2}\left(r-4\right)(e^{-2r}r-2e^{-2r} -r-2)\left(u_{i-1}^n+u_{i+1}^n\right)\\ -\frac{1}{2}r(e^{-2r}{r}^2-5e^{-2r}r-r^2+4e^{-2r}-3r-4)\left(u_{i-2}^n+u_{i+2}^n\right)\\ +\frac{1}{2}{r}^2(e^{-2r}-1)(r-4)\left(u_{i-3}^n+u_{i+3}^n\right)-r^3(e^{-2r}-1)\left(u_{i-4}^n+u_{i+4}^n\right)\end{array}\right]$

S3 (¼, ½, ½, ½, ¾), i is odd:

(A5)

$u_i^{n+2}=\frac{1}{(4+r){(1+r)}^3(4+3r)}\left[\begin{array}{l}(4-r)(r^4+r^3+9r^2+r+4)u_i^n+r(r^4+19r^2+32)\left(u_{i-1}^n+u_{i+1}^n\right)\\ r^2(4-r)(r^2+r+8)\left(u_{i-2}^n+u_{i+2}^n\right)+24r^3\left(u_{i-3}^n+u_{i+3}^n\right)-2r^4(r-4)\left(u_{i-4}^n+u_{i+4}^n\right)\\ +4r^5\left(u_{i-5}^n+u_{i+5}^n\right)\end{array}\right]$

S3 (¼, ½, ½, ½, ¾), i is even:

(A6)

$u_i^{n+2}=\frac{1}{{(1+r)}^3(4+3r)}\left[\begin{array}{l}\left(r^4-r^3+9r^2-r+4\right)u_i^n+( 4-r)(r^3+2r)\left(u_{i-1}^n+u_{i+1}^n\right)\\ +r^2(r^2-r+8)\left(u_{i-2}^n+u_{i+2}^n\right)+r^3(4-r)\left(u_{i-3}^n+u_{i+3}^n\right)+2r^4\left(u_{i-4}^n+u_{i+4}^n\right)\end{array}\right]$

S4 (0, ½, ½, ½, 1), i is odd:

(A7)

$u_i^{n+2}=\frac{1}{{(1+r)}^4}\left[\begin{array}{l}\left(2r^2+1\right)u_i^n+\frac{r}{2}\left(r^2+4\right)\left(u_{i-1}^n+u_{i+1}^n\right)+2r^2\left(u_{i-2}^n+u_{i+2}^n\right)-\frac{r^3}{4}\left(r^2-6\right)\left(u_{i-3}^n+u_{i+3}^n\right)\\ +\frac{r^4}{2}\left(u_{i-4}^n+u_{i+4}^n\right)+\frac{r^5}{4}\left(u_{i-5}^n+u_{i+5}^n\right)\end{array}\right]$

S4 (0, ½, ½, ½, 1), i is even:

(A8)

$u_i^{n+2}=\frac{1}{{(1+r)}^4}\left[\left(2r^2+1\right)u_i^n+r(r^2+2)\left(u_{i-1}^n+u_{i+1}^n\right)+2r^2\left(u_{i-2}^n+u_{i+2}^n\right)+r^3\left(u_{i-3}^n+u_{i+3}^n\right)+r^4\left(u_{i-4}^n+u_{i+4}^n\right)\right]$

S5 (0, ½, ½, C, 1), i is odd:

(A9)

$u_i^{n+2}=\frac{1}{{(1+r)}^3}\left[\begin{array}{l}\frac{1}{2}\left((2r^2-r)e^{-2r}+2r^2-r+2\right)u_i^n+\frac{r}{8}\left((r^3+2r^2-5r+4)e^{-2r}+r^3-3r+12\right)\left(u_{i-1}^n+u_{i+1}^n\right)\\ -\frac{r}{4}\left({(r-1)}^2{e}^{-2r}+r^2-4r-1\right)\left(u_{i-2}^n+u_{i+2}^n\right)+\frac{r^2}{8}\left(2e^{-2r}r-2r^2-3e^{-2r}+4r+3\right)\left(u_{i-3}^n+u_{i+3}^n\right)\\ -\frac{r^3}{4}(e^{-2r}-1)\left(u_{i-4}^n+u_{i+4}^n\right)-\frac{r^4}{8}(e^{-2r}-1)\left(u_{i-5}^n+u_{i+5}^n\right)\end{array}\right]$

S5 (0, ½, ½, C, 1), i is even:

(A10)

$u_i^{n+2}=\frac{1}{{(1+r)}^3}\left[\begin{array}{l}\frac{1}{2}\left((2-r)e^{-2r}+3r\right)u_i^n+\frac{1}{2}\left((2r-1)e^{-2r}+2r^2+1\right)\left(u_{i-1}^n+u_{i+1}^n\right)-\frac{r^2}{2}(e^{-2r}-1)\left(u_{i-3}^n+u_{i+3}^n\right)\\ +\frac{r}{4}\left((r^2+2r-3)e^{-2r}+r^2+3\right)\left(u_{i-2}^n+u_{i+2}^n\right)-\frac{r^3}{4}(e^{-2r}-1)\left(u_{i-4}^n+u_{i+4}^n\right)\end{array}\right]$

We have to show that the zeroth-, first-, and second-order local errors are zero, i.e., the numerical solution is identical to the exact solution of the PDE initial value problem Equation (A1) up to second order. We will use the usual series expansion up to second order

(A11)

$e^{-x}=1-x+\frac{x^2}{2!}+O\left(x^3\right)$

The exact solution at the end of a doubled time step is the following:

${\overrightarrow{u}}^{n+2}=e^{2Mh}{\overrightarrow{u}}^n=\left(1+2Mh+4M^2\frac{h^2}{2}+O\left(h^3\right)\right){\overrightarrow{u}}^n.$

After some simple algebraic calculations, we obtain:

(A12)

$u_i^{n+2}=\left(2r^2\right)u_{i-2}^n+\left(2r-8r^2\right)u_{i-1}^n+\left(1-4r+12r^2\right)u_i^n+\left(2r-8r^2\right)u_{i+1}^n+\left(2r^2\right)u_{i+2}^n+O\left(r^3\right) .$

Now we calculate the coefficients of Equation (A1) up to the second order:

(A13)

$u_i^{n+2}=\left(1-4r+12r^2\right)u_i^n+\left(2r-8r^2\right)\left(u_{i-1}^n+u_{i+1}^n\right)+\left(2r^2\right)\left(u_{i-2}^n+u_{i+2}^n\right)+{\displaystyle \sum _{i=1}^4{O}_i\left(r^3\right)}.$

Substituting these terms into Equation (A1) we obtain

(A14)

$\begin{array}{l}\frac{1}{8}\left(3e^{-8r}-3e^{-7r}+3e^{-6r}+e^{-5r}+e^{-4r}-e^{-3r}+e^{-2r}+3e^{-r}\right)=1-4r+12r^2+O_1\left(r^3\right)\\ \left(-\frac{5}{16}{e}^{-8r}+\frac{1}{4}{e}^{-7r}-\frac{1}{8}{e}^{-6r}-\frac{1}{4}{e}^{-5r}+\frac{1}{4}{e}^{-3r}+\frac{1}{8}{e}^{-2r}-\frac{1}{4}{e}^{-r}+\frac{5}{16}\right)=2r-8r^2+O_2\left(r^3\right)\\ \left(\frac{1}{4}{e}^{-8r}-\frac{1}{4}{e}^{-7r}+\frac{1}{4}{e}^{-5r}-\frac{1}{4}{e}^{-4r}-\frac{1}{4}{e}^{-3r}+\frac{1}{4}{e}^{-r}\right)=2r^2+O_3\left(r^3\right)\\ \left(-\frac{5}{32}{e}^{-8r}+\frac{3}{16}{e}^{-7r}+\frac{1}{16}{e}^{-6r}-\frac{1}{16}{e}^{-5r}+\frac{1}{16}{e}^{-3r}-\frac{1}{16}{e}^{-2r}-\frac{3}{16}{e}^{-r}+\frac{5}{32}\right)=O_4\left(r^3\right)\end{array}$

It can be immediately seen that this expression is identical to the analytical solution Equation (A12) up to second order. Equation (A2) can be handled similarly to obtain the same expression as in Equation (A12) or (A14), which proves the theorem for the S1 algorithm. For the S2…S5 methods we use the power series expansion ${\left(1+x\right)}^{-1}=1-x+x^2-\dots$, so we have

$\frac{1}{{\left(1+r\right)}^2}=1-2r+3r^2+O\left(r^3\right) ,\frac{1}{4+3r}=\frac{1}{4}\left(1-\frac{3}{4}r+\frac{9}{16}{r}^2+O\left(r^3\right) \right),$

etc. These expressions can be substituted back to Equations (A3) and (A4) to obtain Equation (A14) again. In the same manner, we proved the second-order property of algorithms S3…S5, which is not detailed here, being even more trivial than the case of S1 and S2.  □

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. The stencil of the odd-even hopscotch structures for two nodes. (a) The original OEH algorithm. Green arrows represent first-stage operations and data usage while red arrows are for the second stage. (b) The new, shifted OEH algorithm. Yellow, green, purple, brown, and blue arrows indicate operations at Stages 1, 2, 3, 4 and 5, respectively.

Figure 2. Arrangement of the generalized variables. The double-line, red arrows symbolize conductive (heat) transport through the resistances Rij. The blue line symbolizes thermal isolation at the boundaries of the system. We emphasize that, although topologically the grid should be rectangular to keep the explicit nature of the OEH-type methods, the geometrical shape of the cells is not necessarily rectangular.

Figure 3. Energy errors as a function of the effective time step size for the first (moderately stiff) system, in the case of the original OEH method (OEH REF), the original one-stage CNe and UPFD method, the new algorithms S1–S5, and the Crank–Nicholson method. For the detailed description of the algorithms, see Section 2 and Section 3.2.

Figure 4. Errors as a function of the running time for the first (moderately stiff) system, in the case of the original OEH method (OEH REF), one-stage CNe and UPFD methods, the Dufort–Frankel method, the new Algorithms S1–S5, and the Crank–Nicolson with two implementations and different MATLAB routines.

Figure 5. Energy errors as a function of the running time for the first (moderately stiff) system, in the case of the original OEH method (OEH REF), one-stage CNe and UPFD methods, the Dufort–Frankel method, the new Algorithms S1–S5, and the Crank–Nicolson with two implementations and different MATLAB routines.

Figure 6. Energy errors as a function of the time step size for the second (very stiff) system, in the case of the original OEH method (OEH REF), the original one-stage CNe, UPFD and the Dufort-Frankel methods, the new algorithms S1–S5, and the Crank–Nicholson method.

Figure 7. L∞ errors as a function of the running time for the second (very stiff) system, in the case of the original OEH method (OEH REF), one-stage CNe and UPFD methods, the Dufort–Frankel method, the new algorithms S1–S5, and the Crank–Nicolson with two implementations and different MATLAB routines.

Figure 8. The L∞ errors as a function of hEFF for the u2 solutions of the heat equation for α=1.

Figure 9. The graphs of the second solution, where u0 is the initial function u2(x,  t0) and uref is the analytical solution at the final time, while S3 and S5 are the corresponding numerical solutions served by the S3 (¼, ½, ½, ½, ¾) and the S5 (0, ½, ½, C, 1) algorithms for h=2×10−4.

Figure 10. The L∞ errors as a function of hEFF for the numerical solutions of Fisher’s equation for α=1 and β=2.

Figure 11. The graphs of the eigenvalues of H vs. mesh ratio functions for the S4 (0, ½, ½, ½, 1) algorithm for N = 20. (a) r∈[0,  10] (b) r∈[10,  1000].

Table 1

Comparison of different algorithms for the moderately stiff system of 10,000 cells.

Numerical Method	Running Time (s)	$\mathbf{Error}({\mathit{L}}_{\mathit{\infty }})$	$\mathbf{Error}({\mathit{L}}_1)$	$\mathbf{Energy} \mathbf{Error}$
ode15s, $\mathrm{Tol}={10}^3$	$3.97\times {10}^2$	$1.3\times {10}^{-2}$	$1.1\times {10}^{-3}$	$5.62\times {10}^1$
ode23s, $\mathrm{Tol}={10}^3$	$4.346\times {10}^3$	$4.2\times {10}^{-4}$	$3.0\times {10}^{-5}$	$1.5\times {10}^{-1}$
ode23t, $\mathrm{Tol}={10}^{-8}$	$8.49\times {10}^2$	$2.9\times {10}^{-7}$	$2.0\times {10}^{-8}$	$1.0\times {10}^{-4}$
ode23tb, $\mathrm{Tol}={10}^2$	$4.28\times {10}^2$	$4.1\times {10}^{-4}$	$2.9\times {10}^{-5}$	$1.4\times {10}^{-4}$
ode45, $\mathrm{Tol}={10}^{-1}$	$2.1\times {10}^1$	$3.3\times {10}^{-3}$	$6.5\times {10}^{-5}$	$2.7\times {10}^{-3}$
ode23, $\mathrm{Tol}={10}^{-6}$	$2.7\times {10}^1$	$3.7\times {10}^{-7}$	$9.6\times {10}^{-9}$	$4.8\times {10}^{-5}$
ode113, $\mathrm{Tol}={10}^{-6}$	$1.91\times {10}^1$	$6.7\times {10}^{-7}$	$4.2\times {10}^{-10}$	$1.9\times {10}^{-6}$
S1, $h=1.25\times {10}^{-4}$	$1.97\times {10}^{-1}$	$9.06\times {10}^{-6}$	$2.63\times {10}^{-7}$	$2.56\times {10}^{-3}$
S2, $h=1.25\times {10}^{-3}$	$2.02\times {10}^{-2}$	$3.39\times {10}^{-4}$	$6.93\times {10}^{-6}$	$5.08\times {10}^{-2}$
S3, $h=2.5\times {10}^{-4}$	$1.01\times {10}^{-1}$	$1.88\times {10}^{-5}$	$3.64\times {10}^{-7}$	$3.44\times {10}^{-3}$
S4, $h=5\times {10}^{-4}$	$5.03\times {10}^{-2}$	$1.06\times {10}^{-4}$	$1.07\times {10}^{-6}$	$1.42\times {10}^{-3}$
S5, $h=2.5\times {10}^{-5}$	$9.75\times {10}^{-1}$	$2.62\times {10}^{-7}$	$4.44\times {10}^{-9}$	$3.15\times {10}^{-5}$
CrN inv $h=2.5\times {10}^{-3}$	$1.92\times {10}^1$	$5.56\times {10}^{-5}$	$3.93\times {10}^{-6}$	$1.73\times {10}^{-2}$
CrN lins $h=2.5\times {10}^{-3}$	$1.35\times {10}^2$	$5.56\times {10}^{-5}$	$3.93\times {10}^{-6}$	$1.73\times {10}^{-2}$

Table 2

Comparison of different algorithms for the very stiff system of 10,000 cells.

Numerical Method	Running Time (s)	$\mathbf{Error}({\mathit{L}}_{\mathit{\infty }})$	$\mathbf{Error}({\mathit{L}}_1)$	$\mathbf{Energy} \mathbf{Error}$
ode15s, $\mathrm{Tol}={10}^{-6}$	$6.8\times {10}^2$	$4.1\times {10}^{-7}$	$1.5\times {10}^{-8}$	$7.5\times {10}^{-5}$
ode23s, $\mathrm{Tol}={10}^3$	$5.694\times {10}^3$	$4.7\times {10}^{-4}$	$2.4\times {10}^{-4}$	$1.2\times {10}^{-1}$
ode23t, $\mathrm{Tol}={10}^3$	$3.1\times {10}^2$	$8.1\times {10}^{-2}$	$2.1\times {10}^{-3}$	$1.06\times {10}^1$
ode23tb, $\mathrm{Tol}={10}^{-8}$	$2.037\times {10}^3$	$2.3\times {10}^{-7}$	$1.2\times {10}^{-8}$	$5.8\times {10}^{-5}$
ode45, $\mathrm{Tol}={10}^0$	$9.480\times {10}^3$	$8.1\times {10}^{-2}$	$1.5\times {10}^{-5}$	$7.0\times {10}^{-2}$
ode23, $\mathrm{Tol}={10}^{-6}$	$5.317\times {10}^3$	$1.2\times {10}^{-6}$	$2.3\times {10}^{-10}$	$1.1\times {10}^{-6}$
ode113, $\mathrm{Tol}={10}^{-2}$	$6.046\times {10}^3$	$8.9\times {10}^{-4}$	$1.7\times {10}^{-7}$	$7.7\times {10}^{-4}$
S1, $h=1.25\times {10}^{-4}$	$1.98\times {10}^{-1}$	$8.46\times {10}^{-2}$	$4.55\times {10}^{-4}$	$6.72\times {10}^0$
S2, $h=5.0\times {10}^{-6}$	$4.17\times {10}^0$	$4.81\times {10}^{-4}$	$3.69\times {10}^{-6}$	$6.65\times {10}^{-2}$
S3, $h=2.5\times {10}^{-6}$	$9.85\times {10}^0$	$1.99\times {10}^{-4}$	$7.65\times {10}^{-7}$	$1.31\times {10}^{-2}$
S4, $h=1.25\times {10}^{-4}$	$1.95\times {10}^{-1}$	$3.28\times {10}^{-3}$	$8.88\times {10}^{-6}$	$2.68\times {10}^{-3}$
S5, $h=5\times {10}^{-7}$	$4.95\times {10}^1$	$1.55\times {10}^{-6}$	$8.71\times {10}^{-9}$	$1.69\times {10}^{-4}$
CrN inv $h=5\times {10}^{-4}$	$2.66\times {10}^1$	$1.15\times {10}^{-1}$	$2.42\times {10}^{-4}$	$1.26\times {10}^{-2}$
CrN lins $h=5\times {10}^{-4}$	$7.19\times {10}^2$	$1.15\times {10}^{-1}$	$2.42\times {10}^{-4}$	$1.26\times {10}^{-2}$