Consistency and Convergence Properties of 20 Recent and Old Numerical Schemes for the Diffusion Equation

Abstract

We collected 20 explicit and stable numerical algorithms for the one-dimensional transient diffusion equation and analytically examined their consistency and convergence properties. Most of the methods used have been constructed recently and their truncation errors are given in this paper for the first time. The truncation errors contain the ratio of the time and space steps; thus, the algorithms are conditionally consistent. We performed six numerical tests to compare their performance and try to explain the observed accuracies based on the truncation errors. In one of the experiments, the diffusion coefficient is supposed to change strongly in time, where a nontrivial analytical solution containing the Kummer function was successfully reproduced.

1. Introduction

Fickian diffusion or Fourier-type heat conduction is one of the most fundamental mass- or energy-transport processes. In the simplest case, it can be described by a single linear partial differential equation (PDE) of space and time:

$$\frac{\partial u\left(x,t\right)}{\partial t}=\alpha \frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}$$

(1)

where $x,\hspace{0.17em}\hspace{0.17em}t\in ℝ$, $u=u\left(x,t\right)$ is the concentration (temperature in case of heat conduction) as a function of space and time, and α is the constant diffusion coefficient. The generalizations of this equation, such as the advection (or drift)–diffusion equation and different kinds of advection–diffusion–reaction equations [1], can model the transport and propagation of various particles in physical, chemical and biological systems [2]. Moreover, very similar, albeit more complicated, PDEs are used to model fluid flow through porous media [3], such as moisture [4], ground water or crude oil in underground reservoirs [5].

It is well-known that countless numerical methods have been proposed to solve Equation (1) and its generalizations. Most of them belong to the broad family of finite difference schemes (FDM) [6,7], which often means a kind of method of lines [1], where the spatial variables are discretized first to obtain an ordinary differential equation (ODE) system and then an ODE solver is employed. Nevertheless, finite element methods (FEM) are also used [8], especially if the geometry is complex. Most of these methods can be considered either as an explicit or an implicit method, and sometimes they have a mixed character [9]. If the physical properties, i.e., the coefficients of the equation, are highly non-uniform in space, the eigenvalues of the system matrix likely have a range of several orders of magnitude, which implies that one has a highly stiff problem. In this case, the so-called Courant–Friedrichs–Lewy (CFL) limit (which refers to conventional explicit methods, e.g., the Runge–Kutta or Adams–Bashforth schemes) is very small. This means that these algorithms are unstable unless the time step size is smaller than this low threshold, and therefore the simulation can be unacceptably slow.

Since implicit methods have much better stability properties, they are frequently used for these equations [10,11,12,13,14,15,16]. Nevertheless, they involve the solution of a system of algebraic equations at each time step whose parallelization is much less straightforward than that of explicit methods. Thus, the calculations can be time-consuming for very-large-sized and non-tridiagonal matrices, which is the case when the number of space dimensions is larger than one. In these cases, even the simplest explicit Euler time integration can be much faster than the implicit one [10,11,17]. Regrettably, the formerly rapid increase of the clock frequencies of processors disappeared about two decades ago, which reinforced the progression towards increasing parallelism in high-performance computing [18,19].

These are the main reasons why there is a growing interest in explicit and unconditionally stable algorithms [20,21,22,23,24,25,26]. In the last few years, our research group has developed a couple of new algorithms of this kind. In our original papers [14,27,28,29,30,31,32,33], we examined our new methods analytically and proved that they are indeed unconditionally stable for the linear heat equation. However, when we examined the consistency, we always found that they were only conditionally consistent, in accordance with the literature [11,12,25,34] (p. 120). The goal of this paper is to collect a large number of these methods, give their truncation error, and draw some consequences about the accuracy of the methods. Note that in textbooks and monographs, typically the standard methods, such as the explicit and implicit Euler, as well as the Crank–Nicolson method, are analyzed [35,36]. Richtmyer and Morton [12,13,37], as well as Mitchell and Griffiths [13,14,38], give the principal terms of the truncation error of the Dufort–Frankel method, but none of them do the same for the odd–even hopscotch method, albeit the latter book devotes four pages to that method. We emphasize that the truncation errors of some of the methods have never been calculated before and this might constitute the main novelty of our paper. In the case of our methods, the convergence properties were analyzed previously not via the truncation errors, but in a different way. As we will explain in detail later, this means that the numerical solutions given by the methods were compared to the known analytical solution of the ODE system obtained after the space discretization of the PDE (1). If possible, we will provide both kinds of errors here and compare them. We will also perform some numerical tests aiming to underpin the theoretical results. We note that it was demonstrated that these methods can be used for extremely stiff systems in more than one space dimension. They can provide fairly accurate results much faster than the widely used Runge–Kutta algorithms [14,15,39] or the professionally optimized MATLAB ‘ode’ solvers. Nevertheless, in the current paper, the purpose of the numerical examples is not to test them under extreme conditions, but to illustrate the theoretical results and see which method is more competitive for the studied Equation (1). We stress that unconditional stability is not the rule but rather the exception in the huge crowd of explicit methods. It is well known, for example, that no explicit Runge–Kutta method can be A-stable [40] (p. 60). Therefore, unconditionally stable explicit methods are rarely tested against one another, and performing tests with 20 of them is absolutely unique.

The rest of the paper is organized as follows. In Section 2 we briefly present the numerical algorithms. In Section 3 we describe the different methods for the investigation of consistency and convergence, then presents the calculated error terms. The numerical experiments are presented in Section 4, while the conclusions are summarized in Section 5.

2. The Analyzed Numerical Methods

We considered the 1D interval $x\in \left[x_0\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{x}_N\right]$ with length $L=\hspace{0.17em}{x}_N-x_0$ and constructed an equidistant spatial grid with the node coordinates $x_0\hspace{0.17em},\hspace{0.17em} x_1\hspace{0.17em},\hspace{0.17em} \dots ,\hspace{0.17em}\hspace{0.17em}{x}_N$, so $x_j=x_{j-1}+\Delta x ,\hspace{0.17em} \hspace{0.17em}j=1,\dots ,N,\hspace{0.33em}\hspace{0.17em}\Delta x=L/N\hspace{0.17em}$. The time variable is also discretized, so if $t\in \left[t^0,\hspace{0.17em} t^{\mathrm{fin}}\right]$, then $t^j=t^0+j\Delta t\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}j=1,\dots ,\hspace{0.17em} T\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}T\Delta t=t^{\mathrm{fin}}-t^0$. The mesh-ratio is defined as usual, $r=\frac{\alpha \Delta t}{\Delta x^2}$. In this work, only the simplest case of Equation (1) (one space-dimensional equidistant mesh and constant α coefficient) was examined, but we give the original publications where more details can be found about the application of the methods for general meshes as well. Let us start with a brief presentation of the defining formulas of the methods.

• The UPFD method is proposed by Chen-Charpentier and Kojouharov [41] for the linear diffusion–advection–reaction PDE. It is a non-standard combination of the explicit and implicit Euler discretizations, where only the actual node is considered implicitly and the neighbors explicitly. In the case of Equation (1) it reads:

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

  (2)

The fully explicit form is:

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

(3)

2.

The first scheme we invented is the so-called constant neighbor (CNe) algorithm [32,42], which has the following formula:

$$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),$$

(4)

3.

The CpC method [31] is the organization of the CNe scheme into a two-stage algorithm. The first stage is a fractional time step with length $p\Delta t$ with the CNe formula:

$$u_i^{\mathrm{Cp}}=u_i^n{e}^{-2pr}+\frac{u_{i-1}^n+u_{i+1}^n}{2}\left(1-e^{-2pr}\right).$$

These new and the old values are combined linearly to obtain the predictor values:

$$u_i^{\mathrm{pred}}=\left(1-\frac{1}{2p}\right)u_i^n+\frac{1}{2p}{u}_i^{\mathrm{Cp}}.$$

These are used at the second stage, which is a full-time step size corrector step; thus, the final values are:

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

(5)

4.

CpCC: We now make an attempt to improve the accuracy of the CpC method by iterating one more time using the values obtained by the CpC method as predictors. So, after executing Formula (5) in the CpC method, we set $u_i^{\mathrm{pred}}=u_i^{n+1}$ and then calculate Formula (5) again.

5.

The two-stage linear-neighbor (LNe or Lne2) method [32] starts with taking a full-size time step with the CNe method to obtain the predictor values $u_{\mathrm{i}}^{\mathrm{pred}}$, which are valid at the end of the current time step. Using these values, a quantity can be introduced:

$$s_i=u_{i-1}^{\mathrm{pred}}+u_{i+1}^{\mathrm{pred}}-u_{i-1}^{\mathrm{n}}-u_{i+1}^{\mathrm{n}},$$

which is proportional to the aggregated or effective slopes of the neighbors of node i. Now the corrector values of the two-stage LNe method are:

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

(6)

6.

The LNe3 method is a three-stage algorithm [32] whose first two stages are the same as the LNe schemes. At the end of its second stage, based on the corrector values in Equation (6), one can set $u_i^{\mathrm{pred}}=u_i^{n+1}$, recalculate the slopes $s_i$ and then repeat (6) to obtain new corrector results. This procedure gives a three-stage scheme altogether. This algorithm is still second-order, but more accurate than LNe2.

7.

The LNe4 method is a four-stage algorithm, which is obtained by repeating the procedure explained in the previous point after the calculations of the LNe3 method.

8.

The CLL method [33] is very similar to the LNe3 method, but it applies fractional time steps at the first and second stages with step sizes $\Delta t_1=p\Delta t\hspace{0.17em},\hspace{0.33em} \Delta t_2=p_2\Delta t$. Due to these, it achieves third-order convergence in the time step size, but only if $p_2=\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$. If the length of the second stage is not fixed to this value, then low-order expressions with coefficient $3p_2-2$ appear in the truncation errors, which make the method second-order only. In this paper, we fixed $p_2=\raisebox{1ex}{$2$}\!\left/ \!\raisebox{-1ex}{$3$}\right.$ and keep p as a free parameter. Therefore, the first stage applies the following formula:

$$u_i^{\mathrm{C}}=u_i^n\cdot e^{-2pr}+\frac{u_{i-1}^n+u_{i+1}^n}{2}\left(1-e^{-2pr}\right)$$

In the second stage, a formula similar to (6) is used, but with a reduced time step size:

$$u_i^{\mathrm{CL}}=u_i^n{e}^{-4r/3}+\frac{u_{i-1}^n+u_{i+1}^n}{2}\left(1-e^{-4r/3}\right)+\frac{s_i^1}{2p}\left(\frac{2}{3}-\frac{1-e^{-4r/3}}{2r}\right),$$

(7)

where $s_i^1=u_{i-1}^C+u_{i+1}^{\mathrm{C}}-u_{i-1}^n-u_{i+1}^n$. In the third stage, a full time step is taken:

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

(8)

where $s_i^2=u_{i-1}^{\mathrm{CL}}+u_{i+1}^{\mathrm{CL}}-u_{i-1}^n-u_{i+1}^n$.

9.

The pseudo-implicit (PI) two-stage method was developed in our previous publication [14] (Algorithm 5 there) for the conduction–convection–radiation equation. Here, we apply it only to the pure diffusion Equation (1) with parameter $\lambda =1$, which means that a half time step is taken to obtain the predictor values and then a full time step for the corrector values. The formulas are the following:

Stage 1: $u_i^{\mathrm{pred}}=\frac{u_i^n+\raisebox{1ex}{$r$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\left(u_{i-1}^n+u_{i+1}^n\right)}{1+r}$

Stage 2: $u_i^{n+1}=\frac{\left(1-r\right)u_i^n+r\left(u_{i-1}^{\mathrm{pred}}+u_{i+1}^{\mathrm{pred}}\right)}{1+r}$,

Note that the trick of the implicit and explicit treatment of the actual node and its neighbors is the same as in the UPFD method, so this PI method can actually be considered as a generalization of the UPFD scheme.

10.

The alternating direction explicit (ADE) scheme is a known [26,43] but non-conventional method for which the condition of consistency is also known. We include it here for comparison purposes. In a one-dimensional equidistant mesh, one splits the calculation, i.e., first sweeps the mesh from the left to right, and then vice versa. In the case of Dirichlet boundary conditions at nodes 0 and N, one sets:

$$p_i^n=u_i^n\hspace{0.17em},\hspace{0.33em}i=1,\dots N\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{p}_0^{n+1}=u_0^{n+1} \mathrm{and} q_i^n=u_i^n\hspace{0.17em},\hspace{0.33em}i=N-1,\mathrm{...1}\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{q}_N^{n+1}=u_N^{n+1}.$$

Then, the following equations are solved from left to right and from right to left, respectively:

$$\frac{p_i^{n+1}-p_i^n}{\Delta t}=\frac{\alpha }{\Delta x^2}\left(p_{i-1}^{n+1}-p_i^{n+1}-p_i^n+p_{i+1}^n\right) \mathrm{and} \frac{q_i^{n+1}-q_i^n}{\Delta t}=\frac{\alpha }{\Delta x^2}\left(q_{i-1}^n-q_i^{n+1}-q_i^n+q_{i+1}^{n+1}\right).$$

(9)

The explicit expressions are the following:

$$p_i^{n+1}=\frac{\left(1-r\right)p_i^n+r\left(p_{i-1}^{n+1}+p_{i+1}^n\right)}{1+r} \mathrm{and} q_i^{n+1}=\frac{\left(1-r\right)q_i^n+r\left(q_{i-1}^n+q_{i+1}^{n+1}\right)}{1+r}.$$

(10)

The final values are the simple averages of the two half-sided terms: $u_i^{n+1}=\left(p_i^{n+1}+q_i^{n+1}\right)/2$. If one expresses the final value of u, one obtains:

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

which is actually the same formula as in the Crank–Nicolson method, though the solution of this system of equations is completely different in the two methods. We note that for non-uniform meshes, the ADE method loses its fully explicit character and matrix calculations will be necessary.

11.

The Dufort–Frankel (DF) method [44] (p. 313) is the textbook example of explicit and unconditionally stable methods. It is a one-stage but two-step algorithm with the formula:

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

This method is not a self-starter, and thus $u_i^1$ has to be calculated from $u_i^0$ using another method. We use two half-sized UPFD time step for this purpose.

For the application of odd–even hopscotch methods, the space must be discretized by a special, so-called bipartite grid. The cells are labelled as odd and even, and all the nearest neighbors of the odd nodes are even and vice versa. The spatial and temporal structure of the examined algorithms in space and time are presented in Figure 1. Only one odd and one even cell for each method is shown in the figure, where the time flows from the top to the bottom. The stages are symbolized by colored rectangles, while the repeating units of blocks are surrounded by dashed blue lines. In the case of the shifted-hopscotch structure (b), for example, this unit consists of two half and three full time steps which altogether span two full time steps for each node. First, a half-sized time step (symbolized by a green box with the number ‘1’ inside in Figure 1b) is taken for the odd cells using the initial values, then full time steps are taken strictly alternately for the even and odd cells until the end of the last time step (orange boxes), which should be halved for odd cells to reach exactly the same final time as the even nodes do. The main point is that when a new value of u_i is calculated, the latest values of the neighbors $u_{i\pm 1}$ must always be used, which ensures stability and quite fast convergence at the same time.

12.

The original odd–even hopscotch (OOEH) algorithm has been used for half a century [45]. The structure of this algorithm is shown in Figure 1a. It uses the usual FTCS formula (based on explicit Euler time discretization) at the first stage (labels ‘1’ in the orange boxes in the figure) and the backward time central space (BTCS) formula (implicit Euler time discretization, labels ‘2’ in the figure) at the second stage, which are the following:

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

implicit Euler: $u_i^{n+1}=\frac{u_i^n+r\left(u_{i-1}^{n+1}+u_{i+1}^{n+1}\right)}{1+2r}$.

13.

The RH, or reversed (odd–even) hopscotch scheme [28], applies the same structure as the OOEH method, but the order of the formulas are interchanged, so the implicit formula comes first followed by the explicit one. One might think that this algorithm cannot be explicit, since the new values of the neighbors are not known when first-stage calculations start. To resolve this, the implicit formula is applied with the UPFD trick, which means the neighbors are treated not implicitly, but explicitly; thus, the first-stage formulas are:

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

(11)

Since, as we mentioned, the latest values of the neighbors are always used, it is enough to change the formulas of the first and second stages in the code of the OOEH to obtain the code of the RH algorithm. We demonstrated that this RH algorithm has much smaller errors in the case of extremely stiff systems than the OOEH method.

14.

The OEH-CNe (odd–even hopscotch with CNe) applies the same structure as the OOEH method again, but uses only the CNe formulas (4) appropriately.

15.

The shifted-hopscotch (SH) algorithm [29] uses the theta formulas:

$$u_i^{\mathrm{new}}=\frac{\left(1-2r\theta \right)u_i^{\mathrm{old}}+r\left(u_{i-1}^{\mathrm{latest}}+u_{i+1}^{\mathrm{latest}}\right)}{1+2r\left(1-\theta \right)},$$

(12)

where the ‘old’ and ‘latest’ labels refer to the appropriate time levels, which can be seen in Figure 1b. We always use the $\theta$ values which are proven to yield the most effective version (S4 algorithm in [29]). More concretely, $\theta =0$ is used at the first stage and $\theta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ in the second, third and fourth, while $\theta =1$ is used at the fifth stage. When the length of the stage is halved, $\Delta t$ and thus r must be divided by 2 as well in (12).

16.

The shifted-hopscotch–CNe (SH-CNe) algorithm [29] applies the same space-time structure as the SH method, but uses only the CNe formulas (4) appropriately. For example, the calculation at the first stage is the following:

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

(13)

17.

The asymmetric-hopscotch (ASH) method [46] is similar to the SH one, but with only three stages. As is shown in Figure 1c, there are two half and one full time step size stages with the $\theta$ formula, which together span one full time step for all nodes. Based on numerical experiments, the algorithms have optimal performance if $\theta =0$ is used at the first, $\theta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ at the second, and $\theta =1$ at the third stage (A1 algorithm in [46]).

18.

The asymmetric-hopscotch–CNe (ASH-CNe) algorithm applies the same structure as the SH method, but uses the CNe Formulas (13) and (4) appropriately.

19.

The next method is the leapfrog-hopscotch (LH) method [28]; see Figure 1d. As we will see, this is a very powerful method if one uses $\theta =0$ at the first stage and $\theta =\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ in all other stages (L2 algorithm in [30]).

20.

The leapfrog-hopscotch–CNe (LH-CNe) method [30] is obtained if one replaces the θ formulas in the previous LH algorithm with the CNe formula in each stage of the LH structure (Figure 1d), appropriately.

The algorithms and some of their properties are summarized in

Table 1. List of the methods and their abbreviations.

	Abbrev.	Name of the Method	Recent	Known LTE	O	Conv. Comb.
1.	UPFD	Uncond. positive finite difference	-	+	1	+
2.	CNe	Constant neighbor	+	+	1	+
3.	CpC	Two-stage iterated CNe	+	-	2	+
4.	CpCC	Three-stage iterated CNe	new	-	1	+
5.	LNe	Two-stage Linear neighbor	+	+	2	+
6.	LNe3	Three-stage iterated CNe	+	-	2	+
7.	LNe4	Four-stage iterated CNe	+	-	2	+
8.	CLL	Const.-Lin.-Lin. neighbor fract. step	+	+	3	-
9.	PI	Pseudo-implicit	+	-	2	-
10	ADE	Alternating direction explicit	-	+	2	-
11	DF	Dufort-Frankel	-	+	2	-
12	OOEH	Original odd–even hopscotch	-	?	2	-
13	RH	Reversed hopscotch	+	-	2	-
14	OEH-CNe	OEH structure with CNe formulas	+	-	2	+
15	SH	Shifted hopscotch	+	-	2	-
16	SH-CNe	SH structure with CNe formulas	+	-	2	+
17	ASH	Asymmetric hopscotch	+	-	2	-
18	ASH-CNe	ASH structure with CNe formulas	+	-	2	+
19	LH	Leapfrog-hopscotch	+	-	2	-
20	LH-CNe	LH structure with CNe formulas	+	-	2	+

; for example, the order of temporal convergence, marked by ‘O’. All of the twenty algorithms are unconditionally stable for the linear diffusion equation, i.e., the previously mentioned CFL limit is not valid for them. Some of them have the additional, stronger feature that the new $u_i^{n+1}$ values are the convex combination of the initial $u_j^0$ values, which is indicated in the last column of the table. It implies that the maximum and minimum principles [34] (p. 87) are always fulfilled and unphysical oscillations are never produced, even for very large time step sizes. On the other hand, this favorable property restricts the speed of the convergence of these methods, so they are often the least accurate for small and medium time step sizes, as we will observe later. It means they significantly underestimate the speed of the diffusion process, especially for large and medium time step sizes.

All methods containing the CNe or LNe formula, as well as the LH, RH, SH, ASH and the PI methods, have been recently constructed by our research group, which is indicated after the name of the method in the Table. In our original papers, some analytical proofs, verifications using analytical solutions and numerical tests to compare performances are presented. However, in most cases, the truncation error has not been calculated (marked by a ‘-’ sign in the “known LTE” column).

3. Analytical Results

The most widespread definition of the truncation error is the difference between the discretized and the exact equation [47] (p. 31). For example, in the case of Equation (1) and the forward time central space (FTCS) method,

$$\tau =D_t^{+}u-D_x^{2}u,$$

(14)

where

$$D_t^{+}u=\frac{u\left(x_i^n\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{t}^n+\Delta t\right)-u\left(x_i^n\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{t}^n\right)}{\Delta t} \mathrm{and} D_x^{2}u=\frac{u\left(x_i^n\hspace{0.17em}-\Delta x,\hspace{0.17em}\hspace{0.17em}{t}^n\right)-2u\left(x_i^n\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{t}^n\right)+u\left(x_i^n\hspace{0.17em}+\Delta x,\hspace{0.17em}\hspace{0.17em}{t}^n\right)}{\Delta x^2}$$

are the first-order forward and second-order central difference operators for the time and the space variable, respectively. The Taylor-series expansions of u are inserted into Equation (14). The function u is supposed to be the exact solution, thus the substitutions $u_t=\alpha u_{\left(2x\right)}$, $u_{\left(2t\right)}={\left(\alpha u_{\left(2x\right)}\right)}_t={\alpha }^2{u}_{\left(4x\right)}$, $u_{\left(3t\right)}={\alpha }^3{u}_{\left(6x\right)}$ and so on can be used. The discretization error of $D_x^2$ is:

$${\tau }_0={\epsilon }_0=-\frac{\alpha }{12}{u}_{(4x)}\Delta x^2-\frac{\alpha }{360}{u}_{(6x)}\Delta x^4-\frac{\alpha }{\mathrm{20,160}}{u}_{(8x)}\Delta x^6+\mathrm{h}. \mathrm{o}. \mathrm{t}.,$$

(15)

The discretization error of the $D_t^{+}$ operator depends only on $\Delta t$, thus the space- and time-dependent terms in the truncation error can be clearly separated. This is not true in our cases, where there are terms containing the ratio of $\Delta t$ and $\Delta x$. Moreover, for multistage methods, an expression with a similar form as (15) cannot be obtained, or the method is not clearly a FDM at all, since the differentiation with respect to time is not approximated by a finite difference operator. Nevertheless, all of the examined methods are explicit, and thus they can almost always be expressed in the fully explicit form $u_i^{n+1}=f\left(u_i^n,\hspace{0.17em} u_{i\pm 1}^n\hspace{0.17em},\dots \right)$, where there are no quantities at the new n+1 time level present in the right-hand side, but the f function has the parameter r. Note that in the numerical codes, these forms are used.

The Taylor series expansion can be performed at this level as well [30]. In this way, one obtains the difference between the numerical and the exact solution. Let us denote the truncation error obtained by choosing this way by ɛ. To present a simple example, in the case of the UPFD method, τ is obtained when the expansions are applied to Equation (2), while ɛ is obtained when it is applied to Equation (3). In the latter case, the ${\left(1+2r\right)}^{-1}$ coefficient must also be expanded and the result will be different from τ. One of the sources of the difference is that when Equation (3) is obtained from Equation (2), we multiplied it by $\Delta t$, which is quite common; see p. 20 in [37]. Therefore, ɛ must be divided by $\Delta t$ to make the two kinds of errors comparable. In fact, after this division, both τ and ɛ give an estimation for the global error [30]. However, we will see that τ and ɛ are usually still different and it is a nontrivial question which one reflects more the properties of the numerical scheme, i.e., which one is more useful for practical purposes.

The term (15) is present in case of all the examined methods, both in τ and in ɛ. We note that the higher order terms (h. o. t.) will be usually omitted for the sake of brevity. The τ and ɛ errors have been calculated in most cases by both the first and the second author completely independently.

If one follows the so-called method of lines, where Equation (1) is discretized only spatially using the central difference formula, one obtains the ODE system:

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

(16)

where ${\overrightarrow{u}}^0$ is the arbitrary initial vector. The elements of M can be given as:

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

while the first and last row are determined by the boundary condition. The exact solution of (16) at the end of a time step is the following:

$${\overrightarrow{u}}^{n+1}={\mathrm{e}}^{M\Delta t}{\overrightarrow{u}}^n=\left(1+M\Delta t+M^2\frac{\Delta t^2}{2}+M^3\frac{\Delta t^3}{3!}+O\left(\Delta t^4\right)\right){\overrightarrow{u}}^n.$$

Here, as well as in the case of the $e^{-2r}$-type terms in the CNe, LNe, etc., methods, the usual series expansion $e^{\pm x}=1\pm x+x^2/2!\hspace{0.33em}\hspace{0.17em}\pm x^3/3!\hspace{0.33em}\hspace{0.17em}+O\left(x^4\right)$ is used up to the appropriate order. After some simple algebraic calculations, we obtained the expressions of a general element of the exact solution after one and two time steps:

$$u_i^{n+1}=\frac{r^2}{2}{u}_{i\pm 2}^n+\left(r-2r^2\right)u_{i\pm 1}^n+\left(1-2r+3r^2\right)u_i^n+O\left(r^3\right),$$

(17)

and

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

(18)

respectively. The simple notation $u_{i\pm j}^n=u_{i+j}^n+u_{i-j}^n\hspace{0.17em},\hspace{0.17em} \hspace{0.17em}i,j\in {\mathbb{N}}^{+}$ has been introduced for brevity. Formulas (17) and (18) are valid if we are far enough from the boundaries. To ensure this, we always examine the middle element of a large enough matrix. A similar expression can be obtained for each numerical method. More precisely, for the hopscotch methods, we have different formulas for the odd and the even nodes. If the smallest repeating unit spans one or two time steps, then Equation (17) or (18) must be used, respectively. The local error is the difference between the analytical and the numerical expressions, which is divided by $\Delta t$ again to obtain the global error, which is denoted by $\delta$. We say that a method is, e.g., second-order if the zeroth- and the first-order terms vanish in $\delta$, and the second-order terms have the lowest order. This kind of investigation has been performed for some of our methods in the original publications, but only in order to prove the order of the method. The first non-vanishing terms are given for the first time here, which is also a novelty in the current paper.

The calculated errors are presented in the Appendix A.

4. Numerical Experiments

In this section, we reproduce nontrivial analytical solutions in one space dimension. The Dirichlet boundary conditions at the two end points of the interval and the initial $u^0\left(x,\hspace{0.17em} t\right)$ function were obtained by simply substituting the appropriate x and t values into the analytical solution. The MATLAB R2019b software was used for all numerical calculations. We calculated the eigenvalues of the M matrix in Equation (16) for this problem. Using these, the stiffness ratio and the largest time step size where the standard FTCS method is stable can be determined as usual [33].

To measure the accuracy, the widely used $L_{\infty }$ error is calculated, which compares the exact value $u_i^{\mathrm{analytic}}$ and the result $u_i^{\mathrm{num}}$ obtained by the actual numerical method at the final time $t^{\mathrm{fin}}$:

$$Error=\underset{1\le i\le N}{\max }\left|u_{\mathrm{i}}^{\mathrm{analytic}}(t^{\mathrm{fin}})-u_{\mathrm{i}}^{\mathrm{num}}(t^{\mathrm{fin}})\right|.$$

(19)

Experiment 1. We are going to reproduce the following reference solution, which is a quite recent [48] analytical solution:

$$u^{\mathrm{analytic}}\left(x,t\right)=\frac{x}{t^{7/2}}\left(1-\frac{x^2}{3\alpha t}+\frac{x^4}{60{\alpha }^2{t}^2}\right)e^{-\frac{x^2}{4\alpha t}}.$$

(20)

We considered the domain with the boundaries $x_0=-5\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{x}_N =5\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{t}^0=0.3\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{t}^{\mathrm{fin}}=0.56$, while $\alpha =1$. We set $\Delta x=0.02$, which implies that the stiffness ratio is $1.01\cdot {10}^5$, while $\Delta t_{\mathrm{MAX}}^{\mathrm{FTCS}}=3.9\cdot {10}^{-4}$. We examined the error defined in (19) as a function of the time step size $\Delta t$. First, the errors were calculated for a very large $\Delta t$, then with decreased time step sizes until small error values were obtained. In Figure 2 and Figure 3, the errors as a function of the time step size are shown in log-log diagrams. Since we used a fixed-space step size and decreased only the time step size, the errors cannot go to zero, but only to the residual error due to space discretization, which is given by (15). This can be seen in the bottom left of the figures. One can notice that although the CLL method is third-order and it is much more accurate than the LNe method, it converges more slowly than the hopscotch-type methods. In accordance with some experiments in our previous papers, the LH method was the most accurate.

Figure 4, Figure 5 and Figure 6 show how the numerical solutions approximate the analytical solution as the time step size is decreased. In the case of the one-stage CNe method and the three-stage CLL methods, the curves are smooth without unphysical oscillations, while the ASH method produces significant oscillations for larger time step sizes. On the other hand, the CNe method converges very slowly and the CLL converges at a medium rate, while the ASH converges much faster. Some investigation of the behaviour of the errors are presented in our previous paper [39], as well for fixed time step size and space-dependent diffusion coefficient.

Experiment 2. The solution of the semi-discretized Equation (16) was also calculated numerically by the ode15s MATLAB routine with a very stringent tolerance. The spatial domain was extended to $x_0=-10\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{x}_N =10\hspace{0.17em}$, but all other circumstances of Experiment 1, including the analytical solution (20) and the spatial step size $\Delta x=0.02$, still hold. In Figure 7, we present the errors for six methods when not only the analytical, but the numerical reference was used for the error calculation in Formula (19). One can see that if the numerical reference is used, the residual errors disappear, and for very small time step sizes, the error of the third-order CLL method goes below the errors of the second-order methods.

Experiment 3. In this subsection, $x_0=-5\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{x}_N =5\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{t}^0=0.3\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{t}^{\mathrm{fin}}=0.5$. The space step size $\Delta x$ varied, but all other parameters of Experiment 1 were preserved. For comparison purposes, we used the very common implicit (Euler) method:

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

which was implemented by matrix inversion.

In our previous publications, we measured the running times and conveyed plenty of examples where the traditional methods were outperformed by the explicit and unconditionally stable ones. Therefore, our main purpose with these standard schemes was not the comparison of the performances, but to compare the behavior of the errors of the 20 conditionally consistent schemes with an unconditionally consistent case, so that we could explain properly the connection between the truncation errors and the observed errors.

First, the space step is considered as a parameter with six different values of $\Delta x$ (five in the case of the implicit method, since it has much larger running times). The error-curves are presented for the first-order UPFD and implicit methods in Figure 8, while Figure 9 contains the results for the second-order LNe and LH methods. The curves of the UPFD, LNe and LH method cross one another, because if the space step size is smaller, the convergence is slower due to the $\Delta t^2\Delta x^{-4}$ type terms in the truncation error. It is easy to see that for a fixed time step size, decreasing space step sizes does not yield better accuracy. In the case of the implicit method, however, the error decreases if either the time or the space step is decreased. In the case of the UPFD and the LNe method, the horizontal distance between the neighboring curves (belonging to a certain $\Delta x$ and its half $\Delta x/2$) is a factor of 4, while it is only a factor of 2 for the LH method. This can be explained by the truncation error terms of $\Delta t\Delta x^{-2}$ of the UPFD and $\Delta t^2\Delta x^{-4}$ of LNe schemes, and, on the other hand, the term $\Delta t^2\Delta x^{-2}$ of the LH scheme. If the exponent in the denominator is twice as large as in the numerator, then $\Delta t$ must be divided by four to counterbalance the division of $\Delta x$ by two in order to keep the error constant.

Experiment 4. We conducted a case study considering $\Delta t\to 0$ when $\Delta t/{(\Delta x)}^2=0.2$ was fixed. All parameters were the same as in Experiment 3. The error curves can be seen in Figure 10 for all the 20 methods. Two groups of the methods can be clearly distinguished. The errors of the first group tended to constant values, since their truncation errors contained terms where the exponent of the space step size in the denominator was twice as large as that of the time step size in the numerator. These algorithms were not convergent under these circumstances. On the other hand, the errors of the second group tended to zero, which means they did not have truncation error terms of $\Delta t^n\Delta x^{-2n}$. The results of this numerical investigation are consistent with the analytically calculated truncation error expressions which were detailed in Section 2. The exceptions are the SH and the ASH methods, for which we were not able to calculate the τ errors, only the ɛ ones (see Equations (A1) and (A2)), which contain the critical $\Delta t^n\Delta x^{-2n}$ terms. In these cases, the ɛ errors do not properly reflect the observed properties of these methods from this point of view; thus, we think that further investigations would be necessary.

Experiment 5. Here, everything was the same as in Experiment 1, with the exception that the diffusion coefficient and the final time were increased to $\alpha =2$ and $\hspace{0.17em}{t}^{\mathrm{fin}}=5$, while $\Delta x=0.01$. The error-curves are presented in Figure 11 and Figure 12 for the slowly and the quickly converging methods. We also examined how the maximum errors developed as time elapsed for a given medium-sized time step size. Figure 13 and Figure 14 show that after a sharp increase, the accumulated errors themselves diffuse away and the errors tend slowly toward zero, after some small fluctuations in case of a few quickly converging methods.

Experiment 6. We next considered the generalized diffusion equation where the diffusion constant has a power-law time dependence:

$$\frac{\partial u(x,t)}{\partial t}=\widehat{\alpha }\hspace{0.33em}{t}^n\frac{{\partial }^{2}u(x,t)}{\partial x^2},$$

(21)

where $\widehat{\alpha }$ is considered a constant. The function:

$$u\left(x,t\right)=\frac{x}{t^{b-a}}{e}^{-\frac{(1+n)x^2}{4\widehat{\alpha }{t}^{2b}}}{K}_M\left[\frac{1+n-a}{1+n} \hspace{0.17em} , \hspace{0.17em}\frac{3}{2},\hspace{0.17em} \frac{\left(1+n\right)x^2}{4\widehat{\alpha }{t}^{2b}}\right]\hspace{0.17em}\hspace{0.17em}.$$

(22)

is a very recent [30] analytical solution of Equation (21) if $b=\left(1+n\right)/2$ and $K_M$ is the Kummer–M function, which was calculated here by the hypergeom function of MATLAB. We set the parameter values $\widehat{\alpha }=3.1,\hspace{0.17em}\hspace{0.17em}a=39.2\hspace{0.17em},\hspace{0.17em} n=10.4 ,\hspace{0.17em}\hspace{0.17em}{x}_0=-2\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{x}_N =2\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{t}^0=0.9\hspace{0.17em},\hspace{0.33em}\hspace{0.17em}{t}^{\mathrm{fin}}=0.94$ and $\Delta x=0.01$. Figure 15 and Figure 16 show the maximum errors as a function of the time step size for the slowly and quickly converging schemes, respectively. Figure 17 shows the initial temperatures $u^0$, the analytical solution $u^{anal}$ and some numerical results, which were obtained with four different methods using different time step sizes.

5. Discussion and Summary

We have studied 20 numerical methods for the non-steady-state linear diffusion equation. All of the algorithms are unconditionally stable explicit methods; 16 of them have been recently constructed by our research group. We have calculated the truncation errors for each algorithm, which was unknown until now in several cases. Moreover, we have considered another method of investigation, in which the numerical solution is compared to the concrete analytical solution of the ODE system, which is obtained by the spatial discretization of the PDE. We think that this latter method produces errors which explain the observed numerical order of the algorithms more clearly if the spatial discretization is kept fixed and only the time step size is decreased, especially when the stiffness ratio is high.

We have reproduced a recent analytical solution by the studied 20 numerical methods in six numerical experiments. We have made attempts to explain the observed numerical errors based on the analytically calculated truncation errors. We emphasize again that although these algorithms are only conditionally consistent, they give accurate results for time step sizes orders of magnitude larger than the standard explicit methods, and faster (for large systems) than the implicit methods.