Stability analysis and numerical results for some schemes discretising 2D nonconstant coefficient advection–diffusion equations

Appanah Rao Appadu; Hagos Hailu Gidey

doi:10.1515/phys-2023-0195

Article Open Access

Stability analysis and numerical results for some schemes discretising 2D nonconstant coefficient advection–diffusion equations

Appanah Rao Appadu and Hagos Hailu Gidey

Published/Copyright: March 12, 2024

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 22 Issue 1

Abstract

We solve two numerical experiments described by 2D nonconstant coefficient advection–diffusion equations with specified initial and boundary conditions. Three finite difference methods, namely Lax–Wendroff, Du-Fort–Frankel and a nonstandard finite difference scheme, are derived and used to solve the two problems, whereby only the first problem has an exact solution. Stability analysis is performed to obtain a range of values of the time step size at a fixed spatial step size. We obtain the rate of convergence in space when the three methods are used to solve Problem 1. Computational times of the three algorithms are computed for Problem 1. Results are displayed for the two problems using the three methods at times T = 1.0 and T = 5.0 . The main novelty is the stability analysis, which is not straightforward as we are working with numerical methods discretising 2D nonconstant coefficient advection–diffusion equation where many parameters are involved. The second highlight is to determine the most efficient scheme from the three methods. Third, there are very few published studies on analysis and use of numerical methods to solve nonconstant coefficient advection–diffusion equations, and this is one of the very few rare articles treating such topics.

Keywords: nonconstant coefficient; advection–diffusion; Lax–Wendroff; Du-Fort–Frankel; nonstandard finite difference; stability; rate of convergence

Nomenclature

k: temporal step size
Δ x: spatial step size in the direction of x
Δ y: spatial step size in the direction of y
ξ: amplification factor
D 1 , D 2: coefficients of dissipation
t: time
C i , j n: numerical solution at the grid point ( t n , x i , y j )

1 Introduction

The advection–diffusion equation describes the transport of a quantity due to two processes: advection and diffusion. This equation is widely used in various fields of science and engineering such as fluid dynamics, heat transfer, chemical engineering, and atmospheric sciences [1–3]. The numerical solution of advection–diffusion equation is generally a challenging problem due to its nature as it consists of advection and diffusion terms [4].

Finite difference method is one of the most popular classes of numerical methods used to solve advection–diffusion equations. Finite difference methods for advection–diffusion equations with uniform flow and constant coefficients have been extensively studied and developed [5–9].

One of the most common ways of measuring the relative merit of a numerical scheme for advection is to consider the scheme’s dispersion and dissipation [10]. The pioneering work of the theoretical study of finite difference methods was made by Courant et al. [11]. Von Neumann and Richtmyer [12] developed Fourier analysis method of finite difference schemes. Kreiss [13] initiated work on the dispersion and dissipation of finite difference schemes discretising partial differential equations using Fourier method. The dissipative and dispersive features of the Lax–Wendroff (LW) and the MacCormack schemes discretising the linear and non-linear advection equation are discussed in the study of Winnicki et al. [14]. Hirt [15] and Ru-xun and Zhao-hui [16] worked on the remainder analysis approach of finite difference methods. The relative phase error is a measure of the dispersive character of a numerical method [17]. The relative phase error is a ratio that measures the velocity of the computed waves to that of the physical waves [17].

Appadu [18] used the LW, Crank–Nicolson, and a nonstandard finite difference scheme to solve a one-dimensional advection–diffusion equation with constant coefficients. Two optimisation techniques based on minimisation of dispersion error were implemented to find the optimal value of the time step size when the spatial step size is chosen as h = 0.02 , and this is validated using some numerical experiments. Appadu et al. [8] used three numerical methods to solve two problems described by advection–diffusion equations with specified initial and boundary conditions. Two test problems were considered. The first test problem considered has steep boundary layers near x = 1 , and this is a challenging problem as many schemes are plagued by non-physical oscillation near steep boundaries.

The regions of stability of forward-time central space (FTCS) and LW schemes discretising the 1D advection–diffusion equation given by [19]:

(1) ∂ C ∂ t + a ∂ C ∂ x = α ∂ 2 C ∂ x 2 ,

are c 2 2 ≤ s ≤ 1 2 and 0 < s ≤ 1 − c 2 2 , respectively, where s = α Δ t Δ x 2 and c = a Δ t Δ x . We note that a is the coefficient of advection and α is the coefficient of diffusion.

Hutomo et al. [20] used Du-Fort–Frankel scheme in order to solve some nonconstant coefficient advection–diffusion equations on regular and irregular grids. They considered five numerical experiments, out of which only one of them has exact solution. However, no detailed proof of the stability of the Du-Fort–Frankel is given.

In this work, we derive three methods to solve two problems described by nonconstant coefficient advection–diffusion equation. Exact solution is known for only one of the problems considered. Obtaining the stability of numerical methods for advection–diffusion equation is not straightforward and there were several attempts in the past by various authors and some of the previously published studies [15,21–25]. Complication arises for nonconstant coefficient advection–diffusion equations. In this current work, the amplification factor is a function of x i , y j , ω x , ω y , Δ x , Δ y , k , u i , j , and v i , j and the constants D 1 and D 2 . Two approaches are used in order to obtain the range of values of k for stability. In the first approach, we fix D 1 , D 2 , Δ x , Δ y and select some values of x i ∈ [ 0 , 1 ] and y j ∈ [ 0 , 1 ] . We then obtain 3D plots of the modulus of amplification factor ξ vs ω x ∈ [ − π , π ] vs ω y ∈ [ − π , π ] at some selected value of k , starting with a very small value of k and check if ∣ ξ ∣ ≤ 1 for stability. We then increase k gradually until ∣ ξ ∣ ≤ 1 no longer holds. In the second approach, we use the technique of Hindmarsh et al. [25]. We fix Δ x , Δ y , D 1 , D 2 to obtain the amplification factor ξ in terms of k , x i , y j . Then, we obtain 3D plots of ∣ ξ ∣ vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] and obtain range of values of k for which ∣ ξ ∣ ≤ 1 .

This article is organised as follows. We describe the two numerical experiments in Section 2. In Section 3, we derive the three methods, namely Lax–Wendroff, Du-Fort–Frankel and nonstandard finite difference method (NSFD) methods to solve Problem 1. Section 4 is devoted to derivation and study of stability of the methods to solve Problem 2. Sections 5 and 6 provide the numerical results from Problems 1 and 2, respectively. Section 7 highlights the salient features of this article.

2 Numerical experiments

We solve the 2D nonconstant coefficient advection–diffusion equation given by:

(2) ∂ C ∂ t + ∂ ∂ x ( u C ) + ∂ ∂ y ( v C ) = D 1 ∂ 2 C ∂ x 2 + D 2 ∂ 2 C ∂ y 2 ,

where D 1 and D 2 are the constants and u and v are the nonconstants.

2.1 Problem 1 [20]

We solve Eq. (2), where

(3) u ( x ) = D 1 β − α 2 β + p e − β x ,

(4) v ( y ) = D 2 γ − α 2 γ + q e − γ y ,

with the following parameters:

α = − 0.029 , β = 0.5 , γ = 0.5 , p = 0.05 , q = 0.05 , D 1 = 0.004 , and D 2 = 0.004 .

The exact solution is given by [20]:

(5) C ( x , y , t ) = e α t + β x + γ y ,

and the boundary conditions are deduced from exact solution. Dirichlet boundary conditions are used at the boundaries. The domains are x , y ∈ [ 0 , 1 ] and t ∈ ( 0 , T ] . We choose the spatial step sizes in the x and y directions as Δ x = Δ y = 0.05 . For the domain selected, the range of both u and v is 0.06 to 0.08, and this can be easily obtained using any suitable numerical software such as Maple. We have decided to obtain profiles at short and longer time propagation, and this is why we chose T as 1.0 and 5.0 .

2.2 Problem 2 [20]

In this numerical experiment, we solve Eq. (2) where the nonconstant coefficient advection terms are given by:

u ( x , y ) = 0.01 + 0.005 x − 0.005 y , v ( x , y ) = − 0.01 − 0.005 x + 0.005 y ,

with domain x , y ∈ [ 0 , 1 ] and the constants

D 1 = D 2 = 0.0004 . The initial condition is given by:

C ( x , y , t = 0 ) = 10 , if x = 0.5 , y = 0.5 1 , otherwise .

The boundary conditions are as follows:

C ( x = 0 , y , t ) = C ( x = 1 , y , t ) = 1 , C ( x , y = 0 , t ) = C ( x , y = 1 , t ) = 1 .

There is no exact solution for this problem. We display the results at times T = 1 and T = 5 using Δ x = Δ y = 0.05 .

3 Derivation and stability analysis of the three methods for Problem 1

3.1 LW method

We first obtain the family of explicit and implicit schemes in order to discretise the constant coefficient advection–diffusion equation [18]:

(6) ∂ C ∂ t + u ∂ C ∂ x − D ∂ 2 C ∂ x 2 = 0 .

The family of methods is given by:

C i n + 1 − C i n k + u 1 − ϕ Δ x ( 1 − χ ) C i n − C i − 1 n Δ x + χ C i + 1 n − C i − 1 n 2 Δ x + ϕ Δ x ( 1 − χ ) C i n + 1 − C i − 1 n + 1 Δ x + χ C i + 1 n + 1 − C i − 1 n + 1 2 Δ x − D 1 − ϕ ( Δ x ) 2 ( C i + 1 n − 2 C i n + C i − 1 n ) + ϕ ( Δ x ) 2 ( C i + 1 n + 1 − 2 C i n + 1 + C i − 1 n + 1 ) = 0 ,

where ϕ and χ are the temporal and spatial weighting factors. To obtain the LW scheme, we use ϕ = 0 and χ = 1 − u 2 . Thus, the LW scheme when used to discretise Eq. (6) is [6,26]

C i n + 1 − C i n k + u c C i n − C i − 1 n Δ x + ( 1 − c ) C i + 1 n − C i − 1 n 2 Δ x − D C i + 1 n − 2 C i n + C i − 1 n ( Δ x ) 2 = 0 ,

where c = u k Δ x .

Suppose we have a 1D nonconstant coefficient advection–diffusion equation of the form:

∂ C ∂ t + ∂ ∂ x ( u C ) − D ∂ 2 C ∂ 2 x = 0 ,

where u = u ( x ) and D is a constant. This equation can be written as:

(7) ∂ C ∂ t + u ∂ C ∂ x + ∂ u ∂ x C − D ∂ 2 C ∂ 2 x = 0 .

When Eq. (7) is discretised using the LW scheme, the following discretisations are used:

∂ C ∂ t ≈ C i n + 1 − C i n k , ∂ C ∂ x ≈ c i C i n − C i − 1 n Δ x + ( 1 − c i ) C i + 1 n − C i − 1 n 2 Δ x , ∂ 2 C ∂ x 2 ≈ C i + 1 n − 2 C i n + C i − 1 n ( Δ x ) 2 ,

where c i = u i k Δ x . This gives the following scheme:

(8) C i n + 1 − C i n k + u i c i C i n − C i − 1 n Δ x + ( 1 − c i ) C i + 1 n − C i − 1 n 2 Δ x + ∂ u ∂ x i C i n − D C i + 1 n − 2 C i n + C i − 1 n ( Δ x ) 2 = 0 .

Eq. (2) can be rewritten as follows:

(9) ∂ C ∂ t + ∂ u ∂ x C + u ∂ C ∂ x + ∂ v ∂ y C + v ∂ C ∂ y = D 1 ∂ 2 C ∂ x 2 + D 2 ∂ 2 C ∂ y 2 .

The following approximations are used when Eq. (9) is discretised using LW scheme:

∂ C ∂ t ≈ C i , j n + 1 − C i , j n k , ∂ C ∂ x ≈ u i , j k Δ x C i , j n − C i − 1 , j n Δ x + 1 − u i , j k Δ x C i + 1 , j n − C i − 1 , j n 2 Δ x , ∂ 2 C ∂ x 2 ≈ C i + 1 , j n − 2 C i , j n + C i − 1 , j n ( Δ x ) 2 , ∂ C ∂ y ≈ v i , j k Δ y C i , j n − C i , j − 1 n Δ y + 1 − v i , j k Δ y C i , j + 1 n − C i , j − 1 n 2 Δ y , ∂ 2 C ∂ y 2 ≈ C i , j + 1 n − 2 C i , j n + C i , j − 1 n ( Δ y ) 2 .

The LW scheme when used to solve Problem 1 is given by:

(10) C i , j n + 1 − C i , j n k + p ( − β ) e − β x i C i , j n + q ( − γ ) e − γ y j C i , j n + u i u i k Δ x C i , j n − C i − 1 , j n Δ x + 1 − u i k Δ x C i + 1 , j n − C i − 1 , j n 2 Δ x + v j v j k Δ x C i , j n − C i , j − 1 n Δ y + 1 − v j k Δ y C i , j + 1 n − C i , j − 1 n 2 Δ y = D 1 C i + 1 , j n − 2 C i , j n + C i − 1 , j n ( Δ x ) 2 + D 2 C i , j + 1 n − 2 C i , j n + C i , j − 1 n ( Δ y ) 2 ,

where

(11) u i = D 1 β − α 2 β + p e − β x i ,

(12) v j = D 2 γ − α 2 γ + q e − γ y j .

Plots of u ( x ) vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] and v ( x ) vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] when D 1 = D 2 = 0.004 , α = − 0.029 , β = 0.5 , γ = 0.5 , p = 0.05 , and q = 0.05 are displayed in Figure 1, and we observe that both quantities are non-negative.

$Figure 1 Plots of u ( x ) u\left(x) vs x ∈ [ 0 , 1 ] x\in \left[0,1] vs y ∈ [ 0 , 1 ] y\in \left[0,1] and v ( x ) v\left(x) vs x ∈ [ 0 , 1 ] x\in \left[0,1] vs y ∈ [ 0 , 1 ] y\in \left[0,1] : (a) u u and (b) v v .$

Figure 1

Plots of u ( x ) vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] and v ( x ) vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] : (a) u and (b) v .

Using Eq. (10), we obtain

(13) C i , j n + 1 = C i , j n + p β k e − β x i C i , j n + q γ k e − γ y j C i , j n − u i k u i k Δ x C i , j n − C i − 1 , j n Δ x + 1 − u i k Δ x C i + 1 , j n − C i − 1 , j n 2 Δ x − v j k v j k Δ x C i , j n − C i , j − 1 n Δ y + 1 − v j k Δ y C i , j + 1 n − C i , j − 1 n 2 Δ y + D 1 k ( Δ x ) 2 ( C i + 1 , j n − 2 C i , j n + C i − 1 , j n ) + D 2 k ( Δ y ) 2 ( C i , j + 1 n − 2 C i , j n + C i , j − 1 n ) .

To study the stability, we use the ansatz C i , j n = ξ n e I θ x i Δ x e I θ y j Δ y [27], where θ x and θ y are the wave numbers along x - and y - directions, respectively, and I = − 1 . Let ω x = θ x Δ x and ω y = θ y Δ y . The amplification factor is given by:

(14) ξ = 1 + p β k e − β x i + q γ k e − γ y j − u i k u i k ( Δ x ) 2 ( 1 − e − I ω x ) + 1 − u i k Δ x 2 I sin ( ω x ) 2 Δ x − v j k v j k ( Δ y ) 2 ( 1 − e − I ω y ) + 1 − v j k Δ y 2 I sin ( ω y ) 2 Δ y + D 1 k ( Δ x ) 2 ( e I ω x − 2 + e − I ω x ) + D 2 k ( Δ y ) 2 ( e I ω y − 2 + e − I ω y ) .

We note that for Problem 1, the parameters are [20] as follows:

α = − 0.029 , β = 0.5 , γ = 0.5 , p = 0.05 , q = 0.05 , D 1 = 0.004 , and D 2 = 0.004 .

However, different parameters can be used in future studies. We choose Δ x = Δ y = 0.05 and use the parameters which Hutomo et al. [20] used for Problem 1. This gives the amplification factor as follows:

ξ = 1 + 0.0025 k e − 0.5 x i + 0.025 k e − 0.5 y j + 0.124 k 2 e − 0.5 x i + I sin ( ω x ) e − 0.5 x i [ k − 0.124 k 2 − 0.01 k 2 e − 0.5 x i ] + I sin ( ω y ) e − 0.5 y j [ k − 1.24 k 2 − k 2 e − 0.5 y j ] + I sin ( ω x ) [ 3.82 k − 0.3844 k 2 ] − 0.3844 k 2 e − I ω x + I sin ( ω y ) [ 3.82 k − 0.3844 k 2 ] + 0.01 k 2 ( e − 0.5 x i ) 2 + 1.24 k 2 e − 0.5 y j [ 1 − e − I ω y ] + k 2 ( e − 0.5 y j ) 2 [ 1 − e − I ω y ] + 0.7688 k 2 − 0.124 k 2 e − 0.5 x i e − I ω x − 0.01 k 2 ( e − 0.5 x i ) 2 e − I ω x − 0.3844 k 2 e − I ω y − 6.4 k ,

which can be written as follows:

ξ = 1 − 6.4 k + 0.025 k e − 0.5 x i + 0.025 k e − 0.5 y j − 0.7688 k 2 − 1.24 k 2 e − 0.5 x i − k 2 ( e − 0.5 x i ) 2 + 1.24 k 2 e − 0.5 y j cos ( ω y ) − 1.24 k 2 e − 0.5 y j − k 2 ( e − 0.5 y j ) 2 + k 2 ( e − 0.5 y j ) 2 cos ( ω y ) + 3.2 k cos ( ω x ) + 3.2 k cos ( ω y ) + 1.24 k 2 e − 0.5 x i cos ( ω x ) + k 2 ( e − 0.5 x i ) 2 cos ( ω x ) + 0.3844 k 2 cos ( ω x ) + 0.3844 k 2 cos ( ω y ) − 0.62 I k sin ( ω y ) − I k [ e − 0.5 x i sin ( ω x ) + e − 0.5 y j sin ( ω y ) + 0.62 sin ( ω x ) ] .

We note that ξ consists of five parameters: x i , y j , k , ω x , ω y .

Since x i , y j ∈ [ 0 , 1 ] , we can choose x i = y j = 1 . We obtain 3D plots of ∣ ξ ∣ vs ω x ∈ [ − π , π ] vs ω y ∈ [ − π , π ] at some selected values of k , starting with a very small k , say k = 0.001 , and gradually increasing the value. For stability, we need ∣ ξ ∣ ≤ 1 . We observe that the scheme is stable when k ≤ 0.14 , as depicted in Figures 2 and 3.
We now choose x i = y j = 0.5 and repeat the steps in (i). We observe that the range of values of k for stability is again k ≤ 0.14 , as displayed in Figures 4 and 5.

$Figure 2 3D plots of ∣ ξ ∣ | \xi | vs ω x ∈ [ ‒ π , π ] {\omega }_{x}\in \left[‒\pi ,\pi ] vs ω y ∈ [ ‒ π , π ] {\omega }_{y}\in \left[‒\pi ,\pi ] when x i = y j = 1 {x}_{i}={y}_{j}=1 and using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.001 , 0.01 , 0.1 k=0.001,0.01,0.1 : (a) k = 0.001 k=0.001 , (b) k = 0.01 k=0.01 , and (c) k = 0.1 k=0.1 .$

Figure 2

3D plots of ∣ ξ ∣ vs ω x ∈ [ ‒ π , π ] vs ω y ∈ [ ‒ π , π ] when x i = y j = 1 and using Δ x = Δ y = 0.05 and k = 0.001 , 0.01 , 0.1 : (a) k = 0.001 , (b) k = 0.01 , and (c) k = 0.1 .

$Figure 3 3D plots of ∣ ξ ∣ | \xi | vs ω x ∈ [ ‒ π , π ] {\omega }_{x}\in \left[‒\pi ,\pi ] vs ω y ∈ [ ‒ π , π ] {\omega }_{y}\in \left[‒\pi ,\pi ] when x i = y j = 1 {x}_{i}={y}_{j}=1 and using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.12 , 0.14 k=0.12,0.14 , and 0.15: (a) k = 0.12 k=0.12 , (b) k = 0.14 k=0.14 , and (c) k = 0.15 k=0.15 (showing unstability).$

Figure 3

3D plots of ∣ ξ ∣ vs ω x ∈ [ ‒ π , π ] vs ω y ∈ [ ‒ π , π ] when x i = y j = 1 and using Δ x = Δ y = 0.05 and k = 0.12 , 0.14 , and 0.15: (a) k = 0.12 , (b) k = 0.14 , and (c) k = 0.15 (showing unstability).

$Figure 4 3D plots of ∣ ξ ∣ | \xi | vs ω x ∈ [ ‒ π , π ] {\omega }_{x}\in \left[‒\pi ,\pi ] vs ω y ∈ [ ‒ π , π ] {\omega }_{y}\in \left[‒\pi ,\pi ] when x i = y j = 0.5 {x}_{i}={y}_{j}=0.5 and using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.001 , 0.01 k=0.001,0.01 , and 0.1: (a) k = 0.001 k=0.001 , (b) k = 0.01 k=0.01 , and (c) k = 0.1 k=0.1 .$

Figure 4

3D plots of ∣ ξ ∣ vs ω x ∈ [ ‒ π , π ] vs ω y ∈ [ ‒ π , π ] when x i = y j = 0.5 and using Δ x = Δ y = 0.05 and k = 0.001 , 0.01 , and 0.1: (a) k = 0.001 , (b) k = 0.01 , and (c) k = 0.1 .

$Figure 5 3D plots of ∣ ξ ∣ | \xi | vs ω x ∈ [ ‒ π , π ] {\omega }_{x}\in \left[‒\pi ,\pi ] vs ω y ∈ [ ‒ π , π ] {\omega }_{y}\in \left[‒\pi ,\pi ] when x i = y j = 0.5 {x}_{i}={y}_{j}=0.5 and using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.12 , 0.14 k=0.12,0.14 , and 0.15: (a) k = 0.12 k=0.12 , (b) k = 0.14 k=0.14 , (c) k = 0.15 k=0.15 (showing unstability).$

Figure 5

3D plots of ∣ ξ ∣ vs ω x ∈ [ ‒ π , π ] vs ω y ∈ [ ‒ π , π ] when x i = y j = 0.5 and using Δ x = Δ y = 0.05 and k = 0.12 , 0.14 , and 0.15: (a) k = 0.12 , (b) k = 0.14 , (c) k = 0.15 (showing unstability).

3.2 Du-Fort–Frankel method

Du Fort and Frankel [28] proposed a modification of the centred time centred space (CTCS) scheme to ensure second-order accuracy while improving numerical stability [20]. The idea is to replace the middle term in the numerator of CTCS method with an average value. Hence, we have the following scheme to discretise Eq. (9):

(15) C i , j n + 1 − C i , j n − 1 2 k + u i + 1 − u i − 1 2 Δ x C i , j n + u i C i + 1 , j n − C i − 1 , j n 2 Δ x + v j + 1 − v j − 1 2 Δ y C i , j n + v j C i , j + 1 n − C i , j − 1 n 2 Δ y = D 1 C i + 1 , j n − C i , j n + 1 − C i , j n − 1 + C i − 1 , j n ( Δ x ) 2 + D 2 C i , j + 1 n − C i , j n + 1 − C i , j n − 1 + C i , j − 1 n ( Δ y ) 2 .

On rearrangement, we obtain

(16) C i , j n + 1 = 1 1 + 2 B x + 2 B y { ( 1 − 2 B x − 2 B y ) C i , j n − 1 + u i k Δ x + 2 B x C i − 1 , j n + − u i k Δ x + 2 B x C i + 1 , j n + v j k Δ y + 2 B y C i , j − 1 n + − v j k Δ y + 2 B y C i , j + 1 n − k Δ x ( u i + 1 − u i − 1 ) + k Δ y ( v j + 1 n − v j − 1 ) C i , j n ,

where B x = k D 1 ( Δ x ) 2 and B y = k D 2 ( Δ y ) 2 .

The amplification factor ξ satisfies the following equation:

(17) ξ = 1 1 + 2 B x + 2 B y ( 1 − 2 B x − 2 B y ) ξ − 1 + u i k Δ x + 2 B x e − I ω x + − u i k Δ x + 2 B x e I ω x + v j k Δ y + 2 B y e − I ω y + − v j k Δ y + 2 B y e I ω y − k Δ x ( u i + 1 − u i − 1 ) + k Δ y ( v j + 1 − v j − 1 ) .

We can express u i + 1 − u i − 1 as p e − β x i + 1 − p e − β x i − 1 and v j + 1 n − v j − 1 as q e − γ y j + 1 − q e − γ y j − 1 using Eqs. (11) and (12).

Let x i + 1 , x i , x i − 1 be denoted by x 3 , x 2 , x 1 and y j + 1 , y j , y j − 1 be denoted by y 3 , y 2 , y 1 , respectively.

We choose Δ x = Δ y = 0.05 along with the required parameters for Problem 1. This gives

ξ = 1 1 + 6.4 k { ( 1 − 6.4 k ) ξ − 1 + ( 20 ( 0.031 + 0.05 e − 0.5 x 2 ) k + 3.2 k ) e − I ω x + ( 3.2 k − 20 ( 0.031 + 0.05 e − 0.5 x 2 ) k ) e I ω x + ( 3.2 k + 20 ( 0.031 + 0.05 e − 0.5 y 2 ) k ) e − I ω y + ( 3.2 k − 20 ( 0.031 + 0.05 e − 0.5 x 2 ) k ) e I ω y − 20 k [ ( 0.05 e − 0.5 x 3 − 0.05 e − 0.5 x 1 ) + ( 0.05 e − 0.5 y 3 − 0.05 e − 0.5 y 1 ) ] } .

Since Δ x = Δ y = 0.05 and x , y ∈ [ 0 , 1 ] , we can choose, for instance, among the various possibilities: x 3 = 0.55 , x 2 = 0.5 , x 1 = 0.45 and y 3 = 0.55 , y 2 = 0.5 , y 1 = 0.45 . This then gives the following equation involving the amplification factor, ξ :

(18) ξ ( 1 + 6.4 k ) = ( 1 − 6.4 k ) ξ − 1 + 0.0778819120 k + 4.598800783 k ( e − I ω x + e − I ω y ) + 1.801199217 k ( e I ω x + e I ω y ) .

From there, we check if the scheme is stable for some selected values of k . For stability, we need to have ∣ ξ ∣ ≤ 1 for ω x ∈ [ − π , π ] and ω y ∈ [ − π , π ] at those selected values of k . We start with low value of k and increase k gradually. Let us choose k = 0.01 . We obtain the following quadratic equation:

(19) ξ 2 − 0.8796992482 − 0.04322181187 ( cos ( ω x ) − I sin ( ω x ) ) ξ − 0.01692856407 ( cos ( ω x ) + I sin ( ω x ) ) ξ − 0.04322181187 ( cos ( ω y ) − I sin ( ω y ) ) ξ − 0.01692856407 ( cos ( ω y ) + I sin ( ω y ) ) ξ − 0.0007320318722 ξ = 0 .

Let the two roots be ξ 1 and ξ 2 . We obtain plots of ∣ ξ 1 ∣ vs ω x ∈ [ − π , π ] vs ω y ∈ [ − π , π ] and also ∣ ξ 2 ∣ vs ω x ∈ [ − π , π ] vs ω y ∈ [ − π , π ] in Figure 6.

We now repeat the same procedures from Eqs. (18) to (19) and use different values of k each occasion and check if the scheme is stable at that value of k . We will not present all the figures due to restriction in number of pages possible. We present Table 1, which summarises the range of values of k for stability. For stability, we need k ≤ 0.16 .

Table 1

Checking stability of Du-Fort–Frankel scheme for Problem 1 at some values of k when Δ x = Δ y = 0.05

Value of k	Is scheme stable at that value of k ?
0.001	Yes
0.01	Yes
0.1	Yes
0.16	Yes
0.17	No
0.18	No

$Figure 6 Plots of the modulus of amplification factors vs ω x {\omega }_{x} vs ω y {\omega }_{y} when Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.01 k=0.01 : (a) ∣ ξ 1 ∣ | {\xi }_{1}| vs ω x {\omega }_{x} vs ω y {\omega }_{y} and (b) ∣ ξ 2 ∣ | {\xi }_{2}| vs ω x {\omega }_{x} vs ω y {\omega }_{y} .$

Figure 6

Plots of the modulus of amplification factors vs ω x vs ω y when Δ x = Δ y = 0.05 and k = 0.01 : (a) ∣ ξ 1 ∣ vs ω x vs ω y and (b) ∣ ξ 2 ∣ vs ω x vs ω y .

3.3 NSFD

We refer to Kojouharov and Cheni [29] and Mickens and Washington [30], who constructed a nonstandard finite difference scheme for an advection–diffusion reaction equation of the form:

u t + ν 0 u x = D u x x − λ u − ε u 1 ∕ 3 ,

where ν 0 , D , and λ ε are non-negative parameters.

To construct the NSFD scheme for Eq. (9), we use the following approximations:

∂ C ∂ t ≈ C i , j n + 1 − C i , j n ϕ ( k ) , ∂ C ∂ x ≈ C i , j n − C i − 1 , j n ψ ( Δ x ) , ∂ 2 C ∂ x 2 ≈ C i + 1 , j n − 2 C i , j n + C i − 1 , j n ( ψ ( Δ x ) ) 2 , ∂ C ∂ y ≈ C i , j n − C i , j − 1 n ψ ( Δ y ) , ∂ 2 C ∂ y 2 ≈ C i , j + 1 n − 2 C i , j n + C i , j − 1 n ( ψ ( Δ y ) ) 2 ,

where ψ ( Δ x ) = exp ( Δ x ) − 1 and ϕ ( k ) = exp ( k ) − 1 . This gives the NSFD scheme as:

(20) C i , j n + 1 − C i , j n ϕ ( k ) + ∂ u ∂ x i C i , j n + u i C i , j n − C i − 1 , j n ψ ( Δ x ) + ∂ v ∂ y j C i , j n + v j C i , j n − C i , j − 1 n ψ ( Δ y ) = D 1 C i + 1 , j n − 2 C i , j n + C i − 1 , j n ( ψ ( Δ x ) ) 2 + D 2 C i , j + 1 n − 2 C i , j n + C i , j − 1 n ( ψ ( Δ y ) ) 2 .

We choose the functional relation:

(21) ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 1 2 ,

and we obtain the following scheme:

(22) C i , j n + 1 = 1 + p β ϕ ( k ) e − β x i − u i ϕ ( k ) ψ ( Δ x ) + q γ ϕ ( k ) e − γ y j − v j ϕ ( k ) ψ ( Δ y ) − D 1 − D 2 C i , j n + D 1 2 C i + 1 , j n + D 2 2 C i , j + 1 n + u i ϕ ( k ) ψ ( Δ x ) + D 1 2 C i − 1 , j n + v j ϕ ( k ) ψ ( Δ y ) + D 2 2 C i , j − 1 n .

To run experiments, we use (21), choose Δ x = Δ y = 0.05 , and this gives k = ln ( 1 + 0.5 × 0.0 5 2 ) ≃ 1.249 × 1 0 − 3 .

It is clear that the coefficients of C i + 1 , j n and C i , j + 1 n are non-negative. Also, the coefficients of C i − 1 , j n and C i , j − 1 n are non-negative as we see from Figure 1 that 0.06 ≤ u i ≤ 0.08 and 0.06 ≤ v j ≤ 0.08 . We display the variation of coefficient of C i , j n vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] in Figure 7 and see that it is in the range of 0.9882 to 0.9890 for x ∈ [ 0 , 1 ] and y ∈ [ 0 , 1 ] . We therefore conclude that NSFD replicates the positivity of the continuous model when the functional relation in (21) is used. There is no need to obtain the stability of NSFD using von Neumann stability analysis. We can just check if the conditions for which NSFD replicates positivity of the continuous model.

$Figure 7 Coefficient of C i , j n {C}_{i,j}^{n} vs x ∈ [ 0 , 1 ] x\in \left[0,1] vs y ∈ [ 0 , 1 ] y\in \left[0,1] .$

Figure 7

Coefficient of C i , j n vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] .

4 Derivation and stability analysis of the three methods for Problem 2

4.1 LW scheme

The scheme is given by:

(23) C i , j n + 1 − C i , j n k + 0.01 C i , j n + u i , j u i , j k Δ x C i , j n − C i − 1 , j n Δ x + 1 − u i , j k Δ x C i + 1 , j n − C i − 1 , j n 2 Δ x + v i , j v i , j k Δ y C i , j n − C i , j − 1 n Δ y + 1 − v i , j k Δ y C i , j + 1 n − C i , j − 1 n 2 Δ y = D 1 C i + 1 , j n − 2 C i , j n + C i − 1 , j n ( Δ x ) 2 + D 2 C i , j + 1 n − 2 C i , j n + C i , j − 1 n ( Δ y ) 2 ,

where u i , j = 0.01 + 0.005 x i − 0.005 y j , v i , j = − 0.01 − 0.005 x i + 0.005 y j , D 1 = 0.0004 , and D 2 = 0.0004 .

Eq. (23) can be written as:

C i , j n + 1 = C i , j n − 0.01 k C i , j n − u i , j k u i , j k Δ x C i , j n − C i − 1 , j n Δ x + 1 − u i , j k Δ x C i + 1 , j n − C i − 1 , j n 2 Δ x − v i , j k v i , j k Δ y C i , j n − C i , j − 1 n Δ y + 1 − v i , j k Δ y C i , j + 1 n − C i , j − 1 n 2 Δ y + D 1 k ( Δ x ) 2 ( C i + 1 , j n − 2 C i , j n + C i − 1 , j n ) + D 2 k ( Δ y ) 2 ( C i , j + 1 n − 2 C i , j n + C i , j − 1 n ) .

The amplification factor is given by:

(24) ξ = 1 − 0.01 k − u i , j k u i , j k Δ x 1 − e − I ω x Δ x + 1 − u i , j k Δ x 2 I sin ( ω x ) 2 Δ x − v i , j k v i , j k Δ y 1 − e − I ω y Δ y + 1 − v i , j k Δ y 2 I sin ω y 2 Δ y + D 1 k ( Δ x ) 2 ( 2 cos ( ω x ) − 2 ) + D 2 k ( Δ y ) 2 ( 2 cos ( ω y ) − 2 ) .

To find range of values of k when Δ x = Δ y = 0.05 , we use the approach of Hindmarsh et al. [25].

Case 1

We fix ω x = π and ω y = π . We substitute u i , j and v i , j in terms of x i and y j and replace Δ x and Δ y by 0.05 in Eq. (24) to obtain

ξ = 1 − 1.29 k − 1,600 k 2 ( 0.01 + 0.005 x i − 0.005 y j ) 2 .

We obtain a 3D plot of ∣ ξ ∣ vs x i ∈ [ 0 , 1 ] vs y j ∈ [ 0 , 1 ] and increase k gradually. Range of values of k for stability is 0 < k < 1.16 .

Case 2

Here, we use the approximations sin ( ω x ) ≈ ω x and cos ( ω x ) ≈ 1 − ω x 2 2 . Then, the amplification factor is given by:

ξ ≃ 1 − 0.01 k − 0.002 k 2 ( ω x 2 + ω y 2 ) + 0.02 k 2 ω x 2 ( − x i + y j ) + 0.02 k 2 ω y 2 ( − x i + y j ) + 0.1 k ω y ( x i − y j ) − 0.005 k 2 ω x 2 ( x i 2 + y j 2 ) + 0.01 x i y j k 2 ( ω x 2 + ω y 2 ) + 0.2 I k ( − ω x + ω y ) + 0.1 k ω x ( − x i + y j ) − 0.005 k 2 ω y 2 ( x i 2 + y j 2 ) − 0.16 k ( ω x 2 + ω y 2 ) .

This gives

∣ ξ ∣ 2 ≃ 1 + ω x 2 [ 0.0002 k 3 ( x i − y j ) − 0.32 k + 0.0016 k 2 − 0.0001 x i y j k 3 + 0.0002 k 3 + 0.00005 k 3 ( x i 2 + y j 2 ) ] + ω y 2 [ 0.0002 k 3 ( − x i + y j ) − 0.32 k + 0.0016 k 2 − 0.0001 x i y j k 3 + 0.0002 k 3 + 0.00005 k 3 ( x i 2 + y j 2 ) ] − 0.01 k + 0.000025 k 2 ≃ 1 + 0.0001 k 2 − 0.002 k .

Solving for ∣ ξ ∣ ≤ 1 with k > 0 , we obtain 0 < k ≤ 200 .

Hence, the intersection of the two inequalities gives the region of the stability as 0 < k ≤ 1.16 .

4.2 Du-Fort–Frankel method

The scheme is given by:

C i , j n + 1 − C i , j n − 1 2 k + u i + 1 , j − u i − 1 , j 2 Δ x C i , j n + u i , j C i + 1 , j n − C i − 1 , j n 2 Δ x + v i , j + 1 − v i , j − 1 2 Δ y C i , j n + v i , j C i , j + 1 n − C i , j − 1 n 2 Δ y = D 1 C i + 1 , j n − C i , j n + 1 − C i , j n − 1 + C i − 1 , j n ( Δ x ) 2 + D 2 C i , j + 1 n − C i , j n + 1 − C i , j n − 1 + C i , j − 1 n ( Δ y ) 2 ,

where u i , j = 0.01 + 0.005 x i − 0.005 y j , v i , j = − 0.01 − 0.005 x i + 0.005 y j , D 1 = 0.0004 , and D 2 = 0.0004 .

We fix Δ x = Δ y = 0.05 . The amplification factor satisfies the following equation:

(25) ξ − ξ − 1 2 k + ( 0.01 + 0.005 x i + 1 − 0.005 y j ) − ( 0.01 + 0.005 x i − 1 − 0.005 y j ) 2 ( 0.05 ) + 0.01 + 0.005 x i − 0.005 y j 2 ( 0.05 ) ( 2 I sin ( ω x ) ) + ( − 0.01 − 0.005 x i + 0.005 y j + 1 ) − ( − 0.01 − 0.005 x i + 0.005 y j − 1 ) 2 ( 0.05 ) + − 0.01 − 0.005 x i + 0.005 y j 2 ( 0.05 ) ( 2 I sin ( ω y ) ) = 0.0004 ( 0.05 ) 2 ( e I ω x − ξ − ξ − 1 + e − I ω x ) + 0.0004 ( 0.05 ) 2 ( e I ω y − ξ − ξ − 1 + e − I ω y ) .

We let x i + 1 = x 3 , x i = x 2 and x i − 1 = x 1 . and on simplification, we obtain

ξ − ξ − 1 + 0.10 k x 3 − 0.10 k x 1 − 0.4 I k sin ( ω y ) − 0.20 I k sin ( ω x ) y 2 + 0.20 I k sin ( ω x ) x 2 + 0.10 k y 3 − 0.10 k y 1 + 0.4 I k sin ( ω x ) + 0.20 I k sin ( ω y ) y 2 − 0.20 I k sin ( ω y ) x 2 − 0.64 k cos ( ω x ) + 0.64 k ξ + 0.64 k ξ − 1 − 0.64 k cos ( ω y ) = 0 .

Case 1:

We fix ω x = π and ω y = π .

The amplification factor satisfies the following equation:

( 1 + 0.64 k ) ξ 2 − 1 + 0.10 k ξ ( x 3 − x 1 ) + 0.10 k ξ ( y 3 − y 1 ) + 1.28 k ξ + 0.64 k = 0 .

For suitable values of x 1 , y 1 chosen, the difference x 3 − x 1 and y 3 − y 1 will be constant. We have x 3 − x 1 = 2 Δ x = 0.1 and y 3 − y 1 = 2 Δ y = 0.1 . We thus have

( 1 + 0.64 k ) ξ 2 − 1 + 1.3 k ξ + 0.64 k = 0 .

Solving for ∣ ξ ∣ ≤ 1 , we obtain

k ∈ [ 0 , ∞ ) .

Case 2:

We consider the case ω x → 0 and ω y → 0 , and Eq. (25) gives

( 1 + 0.64 k ) ξ 2 − 1 + 0.10 k ξ ( x 3 − x 1 ) − 0.4 I k ω y ξ + 0.20 I k ω y y 2 ξ + 0.20 I k ω x x 2 + 0.10 k ξ ( y 3 − y 1 ) ξ + 0.4 I k ω x ξ − 0.20 I k ω x y 2 ξ − 0.20 I k ω y x 2 ξ − 1.28 k ξ + 0.32 k ω x 2 ξ + 0.64 k + 0.32 k ω y 2 ξ = 0 .

As ω x , ω y ≃ 0 , we therefore obtain

( 1 + 0.64 k ) ξ 2 − 1 + 0.10 k ξ ( 0.1 ) + 0.10 k ξ ( 0.1 ) − 1.28 k ξ + 0.64 k = 0 .

Solving ∣ ξ ∣ ≤ 1 , we obtain k ∈ [ 0 , 8.873565 ] .

Hence, the range of values of k for stability is 0 ≤ k ≤ 8.873565 , when Δ x = Δ y = 0.05 .

4.3 NSFD

When NSFD is used to solve Problem 2, we obtain the following scheme:

C i , j n + 1 − C i , j n ϕ ( k ) + 0.005 C i , j n + u i , j C i , j n − C i − 1 , j n ψ ( Δ x ) + 0.005 C i , j n + v i , j C i , j n − C i , j − 1 n ψ ( Δ y ) = D 1 C i + 1 , j n − 2 C i , j n + C i − 1 , j n [ ψ ( Δ x ) ] 2 + D 2 C i , j + 1 n − 2 C i , j n + C i , j − 1 n [ ψ ( Δ y ) ] 2 ,

where ϕ ( k ) = e k − 1 , ψ ( Δ x ) = e Δ x − 1 , and ψ ( Δ y ) = e Δ y − 1 . A single expression for the scheme is

(26) C i , j n + 1 = 1 − 0.01 ϕ ( k ) − u i , j ϕ ( k ) ψ ( Δ x ) − v i , j ϕ ( k ) ψ ( Δ y ) − 2 D 1 ϕ ( k ) [ ψ ( Δ x ) ] 2 − 2 D 2 ϕ ( k ) [ ψ ( Δ y ) ] 2 C i , j n + D 1 ϕ ( k ) [ ψ ( Δ x ) ] 2 C i + 1 , j n + D 2 ϕ ( k ) [ ψ ( Δ y ) ] 2 C i , j + 1 n + u i ϕ ( k ) ψ ( Δ x ) + D 1 ϕ ( k ) [ ψ ( Δ x ) ] 2 C i − 1 , j n + v j ϕ ( k ) ψ ( Δ y ) + D 2 ϕ ( k ) [ ψ ( Δ y ) ] 2 C i , j − 1 n .

We choose ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 1 4 . Eq. (26) reduces to

C i , j n + 1 = 1 − 0.01 ϕ ( k ) − u i , j ϕ ( k ) ψ ( Δ x ) − v i , j ϕ ( k ) ψ ( Δ y ) − D 1 2 − D 2 2 C i , j n + D 1 4 C i + 1 , j n + D 2 4 C i , j + 1 n + u i , j ϕ ( k ) ψ ( Δ x ) + D 1 4 C i − 1 , j n + v i , j ϕ ( k ) ψ ( Δ y ) + D 2 4 C i , j − 1 n ,

where D 1 = D 2 = 0.0004 . We fix Δ x = Δ y = 0.05 , and this gives k = ln ( 1 + 0.25 × 0.0 5 2 ) ≃ 6.25 × 1 0 − 4 .

Using these values, the coefficient of C i , j n is 0.9992, coefficient of C i − 1 , j n is

3.2190 × 1 0 − 4 + 6.0951 × 1 0 − 5 x i − 6.0951 × 1 0 − 5 y j ,

and the coefficient of C i , j − 1 n is 7.8099 × 1 0 − 5 − 6.0951 × 1 0 − 5 x i + 6.0951 × 1 0 − 5 y j . We obtain the plots of coefficients of C i − 1 , j n and C i , j − 1 n versus x ∈ [ 0 , 1 ] versus y ∈ [ 0 , 1 ] in Figure 8, and we observe that these coefficients are all positive. This confirms that the scheme replicates the positivity of the continuous model when

$Figure 8 Plots of coefficients of C i ‒ 1 , j n {C}_{i‒1,j}^{n} and C i , j ‒ 1 n {C}_{i,j‒1}^{n} vs x ∈ [ 0 , 1 ] x\in \left[0,1] vs y ∈ [ 0 , 1 ] y\in \left[0,1] : (a) coefficient of C i ‒ 1 , j n {C}_{i‒1,j}^{n} and (b) coefficient of C i , j ‒ 1 n {C}_{i,j‒1}^{n} .$

Figure 8

Plots of coefficients of C i ‒ 1 , j n and C i , j ‒ 1 n vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] : (a) coefficient of C i ‒ 1 , j n and (b) coefficient of C i , j ‒ 1 n .

ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 1 4 .

5 Numerical results for Problem 1

We choose to display results at times, say T = 1 and T = 5 . We fix Δ x = Δ y = 0.05 and run experiment with k = 0.01 and k = 0.001 for both LW and Du-Fort–Frankel methods, whereas the time step for the NSFD is fixed at k = ln ( 1 + 0.5 ( 0.0 5 2 ) ) ≃ 1.249 × 1 0 − 3 . We obtain reasonable results as shown in Figures 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19. From those results obtained, it is observed that the LW method gives best approximation followed by NSFD and Du-Fort–Frankel, respectively.

$Figure 9 Plot of numerical solution vs x x vs y y using LW at time T = 1 T=1 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k : (a) k = 0.01 k=0.01 and (b) k = 0.001 k=0.001 .$

Figure 9

Plot of numerical solution vs x vs y using LW at time T = 1 using Δ x = Δ y = 0.05 and some values of k : (a) k = 0.01 and (b) k = 0.001 .

$Figure 10 Plot of numerical solution vs x x vs y y using LW at time T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 some values of k k : (a) k = 0.01 k=0.01 and (b) k = 0.001 k=0.001 .$

Figure 10

Plot of numerical solution vs x vs y using LW at time T = 5 using Δ x = Δ y = 0.05 some values of k : (a) k = 0.01 and (b) k = 0.001 .

$Figure 11 Plot of numerical solution vs x x vs y y using Du-Fort–Frankel at time T = 1 T=1 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k : (a) k = 0.01 k=0.01 and (b) k = 0.001 k=0.001 .$

Figure 11

Plot of numerical solution vs x vs y using Du-Fort–Frankel at time T = 1 using Δ x = Δ y = 0.05 and some values of k : (a) k = 0.01 and (b) k = 0.001 .

$Figure 12 Plot of numerical solution vs x x vs y y using Du-Fort–Frankel at time T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 some values of k k : (a) k = 0.01 k=0.01 and (b) k = 0.001 k=0.001 .$

Figure 12

Plot of numerical solution vs x vs y using Du-Fort–Frankel at time T = 5 using Δ x = Δ y = 0.05 some values of k : (a) k = 0.01 and (b) k = 0.001 .

$Figure 13 Plot of numerical solution vs x x vs y y using NSFD at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 1.249 × 1 0 ‒ 3 k=1.249\times 1{0}^{‒3} : (a) T = 1 T=1 and (b) T = 5 T=5 .$

Figure 13

Plot of numerical solution vs x vs y using NSFD at times T = 1 and T = 5 using Δ x = Δ y = 0.05 and k = 1.249 × 1 0 ‒ 3 : (a) T = 1 and (b) T = 5 .

$Figure 14 Absolute errors vs x x vs y y , using the LW at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.01 k=0.01 : (a) T = 1 T=1 and (b) T = 5 T=5 .$

Figure 14

Absolute errors vs x vs y , using the LW at times T = 1 and T = 5 using Δ x = Δ y = 0.05 and k = 0.01 : (a) T = 1 and (b) T = 5 .

$Figure 15 Absolute errors vs x x vs y y , using Du-Fort–Frankel at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.01 k=0.01 : (a) T = 1 T=1 with k = 0.01 k=0.01 and (b) T = 5 T=5 with k = 0.01 k=0.01 .$

Figure 15

Absolute errors vs x vs y , using Du-Fort–Frankel at times T = 1 and T = 5 using Δ x = Δ y = 0.05 and k = 0.01 : (a) T = 1 with k = 0.01 and (b) T = 5 with k = 0.01 .

$Figure 16 Absolute errors vs x x vs y y , using NSFD at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 1.249 × 1 0 ‒ 3 k=1.249\times 1{0}^{‒3} : (a) T = 1 T=1 and (b) T = 5 T=5 .$

Figure 16

Absolute errors vs x vs y , using NSFD at times T = 1 and T = 5 using Δ x = Δ y = 0.05 and k = 1.249 × 1 0 ‒ 3 : (a) T = 1 and (b) T = 5 .

$Figure 17 Relative errors vs x x vs y y , using LW at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.01 k=0.01 : (a) T = 1 T=1 and (b) T = 5 T=5 .$

Figure 17

Relative errors vs x vs y , using LW at times T = 1 and T = 5 using Δ x = Δ y = 0.05 and k = 0.01 : (a) T = 1 and (b) T = 5 .

$Figure 18 Relative errors vs x x vs y y , using Du-Fort–Frankel at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 0.01 k=0.01 : (a) T = 1 T=1 and (b) T = 5 T=5 .$

Figure 18

Relative errors vs x vs y , using Du-Fort–Frankel at times T = 1 and T = 5 using Δ x = Δ y = 0.05 and k = 0.01 : (a) T = 1 and (b) T = 5 .

$Figure 19 Relative errors vs x x vs y y , using NSFD at times T = 1 T=1 and T = 5 T=5 using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 with k = 1.249 × 1 0 − 3 k=1.249\times 1{0}^{-3} : (a) T = 1 T=1 and (b) T = 5 T=5 .$

Figure 19

Relative errors vs x vs y , using NSFD at times T = 1 and T = 5 using Δ x = Δ y = 0.05 with k = 1.249 × 1 0 − 3 : (a) T = 1 and (b) T = 5 .

Plots of numerical profiles vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] at T = 1 and T = 5 are displayed in Figures 9–13. Plots of absolute errors and relative errors vs x vs y are depicted in Figures 14 and 19. We also show numerically based on Figure 20 that LW is unstable at T = 40 when k = 0.15625 , and Δ x = Δ y = 0.05 , and this provides some validation to the stability region obtained on page 5, subsection 3.1 ( k ≤ 0.14 ).

$Figure 20 Plots of exact and numerical solutions vs x x vs y y using LW scheme with k = 0.15625 , Δ x = Δ y = 0.05 k=0.15625,\Delta x=\Delta y=0.05 at time T = 40 T=40 : (a) Exact solution and (b) numerical solution.$

Figure 20

Plots of exact and numerical solutions vs x vs y using LW scheme with k = 0.15625 , Δ x = Δ y = 0.05 at time T = 40 : (a) Exact solution and (b) numerical solution.

Information about the relative errors is obtained from Figures 17–19. At T = 1 , all the three schemes are efficient. Maximum relative errors using LW, NSFD, and Du-Fort–Frankel are 0.025, 0.3, and 2.5%, respectively. At T = 5 , only LW and NSFD have maximum relative error less than 5%. Maximum relative errors using LW, NSFD and Du-Fort–Frankel are 0.008, 1.5, and 10, respectively, at T = 5 .

To compute numerically the rate of convergence, we choose time T = 0.1 . For the LW scheme, we chose to work with k h 2 as constant. Thus, starting with k = 0.05 and Δ x = Δ y = 0.05 , k is divided by 2 2 whenever the spatial step size is divided by 2. From Table 2, we deduce that the numerical rate of convergence in space for the LW and Du-Fort–Frankel methods is two, while that for NSFD is one.

Table 2

Rate of convergence in space when the three methods are used to solve Problem 1 at time 0.1

Schemes	h	l 2 error	Rate ( l 2 )	CPU time
LW	0.0500	5.1747 × 1 0 ‒ 6	—	0.053
	0.0250	1.3084 × 1 0 ‒ 6	1.9837	0.071
	0.0125	3.2793 × 1 0 ‒ 7	1.9963	0.175
	0.00625	8.2028 × 1 0 ‒ 8	1.9992	1.534
Du-Fort–Frankel	0.0500	2.3526 × 1 0 ‒ 3	—	0.073
	0.0250	1.4019 × 1 0 ‒ 3	0.7468	0.076
	0.0125	5.3364 × 1 0 ‒ 4	1.3934	0.086
	0.00625	1.5811 × 1 0 ‒ 4	1.7549	0.122
NSFD	0.0500	4.0131 × 1 0 ‒ 4	—	0.070
	0.0250	2.0688 × 1 0 ‒ 4	0.9559	0.368
	0.0125	1.0469 × 1 0 ‒ 4	0.9827	12.371

6 Numerical results for Problem 2

We present results of Problem 2 at times T = 1 and T = 5 in Figures 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31.

$Figure 21 Surface plots of numerical solution using LW vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 1 T=1 : (a) initial, (b) LW using k = 0.01 k=0.01 , and (c) LW using k = 0.1 k=0.1 .$

Figure 21

Surface plots of numerical solution using LW vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 1 : (a) initial, (b) LW using k = 0.01 , and (c) LW using k = 0.1 .

$Figure 22 Contour plots of numerical solution using LW vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 1 T=1 : (a) initial, (b) LW using k = 0.01 k=0.01 , and (c) LW using k = 0.1 k=0.1 .$

Figure 22

Contour plots of numerical solution using LW vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 1 : (a) initial, (b) LW using k = 0.01 , and (c) LW using k = 0.1 .

$Figure 23 Surface plots of numerical solution using LW vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 5 T=5 : (a) LW using k = 0.01 k=0.01 , (b) LW using k = 0.1 k=0.1 , and (c) LW using k = 1 k=1 .$

Figure 23

Surface plots of numerical solution using LW vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 5 : (a) LW using k = 0.01 , (b) LW using k = 0.1 , and (c) LW using k = 1 .

$Figure 24 Contour plots of numerical solution using LW vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 5 T=5 : (a) LW using k = 0.01 k=0.01 , (b) LW using k = 0.1 k=0.1 , and (c) LW using k = 1 k=1 .$

Figure 24

Contour plots of numerical solution using LW vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 5 : (a) LW using k = 0.01 , (b) LW using k = 0.1 , and (c) LW using k = 1 .

$Figure 25 Surface plots of numerical solution using Du-Fort–Frankel vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 1 T=1 : (a) initial, (b) k = 0.01 k=0.01 , and (c) k = 0.1 k=0.1 .$

Figure 25

Surface plots of numerical solution using Du-Fort–Frankel vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 1 : (a) initial, (b) k = 0.01 , and (c) k = 0.1 .

$Figure 26 Surface and contour plots of numerical solution using Du-Fort–Frankel vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 1 T=1 : (a) initial, (b) k = 0.01 k=0.01 , and (c) k = 0.1 k=0.1 .$

Figure 26

Surface and contour plots of numerical solution using Du-Fort–Frankel vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 1 : (a) initial, (b) k = 0.01 , and (c) k = 0.1 .

$Figure 27 Surface plots of numerical solution using Du-Fort–Frankel vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 5 T=5 : (a) k = 0.01 k=0.01 , (b) k = 0.1 k=0.1 , and (c) k = 1 k=1 .$

Figure 27

Surface plots of numerical solution using Du-Fort–Frankel vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 5 : (a) k = 0.01 , (b) k = 0.1 , and (c) k = 1 .

$Figure 28 Contour plots of numerical solution using Du-Fort–Frankel vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and some values of k k at time T = 5 T=5 : (a) k = 0.01 k=0.01 , (b) k = 0.1 k=0.1 , and (c) k = 1 k=1 .$

Figure 28

Contour plots of numerical solution using Du-Fort–Frankel vs x vs y using Δ x = Δ y = 0.05 and some values of k at time T = 5 : (a) k = 0.01 , (b) k = 0.1 , and (c) k = 1 .

$Figure 29 Surface plots of numerical solution using NSFD vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k ≃ 1.249 × 1 0 ‒ 3 k\simeq 1.249\times 1{0}^{‒3} at time T = 1 T=1 and T = 5 T=5 : (a) initial, (b) NSFD at T = 1 T=1 , and (c) NSFD at T = 5 T=5 .$

Figure 29

Surface plots of numerical solution using NSFD vs x vs y using Δ x = Δ y = 0.05 and k ≃ 1.249 × 1 0 ‒ 3 at time T = 1 and T = 5 : (a) initial, (b) NSFD at T = 1 , and (c) NSFD at T = 5 .

$Figure 30 Contour plots of numerical solution using NSFD vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k ≃ 1.249 × 1 0 ‒ 3 k\simeq 1.249\times 1{0}^{‒3} at time T = 1 T=1 and T = 5 T=5 : (a) initial, (b) NSFD at T = 1 T=1 , and (c) NSFD at T = 5 T=5 .$

Figure 30

Contour plots of numerical solution using NSFD vs x vs y using Δ x = Δ y = 0.05 and k ≃ 1.249 × 1 0 ‒ 3 at time T = 1 and T = 5 : (a) initial, (b) NSFD at T = 1 , and (c) NSFD at T = 5 .

$Figure 31 Plots of numerical results vs x ∈ [ 0 , 1 ] x\in \left[0,1] vs y ∈ [ 0 , 1 ] y\in \left[0,1] using LW when Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 at time T = 10 T=10 : (a) LW when k = 0.5 k=0.5 , (b) LW when k = 1 k=1 , and (c) LW when k = 2 k=2 .$

Figure 31

Plots of numerical results vs x ∈ [ 0 , 1 ] vs y ∈ [ 0 , 1 ] using LW when Δ x = Δ y = 0.05 at time T = 10 : (a) LW when k = 0.5 , (b) LW when k = 1 , and (c) LW when k = 2 .

We display results at time T = 1 using the three methods in Figures 21, 22, 25, 26, 29, and 30. The range of the initial concentration is 1–10. At T = 1 , the range of concentration is 1–5.5 from all the three schemes used based on Figures 21, 25, and 29. Some dispersive oscillations are seen when Du-Fort–Frankel is used.

Results at time T = 5 are displayed in Figures 23, 24, 27, 28, 29, and 30. The range of concentration using LW and NSFD is 1 to 1.8, while the range using Du-Fort–Frankel is 1 to 1.9. We observe some unbounded values in the numerical solution of LW depicted by Figure 31(c). This validates the result from stability analysis where we deduced that range of k when Δ x = Δ y = 0.05 for LW is 0 < k ≤ 1.16 .

7 Conclusion

In this article, three finite difference schemes, LW, Du-Fort–Frankel, and NSFD, were constructed and used to solve a 2D advection–diffusion equation with nonconstant coefficients. Two problems were considered whereby one has an exact solution. First, for the problem with exact solution, we deduce that the LW is the most efficient scheme followed by NSFD and least performing scheme is Du-Fort–Frankel. Secondly, it is known that the stability analysis of standard finite difference schemes discretising 1D advection–diffusion equation is not straightforward and here, complication arises as we are considering 2D nonconstant coefficient advection–diffusion equation. However, we have managed to obtain the range of values of k for stability of LW and Du-Fort–Frankel schemes using von Neumann stability analysis coupled with other techniques. Third, in the case of NSFD, since we are dealing with non-negative variables, we obtain the condition for which the scheme replicates positivity of the solution of the continuous model and this is not too time-consuming to achieve. Fourth, the numerical rate of convergence of the three methods is obtained. We compare the CPU times of the three algorithms when used to solve Problem 1. Moreover, we display the results of Problem 2 using the three methods and observe that some dispersive oscillations are present in the solutions from the Du-Fort–Frankel scheme.

Future work will involve the use of LW and NSFD schemes to solve some real-life problems dealing with irregular grids and other coefficients of dissipation.

Acknowledgments

A.R. Appadu is grateful to Nelson Mandela University to be able to carry out this work during his sabbatical leave from January 2023 to June 2023. The authors are grateful to the three anonymous reviewers for the feedback, which enabled them to significantly improve this article in terms of both content and presentation.

Funding information: AR Appadu is grateful to NMU for allowing him to use his publication funds to pay for open access fees.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission.
Conflict of interest: The authors state no conflict of interest.
Data availability statement: Data is generated from initial conditions.

References

[1] Parlange JY. Water transport in soils. Annual Rev Fluid Mechanics. 1980;12(1):77–102. 10.1146/annurev.fl.12.010180.000453Search in Google Scholar

[2] Tian ZF and Yu P. A high-order exponential scheme for solving 1D unsteady convection-diffusion equations. J Comput Appl Math. 2011;235(8):2477–91. 10.1016/j.cam.2010.11.001Search in Google Scholar

[3] Kumar RK, Price JE, Sargent NS, and Tarin D. Binding of 125I-labelled concanavalin A by cells from spontaneously arising murine mammary carcinomas and experimentally induced metastases. Anticancer Res. 1983;3(5):343–6. Search in Google Scholar

[4] Szymkiewicz R. Numerical modeling in open channel hydraulics. Vol. 83. Dordrecht: Springer Science & Business Media, 2010. 10.1007/978-90-481-3674-2Search in Google Scholar

[5] Sommeijer BP and Kok J. Implementation and performance of the time integration of a 3D transport model. Int J Numer Meth Fluids. 1995;20(3):213–31. 10.1002/fld.1650200303Search in Google Scholar

[6] Dehghan M. Numerical solution of the three-dimensional advection-diffusion equation. Appl Math Comput. 2004;150(1):5–19. 10.1016/S0096-3003(03)00193-0Search in Google Scholar

[7] Thongmoon M, McKibbin R. A comparison of some numerical methods for the advection-diffusion equation. Res Lett Inf Math Sci. 2006;10:49–62.Search in Google Scholar

[8] Appadu AR, Djoko JK, Gidey HH. A computational study of three numerical methods for some advection-diffusion problems. Appl Math Comput. 2016;272:629–47. 10.1016/j.amc.2015.03.101Search in Google Scholar

[9] Appadu AR, Djoko JK, Gidey HH. Performance of some finite difference methods for a 3D advection-diffusion equation. Revista de la Real Academia de Ciencias Exactas, Físicas y Naturales. Serie A. Matemáticas. 2018;112(4):1179–210. 10.1007/s13398-017-0414-7Search in Google Scholar

[10] Takacs LL. A two-step scheme for the advection equation with minimized dissipation and dispersion errors. Monthly Weather Rev. 1985;113(6):1050–65. 10.1175/1520-0493(1985)113<1050:ATSSFT>2.0.CO;2Search in Google Scholar

[11] Courant R, Friedrichs K, Lewy H. Über die partiellen Differenzengleichungen der mathematischen Physik. Math Annalen. 1928;100(1):32–74. 10.1007/BF01448839Search in Google Scholar

[12] Von Neumann J, Richtmyer RD. A method for the numerical calculation of hydrodynamic shocks. J Appl Phys. 1950;21(3):232–7. 10.1063/1.1699639Search in Google Scholar

[13] Kreiss H-O. Stability theory for difference approximations of mixed initial boundary value problems. I. Math Comput. 1968;22(104):703–14. 10.1090/S0025-5718-1968-0241010-7Search in Google Scholar

[14] Winnicki I, Jasinski J, Pietrek S. New approach to the Lax–Wendroff modified differential equation for linear and nonlinear advection. Numer Meth Partial Differ Equ. 2019;35(6):2275–304. 10.1002/num.22412Search in Google Scholar

[15] Hirt CW. Heuristic stability theory for finite-difference equations. J Comput Phys. 1968;2(4):339–55. 10.1016/0021-9991(68)90041-7Search in Google Scholar

[16] Ru-xun L, Zhao-hui Z. The remainder-effect analysis of finite difference schemes and the applications. Appl Math Mechanics. 1995;16:87–96. 10.1007/BF02453778Search in Google Scholar

[17] Morton KW, Mayers DF. Numerical solution of partial differential equations. Cambridge: Cambridge University Press; 1994. Search in Google Scholar

[18] Appadu AR. Numerical solution of the 1D advection-diffusion equation using standard and nonstandard finite difference schemes. J Appl Math. 2013;2013:734374. 10.1155/2013/734374Search in Google Scholar

[19] Appadu AR. Optimized composite finite difference schemes for atmospheric flow modeling. Numer Meth Partial Differ Equ. 2019;35(6):2171–92. 10.1002/num.22407Search in Google Scholar

[20] Hutomo GD, Kusuma J, Ribal A., Mahie AG, Aris N. Numerical solution of 2D advection-diffusion equation with variable coefficient using Du-Fort–Frankel method. in: Journal of Physics : Conference Series. Vol. 1180. Bristol, United Kingdom: IOP Publishing; 2019. p. 012009. 10.1088/1742-6596/1180/1/012009Search in Google Scholar

[21] Fromm J. The time dependent flow of an incompressible viscous fluid. Methods Comput Phys. 1964;3:346–82. Search in Google Scholar

[22] Siemieniuch J, Gladwell I. Analysis of explicit difference methods for a diffusion-convection equation. Int J Numer Meth Eng. 1978;12(6):899–916. 10.1002/nme.1620120603Search in Google Scholar

[23] Rigal A. Stability analysis of explicit finite difference schemes for the Navier-Stokes equations. Int J Numer Meth Eng. 1979;14(4):617–20. 10.1002/nme.1620140411Search in Google Scholar

[24] Clancy RM. A note on finite differencing of the advection-diffusion equation. Monthly Weather Review. 1981;109(8):1807–9. 10.1175/1520-0493(1981)109<1807:ANOFDO>2.0.CO;2Search in Google Scholar

[25] Hindmarsh AC, Gresho PM, Griffiths DF. The stability of explicit euler time-integration for certain finite difference approximations of the multi-dimensional advection-diffusion equation. Int J Numer Meth Fluids. 1984;4(9):853–97. 10.1002/fld.1650040905Search in Google Scholar

[26] Lax PD. Difference schemes with high-order of accuracy for solving hyperbolic equations. Commun Pure Appl Math. 1964;17:381. 10.1002/cpa.3160170311Search in Google Scholar

[27] Smith R, Tang Y. Optimal and near-optimal advection-diffusion finite-difference schemes. V. Error propagation. Proc R Soc London Ser A Math Phys Eng Sci. 2001;457(2008):803–16. 10.1098/rspa.2000.0693Search in Google Scholar

[28] Du Fort EC, Frankel SP. Stability conditions in the numerical treatment of parabolic differential equations. Math. Comp. 1953;17:35–152. 10.1090/S0025-5718-1953-0059077-7 Search in Google Scholar

[29] Kojouharov HV, Cheni BM. Nonstandard methods for advection-diffusion-reaction. Applications of Nonstandard Finite Difference Schemes. Singapore: World Scientific; 2000. p. 55. 10.1142/9789812813251_0002Search in Google Scholar

[30] Mickens RE, Washington TM. A note on a positivity preserving nonstandard finite difference scheme for a modified parabolic reaction-advection-diffusion PDE. J Differ Equ Appl. 2020;26(11–12):1423–7. 10.1080/10236198.2020.1841755Search in Google Scholar

Received: 2023-11-22

Revised: 2024-01-10

Accepted: 2024-01-17

Published Online: 2024-03-12

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/phys-2023-0195

Keywords for this article

nonconstant coefficient; advection–diffusion; Lax–Wendroff; Du-Fort–Frankel; nonstandard finite difference; stability; rate of convergence

Creative Commons

BY 4.0