Numerical solution of a nonconstant coefficient advection diffusion equation in an irregular domain and analyses of numerical dispersion and dissipation

Appanah Rao Appadu; Hagos Hailu Gidey

doi:10.1515/phys-2025-0137

Article Open Access

Numerical solution of a nonconstant coefficient advection diffusion equation in an irregular domain and analyses of numerical dispersion and dissipation

Appanah Rao Appadu and Hagos Hailu Gidey

Published/Copyright: April 29, 2025

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 23 Issue 1

Abstract

This work is a major extension of our previous work in which we have solved a 2D nonconstant coefficient advection diffusion equation with nonconstant advection and constant diffusion terms on a square domain using the coefficient of dissipation D 1 = D 2 = 0.0004 using three finite difference methods, namely, Lax–Wendroff, Du Fort–Frankel and nonstandard finite difference methods. In this current work, the first novelty is that we solve a 2D nonconstant advection diffusion equation on an irregular domain with a more complicated initial profile and considered five combinations for values of D 1 and D 2 . Moreover, the second novelty is the study of numerical dispersion and dissipation of Lax Wendroff scheme for the five combinations of D 1 and D 2 . Third, we present some numerical profiles from the three methods for the five scenario at two times: T = 0.1 , 1. The fourth novelty is the plot of the modulus of the exact amplification factor, modulus of amplification factor, and relative phase error vs phase angle along x direction vs phase angle along y direction for the Lax–Wendroff scheme at x = y = 0.5 for the five scenarios.

Keywords: advection diffusion; Lax–Wendroff; Du Fort–Frankel; nonstandard finite difference; nonconstant coefficient; stability; irregular domain

1 Introduction

The advection diffusion equation is one of the most important partial differential equations in science and it represents a superposition of two different transport processes: advection and diffusion [1,2]. This model is used to describe several physical phenomena such as transport of pollutants [1], flow in porous media [3], water transport in soils [4], mass transfer [5], and heat transfer in a draining film [6].

Several numerical methods have been developed to solve advection diffusion equation with constant coefficients in one, two, and three dimensions (refer [7–14] and references therein).

A novel finite difference scheme following the work of Dehghan [7] is used along with the Crank–Nicolson and Implicit Chapeau function to solve 3D advection diffusion equation with given initial and boundary conditions [13,15]. The authors compare the performance of the three methods by comparing L 2 -error, L ∞ -error, and some performance indices. Appadu et al. [10] used three numerical methods to solve two test problems described by advection diffusion equations. The first test problem considered has a steep boundary layers near x = 1 and this is a challenging problem as many methods are affected by non-physical oscillation near steep boundaries.

Nonstandard finite difference methods (NSFDs) are an established class of methods to solve reaction diffusion, advection diffusion partial differential equations in particular. Verma and Kayenat [16] derived an exact finite difference using a solitary wave solution and also proposed a nonstandard finite difference schemes for the generalised Burgers–Huxley equation subject to certain initial and boundary conditions. Some work on comparison of exact finite difference methods and NSFDs for two classes of nonlinear advection diffusion reaction equations and for the nonlinear generalised advection diffusion reaction equations are detailed in the study by Kayenat and Verma [17,18]. Verma and Kayenat constructed NSFD schemes for generalised Burgers–Fisher equation [19]. Nonconstant coefficient partial differential equations have applications in many fields such as: engineering science, quantum mechanics, financial mathematics, isomonodromic deformations, Seiberg–Witten invariants and quantum field theory, shallow water waves, plasmas, solitons, and optics [20]. Using the diffusion damping and the variable coefficient advection methodology, El-Nabulsi [21] concluded that the problem of non-uniform Uranium burnup in nuclear reactor may be reduced.

Hutomo et al. [22] solved an advection diffusion equation using the Du Fort–Frankel method on square and irregular domains. The irregular domain considered is based on the lake of Hasanuddin University. The lake is located at Tamalanrea Campus in Makassar and is of length ± 41 m and width ± 34 m .

Hutomo et al. [22] solved the 2D advection diffusion equation

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

where C is the concentration of certain species, x ∈ R is a space variable, t represents the time with t ≥ 0 , u and v are velocity coefficients in x and y directions, D 1 and D 2 are the diffusion coefficients in x and y directions, respectively. They used D 1 = D 2 = 0.00013 and

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

The vector field of the domain of Hasanuddin University lake is shown in Figure 1. To run the numerical experiment, they used the spatial step sizes, Δ x = 0.025 , Δ y = 0.03125 and temporal step size, k = 0.005 . We note that no stability analysis of the numerical method was performed by Hutomo et al. [22].

Figure 1

The vector field of velocity flow of the domain [22].

This study is organised as follows: in Section 2, we describe the numerical experiment considered. In Sections 3, 4, and 5, we construct the three methods namely Lax–Wendroff, Du Fort–Frankel and NSFD and study the stability using approach of Hindmarsch et al. [23] or obtain condition for which the NSFD method replicates positivity of the continuous model. In Section 6, we study the dissipation and dispersion characteristics of Lax–Wendroff scheme. Numerical results are displayed in Section 7. Section 8 highlights the salient features of the study.

2 Numerical experiment

Based on the problem considered by Hutomo et al. [22], which was described in Section 1, we propose to solve the following problem. We solve

(1) ∂ 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 given diffusion coefficients,

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

for x , y ∈ [0, 1] and t ∈ [ 0 , T ] . The initial conditions are given by

C ( x , y , 0 ) = 4 if 0.65 ≤ x ≤ 0.85 , 0.15 ≤ y ≤ 0.35 , and ( x − 0.75 ) 2 + ( y − 0.25 ) 2 ≤ 0.01 0 if 0 ≤ x < 0.5 and y > 0.25 , 0 if 0.25 ≤ x < 0.5 , and y > x , 0 if 0.75 < x ≤ 1 and y ≥ 1.75 − x , 1 otherwise .

We note that C ( x , y , t ) = 1 along the boundary of the domain. The irregular domain is given in Figure 2 and the initial profile is shown in Figure 3.

Figure 2

Domain considered for our study.

Figure 3

Initial profile.

We consider five scenarios.

Scenario 1: D 1 = D 2 = 0.0004 .

Scenario 2: D 1 = D 2 = 0.04 .

Scenario 3: D 1 = D 2 = 0.4 .

Scenario 4: D 1 = 0.04 and D 2 = 0.4 .

Scenario 5: D 1 = 0.4 and D 2 = 0.04 .

For stability, we use the approach of Hindmarsch et al. [23] to find the range of values of k for two cases:

Case 1 ̲

Phase angle along x -direction, ω x = π and phase angle along y -direction, ω y = π . We also fix Δ x = Δ y = 0.05 .

Case 2 ̲

When ω x and ω y tend to zero, we also choose Δ x = Δ y = 0.05 .

3 Derivation and stability of Lax–Wendroff method

The Lax–Wendroff method when used to discretise Eq. (1) is given by [24]

(2) 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 Δ x 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 .

Eq. (2) can be written as

(3) 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 Δ x 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 ) ,

and the amplification factor is

(4) ξ = 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 Δ x 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 ) .

3.1 Scenario 1

A detailed explanation on the stability of Lax–Wendroff scheme for scenario 1 is provided in [24]. We describe briefly some steps involved for the first scenario.

The amplification factor is obtained and we replace D 1 and D 2 by 0.0004. Then, case 1 is considered whereby we fix ω x = π and ω y = π . We replace u i , j and v i , j by 0.01 + 0.005 x i − 0.005 y j and − 0.01 − 0.005 x i + 0.005 y j , respectively. We also fix Δ x = Δ y = 0.05 in Eq. (4) and obtain the amplification factor as

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

We obtain 3D plots of ∣ ξ ∣ vs x ∈ [0, 1] vs y ∈ [0, 1] in Figures 4 and 5. We increase k gradually, starting from a small k say k = 0.001 until ∣ ξ ∣ ≤ 1 is no longer satisfied. The range of values of k is

(5) 0 < k ≤ 1.16 .

$Figure 4 3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for k = 0.001 k=0.001 , 0.01, 0.1 for scenario 1 and case 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 the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ [0, 1] vs y ∈ [0, 1] for k = 0.001 , 0.01, 0.1 for scenario 1 and case 1. (a) k = 0.001 , (b) k = 0.01 , and (c) k = 0.1 .

$Figure 5 3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for k = 1.0 k=1.0 , 1.16, 1.18 for scenario 1 and case 1. (a) k = 1.0 k=1.0 , (b) k = 1.16 k=1.16 , and (c) k = 1.18 k=1.18 .$

Figure 5

3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ [0, 1] vs y ∈ [0, 1] for k = 1.0 , 1.16, 1.18 for scenario 1 and case 1. (a) k = 1.0 , (b) k = 1.16 , and (c) k = 1.18 .

We then consider case 2 which is the situation when ω x and ω y tend to zero. We consider (4), replace u i , j and v i , j by required expressions and fix Δ x = Δ y = 0.05 and use the approximations: sin ( ω x ) ≈ ω x and cos ( ω x ) ≈ 1 − ω x 2 2 . We also fix D 1 = D 2 = 0.0004 and solving for ∣ ξ ∣ 2 ≤ 1 with k > 0 , gives

(6) 0 < k ≤ 100 .

Combining inequalities (5) and (6) gives the range of values of k for stability as 0 < k ≤ 1.16 .

3.2 Scenario 2

We consider Eq. (4) and substitute D 1 , D 2 , by 0.04 and Δ x , Δ y by 0.05. Case 1 is considered where we fix ω x = ω y = π . We next replace u i , j and v i , j in terms of x i and y j . This gives

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

3D plots of ∣ ξ ∣ vs x ∈ [0, 1] vs y ∈ [0, 1] are obtained in Figures 6 and 7. Range of values of k for stability is 0 < k ≤ 0.015 .

$Figure 6 3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] when k = 0.001 k=0.001 , 0.01 for scenario 2 and case 1. (a) k = 0.001 k=0.001 and (b) k = 0.01 k=0.01 .$

Figure 6

3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ [0, 1] vs y ∈ [0, 1] when k = 0.001 , 0.01 for scenario 2 and case 1. (a) k = 0.001 and (b) k = 0.01 .

$Figure 7 3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] when k = 0.015 k=0.015 , 0.016 for scenario 2 and case 1. (a) k = 0.015 k=0.015 and (b) k = 0.016 k=0.016 .$

Figure 7

3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ [0, 1] vs y ∈ [0, 1] when k = 0.015 , 0.016 for scenario 2 and case 1. (a) k = 0.015 and (b) k = 0.016 .

We next make use of Eq. (4) and study stability of the method for case 2. This gives ∣ ξ ∣ 2 ≈ ( 1 − 0.01 k ) 2 and solving ∣ ξ ∣ 2 ≤ 1 gives 0 < k ≤ 100 .

Hence, the range of values of k for stability for scenario 2 with Δ x = Δ y = 0.05 is 0 < k ≤ 0.015 .

3.3 Scenario 3

We consider Eq. (4) and work with case 1, i.e. ω x = ω y = π . We then substitute D 1 , D 2 by 0.4 and Δ x , Δ y by 0.05 and replace u i , j , v i , j in terms of x i and y j . This gives

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

The range of values of k for stability is 0 < k ≤ 0.0015 as shown in Figure 8.

$Figure 8 3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for some values of k k for scenario 3 and case 1. (a) k = 0.001 k=0.001 , (b) k = 0.0015 k=0.0015 , and (c) k = 0.0016 k=0.0016 .$

Figure 8

3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ [0, 1] vs y ∈ [0, 1] for some values of k for scenario 3 and case 1. (a) k = 0.001 , (b) k = 0.0015 , and (c) k = 0.0016 .

We now consider Eq. (4) for case 2. We obtain ∣ ξ ∣ 2 ≈ ( 1 − 0.01 k ) 2 and solving for ∣ ξ ∣ 2 ≤ 1 gives 0 < k ≤ 100 .

Hence, the range of values of k for stability for scenario 3 is 0 < k ≤ 0.0015 .

3.4 Scenario 4

Starting with Eq. (4), we substitute D 1 , D 2 by 0.04 and 0.4, respectively, and fix Δ x = Δ y = 0.05 . Case 1 is considered whereby we fix ω x = ω y = π . Replacing u i , j , v i , j in terms of x i and y j gives the amplification factor

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

3D plots of ∣ ξ ∣ vs x ∈ [0, 1] vs y ∈ [ 0,1 ] are obtained and range of values of k for stability is 0 < k ≤ 0.00284 , as shown in Figure 9.

$Figure 9 3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for some values of k k for scenario 4 and case 1. (a) k = 0.001 k=0.001 , (b) k = 0.0028 k=0.0028 , (c) k = 0.0029 k=0.0029 .$

Figure 9

3D plots of the modulus of the amplification factor of Lax–Wendroff scheme vs x ∈ [0, 1] vs y ∈ [0, 1] for some values of k for scenario 4 and case 1. (a) k = 0.001 , (b) k = 0.0028 , (c) k = 0.0029 .

We then consider case 2. We obtain ∣ ξ ∣ 2 ≈ ( 1 − 0.01 k ) 2 and solving for ∣ ξ ∣ 2 ≤ 1 gives 0 < k ≤ 200 .

Hence, the range of values of k for stability for scenario 4 is 0 < k ≤ 0.00284 .

3.5 Scenario 5

Starting with Eq. (4), we substitute D 1 , D 2 by 0.4 and 0.04, respectively, and fix Δ x = Δ y = 0.05 . We obtain exactly the same amplification factor as in scenario 4. Hence, the range of values of k for stability for scenario 5 with Δ x = Δ y = 0.05 is 0 < k ≤ 0.00284 .

4 Derivation and stability of Du Fort–Frankel method

The Du Fort–Frankel method is a modification of the centred time centred space scheme [22,25] and when used to discretise Eq. (1) is given by [22]

(7) 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 .

The amplification factor satisfies the following equation:

(8) ξ − ξ − 1 2 k + 0.005 ( x i + 1 − x i − 1 ) 2 Δ x + 0.005 ( y j + 1 − y j − 1 ) 2 Δ y + 0.01 + 0.005 x i − 0.005 y j 2 Δ x ( 2 I sin ( ω x ) ) + − 0.01 − 0.005 x i + 0.005 y j 2 Δ y ( 2 I sin ( ω y ) ) = D 1 ( Δ x ) 2 ( e I ω x − ξ − ξ − 1 + e − I ω x ) + D 2 ( Δ y ) 2 ( e I ω y − ξ − ξ − 1 + e − I ω y ) .

Since Δ x = Δ y = 0.05 , Eq. (8) can be rewritten as

(9) ξ − ξ − 1 2 k + 0.01 + 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 2 ( 0.05 ) ( 2 I sin ( ω y ) ) = D 1 ( 0.05 ) 2 ( e I ω x − ξ − ξ − 1 + e − I ω x ) + D 2 ( 0.05 ) 2 ( e I ω y − ξ − ξ − 1 + e − I ω y ) .

We next use the approach of Hindmarsh et al. [23] to find the range of values of k when Δ x = Δ y = 0.05 for the five scenarios described in Section 2.

4.1 Scenario 1

We make use of Eq. (9) and replace D 1 , D 2 by 0.0004.

For case 1, fixing ω x = π and ω y = π , we obtain

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

Solving for ∣ ξ ∣ ≤ 1 , we obtain

(10) 0 < k < ∞ .

For case 2, we obtain the following equation:

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

For ω x , ω y → 0 and on solving ∣ ξ ∣ ≤ 1 , we obtain

(11) 0 < k ≤ 8.873565 .

Combining the two inequalities (10) and (11) gives the range of values of k for stability as 0 ≤ k ≤ 8.873565 .

4.2 Scenario 2

We use (9) and consider case 1. We also replace D 1 and D 2 by 0.04.

Fixing ω x = π and ω y = π gives

( 1 + 64 k ) ξ 2 + 128.02 k ξ + ( 64 k − 1 ) = 0 .

Solving for ∣ ξ ∣ ≤ 1 gives 0 < k ≤ 100 .

For case 2, we obtain

( 1 + 64 k ) ξ 2 − 127.98 k ξ + 0.4 I k ξ ( ω x − ω y ) − 0.2 I k ξ y ( ω x − ω y ) + 0.2 I k ξ x ( ω x − ω y ) + 32 k ξ ( ω x 2 + ω y 2 ) + 64 k − 1 = 0 .

Solving for ∣ ξ ∣ 2 ≤ 1 when ω x → 0 and ω y → 0 gives 0 ≤ k ≤ 0.883918 . Combining the range of values of k for cases 1 and 2 gives 0 < k ≤ 0.883918 for stability.

4.3 Scenario 3

For case 1, replacing D 1 , D 2 by 0.4 in Eq. (9) and fixing ω x = π , ω y = π , we obtain the following equation:

( 1 + 640 k ) ξ 2 + 1280.02 k ξ + 640 k − 1 = 0 .

Solving ∣ ξ ∣ ≤ 1 gives 0 < k ≤ 1 .

For case 2, replacing D 1 , D 2 by 0.4, we obtain

( 1 + 640 k ) ξ 2 − 1279.98 k ξ + 0.4 I k ξ ( ω x − ω y ) − 0.2 I k ξ y ( ω x − ω y ) + 0.2 I k ξ x ( ω x − ω y ) + 320 k ξ ( ω x 2 + ω y 2 ) + 640 k − 1 = 0 .

Solving for ∣ ξ ∣ ≤ 1 when ω x , ω y → 0 gives 0 ≤ k ≤ 0.279510 . Hence, the range of values of k for stability for scenario 3 is 0 < k ≤ 0.279510 .

4.4 Scenario 4

We consider Eq. (9) and replace D 1 , D 2 by 0.04 and 0.4, respectively. Fixing ω x = π and ω y = π gives the quadratic equation

( 1 + 352 k ) ξ 2 + 704.02 k ξ + 352 k − 1 = 0 .

Solving ∣ ξ ∣ 2 ≤ 1 gives 0 < k < 1 .

For case 2, we obtain

( 1 + 352 k ) ξ 2 − 703.98 k ξ + 352 k − 1 = 0 ,

and solving for ∣ ξ ∣ 2 ≤ 1 gives 0 ≤ k ≤ 0.376892 . Hence, the range of values of k for stability for scenario 4 is 0 < k ≤ 0.376891 .

4.5 Scenario 5

We consider Eq. (9) and replace D 1 , D 2 by 0.4 and 0.04, respectively. Fixing ω x = π and ω y = π gives the quadratic equation

( 1 + 352 k ) ξ 2 + 754.02 k ξ + 352 k − 1 = 0 .

Solving for ∣ ξ ∣ 2 ≤ 1 gives 0 < k < 1 .

For case 2, we obtain

( 1 + 352 k ) ξ 2 − 703.98 k ξ + 352 k − 1 = 0 ,

and solving for ∣ ξ ∣ 2 ≤ 1 gives 0 ≤ k ≤ 0.376892 . Hence, the range of values of k for stability for scenario 5 is 0 < k ≤ 0.376892 .

5 Derivation and stability of NSFD

To construct an NSFD for Eq. (1), we use the following approximations [24]:

∂ 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 ϕ ( k ) = e k − 1 , ψ ( Δ x ) = e Δ x − 1 , and ψ ( Δ y ) = e Δ y − 1 . When NSFD is used to discretise Eq. (1), we obtain the following scheme [24]:

(12) C i , j n + 1 − C i , j n ϕ ( k ) + ∂ u ∂ x i C i , j n + u i , j C i , j n − C i − 1 , j n ψ ( Δ x ) + ∂ v ∂ y j 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 .

A single expression for Eq. (12) is given by

C i , j n + 1 = u i , j ϕ ( k ) ψ ( Δ x ) + D 1 ϕ ( k ) ( Δ x ) 2 C i − 1 , j n + v i , j ϕ ( k ) ψ ( Δ y ) + D 2 ϕ ( k ) ( Δ y ) 2 C i , j − 1 n + 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 .

We choose the functional relation

ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 0.5 and obtain

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

We choose Δ x = Δ y = 0.05 . Since ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 0.5 , we obtain k ≈ 1.31350 × 1 0 − 3 .

The coefficients of C i , j n for scenarios 1–5 are 0.999187, 0.919987, 0.199987, 0.559987, and 0.559987, respectively. We obtain plots the coefficients of C i − 1 , j n and C i , j − 1 n vs x ∈ [0, 1] vs y ∈ [0, 1] for the five scenarios in Figures 10 and 11 in order to check if NSFD scheme preserves positivity of the continuous model. Here, we mean that the numerical solutions remain non-negative at any time given non-negative values.

$Figure 10 3D plots of the 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 ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for scenario 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 10

3D plots of the coefficients of C i − 1 , j n and C i , j − 1 n vs x ∈ [0, 1] vs y ∈ [0, 1] for scenario 1. (a) Coefficient of C i − 1 , j n and (b) coefficient of C i , j − 1 n .

$Figure 11 3D plots of the 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 ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for scenario 2. (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 11

3D plots of the coefficients of C i − 1 , j n and C i , j − 1 n vs x ∈ [0, 1] vs y ∈ [0, 1] for scenario 2. (a) Coefficient of C i − 1 , j n and (b) coefficient of C i , j − 1 n .

NSFD preserves positivity of the continuous model for scenarios 2–5 when ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 0.5 and Δ x = Δ y = 0.05 , as shown in Figures 10, 11, 12, 13, to 14. However, it does not satisfy positivity of the continuous model for some values of x and y for scenario 1 as depicted in Figure 10(b). Hence, the NSFD is not useful for scenario 1 as it is not positivity preserving in that situation (Table 1).

$Figure 12 3D plots of the 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 ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for scenario 3. (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 12

3D plots of the coefficients of C i − 1 , j n and C i , j − 1 n vs x ∈ [0, 1] vs y ∈ [0, 1] for scenario 3. (a) Coefficient of C i − 1 , j n and (b) coefficient of C i , j − 1 n .

$Figure 13 3D plots of the 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 ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for scenario 4. (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 13

3D plots of the coefficients of C i − 1 , j n and C i , j − 1 n vs x ∈ [0, 1] vs y ∈ [0, 1] for scenario 4. (a) Coefficient of C i − 1 , j n and (b) coefficient of C i , j − 1 n .

$Figure 14 3D plots of the 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 ∈ x\hspace{0.33em}\in [0, 1] vs y ∈ y\hspace{0.33em}\in [0, 1] for scenario 5. (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 14

3D plots of the coefficients of C i − 1 , j n and C i , j − 1 n vs x ∈ [0, 1] vs y ∈ [0, 1] for scenario 5. (a) Coefficient of C i − 1 , j n and (b) coefficient of C i , j − 1 n .

Table 1

Range of values of k for stability when Δ x = Δ y = 0.05 for the three numerical methods for the five scenarios

Scenario	Lax–Wendroff	Du Fort–Frankel	NSFD
1	0 < k ≤ 1.16	0 < k ≤ 8.873565	Not positive definite
2	0 < k ≤ 0.015	0 < k ≤ 0.883918	k = 1.3135 × 1 0 − 3
3	0 < k ≤ 0.0015	0 < k ≤ 0.279509	k = 1.3135 × 1 0 − 3
4	0 < k ≤ 0.00284	0 < k ≤ 0.376892	k = 1.3135 × 1 0 − 3
5	0 < k ≤ 0.00284	0 < k ≤ 0.376892	k = 1.3135 × 1 0 − 3

6 Numerical dispersion and dissipation of Lax–Wendroff

In this section, we obtain plots of the following quantities vs ω x vs ω y for the five scenarios when the Lax–Wendroff scheme is used to solve Eq. (1)

modulus of the amplification factor.
relative phase error (RPE).

We note that phase angles along the x direction and y direction are denoted by ω x and ω y , respectively. We also plot the modulus of the exact amplification factor vs ω x vs ω y .

6.1 Exact amplification factor

We consider Eq. (1) with u ( x , y ) = 0.01 + 0.005 x − 0.005 y and v ( x , y ) = − 0.01 − 0.005 x + 0.005 y . We use the perturbation for C ( x , y , t ) as e α t e I θ x x e I θ y y [13,26], where α is the dispersion relation. Using this perturbation for C ( x , y , t ) in Eq. (1) gives

α + u ( x , y ) I θ x + 0.01 + v ( x , y ) I θ y = − D 1 θ x 2 − D 2 θ y 2 .

Hence,

α = − I θ x u ( x , y ) − 0.01 − I θ y v ( x , y ) − D 1 θ x 2 − D 2 θ y 2 .

We now obtain the exact amplification factor denoted as ξ exact which is the perturbation for C ( x , y , t + k ) divided by perturbation for C ( x , y , t ) [26,27].

ξ exact = e α k = e ( − I θ x u ( x , y ) ) k e ( − I θ y v ( x , y ) ) k e ( − D 1 θ x 2 − D 2 θ y 2 − 0.01 ) k .

The modulus of the exact amplification factor is given by

∣ ξ exact ∣ = e ( − D 1 θ x 2 − D 2 θ y 2 − 0.01 ) k ,

where θ x = ω x Δ x and θ y = ω y Δ y . Since in this work, we choose Δ x = Δ y = 0.05 , we therefore have

∣ ξ exact ∣ = e ( − 400 D 1 ω x 2 − 400 D 2 ω y 2 − 0.01 ) k .

6.2 RPE

The RPE is a measure of the dispersive characteristics of a scheme [28]. The RPE is calculated as [28]

RPE = arg ( ξ ) arg ( ξ exact ) ,

where ξ is the amplification factor of the numerical method. We can rewrite ξ exact as

ξ exact = e ( − D 1 θ x 2 − D 2 θ y 2 − 0.01 ) k ( cos A + I sin A ) ( cos B + I sin B ) ,

where

A = − θ x u ( x , y ) k , B = − θ y v ( x , y ) k .

The argument of ξ exact is given by

arg ( ξ exact ) = arctan cos A sin B + sin A cos B cos A cos B − sin A sin B .

The resulting expressions for ξ exact and RPE consist of the parameters ω x , ω y , x , and y . We can fix x = y = 0.5 and obtain 3D plots of ∣ ξ exact ∣ , ∣ ξ ∣ , and RPE vs ω x vs ω y for the Lax–Wendroff scheme for

Scenario 1 using k = 0.1 , 0.01.
Scenario 2 using k = 0.01 , 0.001.
Scenario 3 using k = 0.001 , 0.0001.
Scenario 4 using k = 0.001 , 0.0001.
Scenario 5 using k = 0.001 , 0.0001.

6.3 3D plots of ∣ ξ exact ∣ , ∣ ξ ∣ , RPE vs ω x vs ω y

Figures 15, 16, 17, 18, 19, 20, 21, 22, 23, 24 display the plots of modulus of exact amplification factor, modulus of the amplification factor of Lax–Wendroff scheme, RPE of the Lax–Wendroff scheme vs ω x vs ω y . The modulus of exact amplification factor and modulus of amplification factor of Lax–Wendroff are relatively close to each other for the five scenarios considered.

$Figure 15 3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] {\omega }_{x}\in \left[0,\pi ] vs ω y ∈ [ 0 , π ] {\omega }_{y}\in \left[0,\pi ] for the Lax–Wendroff scheme for scenario 1 using k = 0.1 k=0.1 and k = 0.01 k=0.01 . (a) k = 0.1 k=0.1 and (b) k = 0.01 k=0.01 .$

Figure 15

3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] vs ω y ∈ [ 0 , π ] for the Lax–Wendroff scheme for scenario 1 using k = 0.1 and k = 0.01 . (a) k = 0.1 and (b) k = 0.01 .

$Figure 16 3D plots of RPE vs ω x ∈ {\omega }_{x}\hspace{0.33em}\in [0, 1] vs ω y ∈ {\omega }_{y}\hspace{0.33em}\in [0, 1] for the Lax–Wendroff scheme for scenario 1 using k = 0.1 k=0.1 and k = 0.01 k=0.01 . (a) k = 0.1 k=0.1 and (b) k = 0.01 k=0.01 .$

Figure 16

3D plots of RPE vs ω x ∈ [0, 1] vs ω y ∈ [0, 1] for the Lax–Wendroff scheme for scenario 1 using k = 0.1 and k = 0.01 . (a) k = 0.1 and (b) k = 0.01 .

$Figure 17 3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] {\omega }_{x}\in \left[0,\pi ] vs ω y ∈ [ 0 , π ] {\omega }_{y}\in \left[0,\pi ] for the Lax–Wendroff scheme for scenario 2 using k = 0.01 k=0.01 and k = 0.001 k=0.001 . (a) k = 0.01 k=0.01 and (b) k = 0.001 k=0.001 .$

Figure 17

3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] vs ω y ∈ [ 0 , π ] for the Lax–Wendroff scheme for scenario 2 using k = 0.01 and k = 0.001 . (a) k = 0.01 and (b) k = 0.001 .

$Figure 18 3D plots of RPE vs ω x ∈ {\omega }_{x}\hspace{0.33em}\in [0, 1] vs ω y ∈ {\omega }_{y}\hspace{0.33em}\in [0, 1] for the Lax–Wendroff scheme for scenario 2 using k = 0.01 k=0.01 and k = 0.001 k=0.001 . (a) k = 0.01 k=0.01 and (b) k = 0.001 k=0.001 .$

Figure 18

3D plots of RPE vs ω x ∈ [0, 1] vs ω y ∈ [0, 1] for the Lax–Wendroff scheme for scenario 2 using k = 0.01 and k = 0.001 . (a) k = 0.01 and (b) k = 0.001 .

$Figure 19 3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] {\omega }_{x}\in \left[0,\pi ] vs ω y ∈ [ 0 , π ] {\omega }_{y}\in \left[0,\pi ] for the Lax–Wendroff scheme for scenario 3 using k = 0.001 k=0.001 and k = 0.0001 k=0.0001 . (a) k = 0.001 k=0.001 and (b) k = 0.0001 k=0.0001 .$

Figure 19

3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] vs ω y ∈ [ 0 , π ] for the Lax–Wendroff scheme for scenario 3 using k = 0.001 and k = 0.0001 . (a) k = 0.001 and (b) k = 0.0001 .

$Figure 20 3D plots of RPE vs ω x ∈ {\omega }_{x}\hspace{0.33em}\in [0, 1] vs ω y ∈ {\omega }_{y}\hspace{0.33em}\in [0, 1] for the Lax–Wendroff scheme for scenario 3 using k = 0.001 k=0.001 and k = 0.0001 k=0.0001 . (a) k = 0.001 k=0.001 and (b) k = 0.0001 k=0.0001 .$

Figure 20

3D plots of RPE vs ω x ∈ [0, 1] vs ω y ∈ [0, 1] for the Lax–Wendroff scheme for scenario 3 using k = 0.001 and k = 0.0001 . (a) k = 0.001 and (b) k = 0.0001 .

$Figure 21 3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] {\omega }_{x}\in \left[0,\pi ] vs ω y ∈ [ 0 , π ] {\omega }_{y}\in \left[0,\pi ] for the Lax–Wendroff scheme for scenario 4 using k = 0.001 k=0.001 and k = 0.0001 k=0.0001 . (a) k = 0.001 k=0.001 and (b) k = 0.0001 k=0.0001 .$

Figure 21

3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] vs ω y ∈ [ 0 , π ] for the Lax–Wendroff scheme for scenario 4 using k = 0.001 and k = 0.0001 . (a) k = 0.001 and (b) k = 0.0001 .

$Figure 22 3D plots of RPE vs ω x ∈ {\omega }_{x}\hspace{0.33em}\in [0, 1] vs ω y ∈ {\omega }_{y}\hspace{0.33em}\in [0, 1] for the Lax–Wendroff scheme for scenario 4 using k = 0.001 k=0.001 and k = 0.0001 k=0.0001 . (a) k = 0.001 k=0.001 and (b) k = 0.0001 k=0.0001 .$

Figure 22

3D plots of RPE vs ω x ∈ [0, 1] vs ω y ∈ [0, 1] for the Lax–Wendroff scheme for scenario 4 using k = 0.001 and k = 0.0001 . (a) k = 0.001 and (b) k = 0.0001 .

$Figure 23 3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] {\omega }_{x}\in \left[0,\pi ] vs ω y ∈ [ 0 , π ] {\omega }_{y}\in \left[0,\pi ] for the Lax–Wendroff scheme for scenario 5 using k = 0.001 k=0.001 and k = 0.0001 k=0.0001 . (a) k = 0.001 k=0.001 and (b) k = 0.0001 k=0.0001 .$

Figure 23

3D plots of exact amplification factor and amplification factor vs ω x ∈ [ 0 , π ] vs ω y ∈ [ 0 , π ] for the Lax–Wendroff scheme for scenario 5 using k = 0.001 and k = 0.0001 . (a) k = 0.001 and (b) k = 0.0001 .

$Figure 24 3D plots of RPE vs ω x ∈ {\omega }_{x}\hspace{0.33em}\in [0, 1] vs ω y ∈ {\omega }_{y}\hspace{0.33em}\in [0, 1] for the Lax–Wendroff scheme for scenario 5 using k = 0.001 k=0.001 and k = 0.0001 k=0.0001 . (a) k = 0.001 k=0.001 and (b) k = 0.0001 k=0.0001 .$

Figure 24

3D plots of RPE vs ω x ∈ [0, 1] vs ω y ∈ [0, 1] for the Lax–Wendroff scheme for scenario 5 using k = 0.001 and k = 0.0001 . (a) k = 0.001 and (b) k = 0.0001 .

7 Numerical results

All numerical experiments are done in MATLAB platform using Dell core i7 machine.

7.1 Scenario 1

From the analysis of stability for the three methods for scenario 1, we find that the Lax–Wendroff and Du Fort–Frankel schemes are stable when 0 < k ≤ 1.16 and 0 < k ≤ 8.873565 , respectively. We obtain numerical profiles at T = 0.1 and T = 1 using Lax–Wendroff when k = 0.1 and k = 0.01 and using Du Fort–Frankel when k = 1,0.1 and k = 0.01 . We note that for the functional relationship ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 0.5 with Δ x = Δ y = 0.05 and k = 1.3135 × 1 0 − 3 , NSFD is not positively preserving. The contour plots are shown in Figures 25, 26, 27, 28. For scenario 1, Lax–Wendroff and Du Fort–Frankel schemes give quite similar profiles at k = 0.1 and k = 0.01 at times T = 0.1 and T = 1 . The range of the numerical solution is between 1 and 3.5 in these cases. Some dispersive oscillations are seen when Du Fort–Frankel is used at k = 1 at time T = 1 . NSFD is not positive definite for scenario 1 and we observe some overshooting with range of numerical solution being 1–4 at times 0.1 and 1.

$Figure 25 Contour plots of numerical solution vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Du Fort–Frankel scheme for scenario 1 at k = 1 k=1 at time T = 1 T=1 .$

Figure 25

Contour plots of numerical solution vs x vs y using Δ x = Δ y = 0.05 using Du Fort–Frankel scheme for scenario 1 at k = 1 at time T = 1 .

$Figure 26 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff, Du Fort–Frankel schemes for scenario 1 at some values of k k at time T = 0.1 T=0.1 . (a) Lax–Wendroff when k = 0.1 k=0.1 , (b) Lax–Wendroff when k = 0.01 k=0.01 , (c) Du Fort–Frankel when k = 0.1 k=0.1 , (d) Du Fort–Frankel when k = 0.01 k=0.01 .$

Figure 26

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff, Du Fort–Frankel schemes for scenario 1 at some values of k at time T = 0.1 . (a) Lax–Wendroff when k = 0.1 , (b) Lax–Wendroff when k = 0.01 , (c) Du Fort–Frankel when k = 0.1 , (d) Du Fort–Frankel when k = 0.01 .

$Figure 27 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff, Du Fort–Frankel schemes for scenario 1 at some values of k k at time T = 1 T=1 . (a) Lax–Wendroff when k = 0.1 k=0.1 , (b) Lax–Wendroff when k = 0.01 k=0.01 , (c) Du Fort–Frankel when k = 0.1 k=0.1 , (d) Du Fort–Frankel when k = 0.01 k=0.01 .$

Figure 27

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff, Du Fort–Frankel schemes for scenario 1 at some values of k at time T = 1 . (a) Lax–Wendroff when k = 0.1 , (b) Lax–Wendroff when k = 0.01 , (c) Du Fort–Frankel when k = 0.1 , (d) Du Fort–Frankel when k = 0.01 .

$Figure 28 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 1.3135 × 1 0 − 3 k=1.3135\times 1{0}^{-3} using the NSFD scheme for scenario 1 at time T = 0.1 T=0.1 and T = 1 T=1 . (a) NSFD when T = 0.1 T=0.1 and (b) NSFD when T = 1 T=1 .$

Figure 28

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 and k = 1.3135 × 1 0 − 3 using the NSFD scheme for scenario 1 at time T = 0.1 and T = 1 . (a) NSFD when T = 0.1 and (b) NSFD when T = 1 .

7.2 Results for scenario 2

From Section 4, we find that Lax–Wendroff and Du Fort–Frankel schemes are stable when 0 < k ≤ 0.015 and 0 < k ≤ 0.883918 , respectively. We present the profiles at T = 0.1 and T = 1 , using Lax–Wendroff when k = 0.01 and k = 0.001 and using Du Fort–Frankel when k = 0.1,0.01 and k = 0.001 . NSFD is positivity preserving when k = 1.31350 × 1 0 − 3 . The contour plots of the numerical profiles are shown in Figures 29, 30, 31, 32.

$Figure 29 Contour plots of numerical solution vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 for Du Fort–Frankel scheme for scenario 2 when k = 0.1 k=0.1 . (a) Du Fort–Frankel when T = 0.1 T=0.1 , (b) Du Fort–Frankel when T = 1 T=1 .$

Figure 29

Contour plots of numerical solution vs x vs y using Δ x = Δ y = 0.05 for Du Fort–Frankel scheme for scenario 2 when k = 0.1 . (a) Du Fort–Frankel when T = 0.1 , (b) Du Fort–Frankel when T = 1 .

$Figure 30 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 1.3135 × 1 0 − 3 k=1.3135\times 1{0}^{-3} using NSFD schemes for scenario 2 at T = 0.1 T=0.1 and T = 1 T=1 . (a) NSFD when T = 0.1 T=0.1 and (b) NSFD when T = 1 T=1 .$

Figure 30

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 and k = 1.3135 × 1 0 − 3 using NSFD schemes for scenario 2 at T = 0.1 and T = 1 . (a) NSFD when T = 0.1 and (b) NSFD when T = 1 .

$Figure 31 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 2 at some values of k k at time T = 0.1 T=0.1 . (a) Lax–Wendroff when k = 0.01 k=0.01 , (b) Lax–Wendroff when k = 0.001 k=0.001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 31

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 2 at some values of k at time T = 0.1 . (a) Lax–Wendroff when k = 0.01 , (b) Lax–Wendroff when k = 0.001 , (c) Du Fort–Frankel when k = 0.01 , (d) Du Fort–Frankel when k = 0.001 .

$Figure 32 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 2 at some values of k k at time T = 1 T=1 . (a) Lax–Wendroff when k = 0.01 k=0.01 , (b) Lax–Wendroff when k = 0.001 k=0.001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , and (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 32

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 2 at some values of k at time T = 1 . (a) Lax–Wendroff when k = 0.01 , (b) Lax–Wendroff when k = 0.001 , (c) Du Fort–Frankel when k = 0.01 , and (d) Du Fort–Frankel when k = 0.001 .

For scenario 2, NSFD is positive definite and range of numerical solution is 1–3 at time 0.1 and 1–1.25 at time 1. The profiles using Lax–Wendroff and Du Fort–Frankel at k = 0.01 and 0.001 are quite similar at time 0.1 and the range of the numerical solution is 1–2.2. However, at time T = 1 , there are some minor oscillations when Du Fort–Frankel is used with k = 0.01 , 0.001 while the corresponding profiles from Lax–Wendroff are relatively smooth and range of numerical solution is 1–1.1. There are considerable non-physical oscillations when Du Fort–Frankel is used with k = 0.1 at times 0.1 and 1.0 with some values of the numerical solution being negative.

7.3 Results for scenario 3

From the stability analysis, we find that Lax–Wendroff and Du Fort–Frankel are stable when 0 < k ≤ 0.0015 and 0 < k ≤ 0.279510 , respectively. We present the profiles at T = 0.1 , using Lax–Wendroff when k = 0.001 and k = 0.0001 and using Du Fort–Frankel when k = 0.01 and k = 0.001 . NSFD is positivity preserving when k = 1.31350 × 1 0 − 3 . The contour plots of the numerical profiles are shown in Figures 33 and 34. There are massive dispersive oscillation when Du Fort–Frankel is used at k = 0.01 but reasonable profiles are obtained at k = 0.001 .

$Figure 33 Contour plots of numerical solution vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using NSFD scheme when k = 1.3135 × 1 0 − 3 k=1.3135\times 1{0}^{-3} at time T = 0.1 T=0.1 .$

Figure 33

Contour plots of numerical solution vs x vs y using Δ x = Δ y = 0.05 using NSFD scheme when k = 1.3135 × 1 0 − 3 at time T = 0.1 .

$Figure 34 Contour plots of numerical solution vs x x vs y y using Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 for the Lax–Wendroff and Du Fort–Frankel schemes for some values of k k at time T = 0.1 T=0.1 . (a) Lax–Wendroff when k = 0.001 k=0.001 , (b) Lax–Wendroff when k = 0.0001 k=0.0001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , and (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 34

Contour plots of numerical solution vs x vs y using Δ x = Δ y = 0.05 for the Lax–Wendroff and Du Fort–Frankel schemes for some values of k at time T = 0.1 . (a) Lax–Wendroff when k = 0.001 , (b) Lax–Wendroff when k = 0.0001 , (c) Du Fort–Frankel when k = 0.01 , and (d) Du Fort–Frankel when k = 0.001 .

7.4 Results for scenario 4

From the stability analysis, we find that Lax–Wendroff and Du Fort–Frankel are stable when 0 < k ≤ 0.00284 and 0 < k ≤ 0.376892 , respectively. We present the profiles at T = 0.1 , using Lax–Wendroff when k = 0.001 and k = 0.0001 and using Du Fort–Frankel when k = 0.01 and k = 0.001 . NSFD is positivity preserving when k = 1.31350 × 1 0 − 3 . The contour plots of the numerical profiles are shown in Figures 35, 36, to 37. There are massive dispersive oscillations when Du Fort–Frankel is used at k = 0.01 but reasonable profiles are obtained at k = 0.001 . The shape of the circular profile becomes an ellipse as time progresses in scenarios 4 and 5 as in these cases D 1 ≠ D 2 .

$Figure 35 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 1.3135 × 1 0 − 3 k=1.3135\times 1{0}^{-3} using NSFD schemes for scenario 4 at T = 0.1 T=0.1 and T = 1 T=1 . (a) NSFD when T = 0.1 T=0.1 and (b) NSFD when T = 1 T=1 .$

Figure 35

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 and k = 1.3135 × 1 0 − 3 using NSFD schemes for scenario 4 at T = 0.1 and T = 1 . (a) NSFD when T = 0.1 and (b) NSFD when T = 1 .

$Figure 36 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 4 at some values of k k at time T = 0.1 T=0.1 . (a) Lax–Wendroff when k = 0.001 k=0.001 , (b) Lax–Wendroff when k = 0.0001 k=0.0001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , and (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 36

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 4 at some values of k at time T = 0.1 . (a) Lax–Wendroff when k = 0.001 , (b) Lax–Wendroff when k = 0.0001 , (c) Du Fort–Frankel when k = 0.01 , and (d) Du Fort–Frankel when k = 0.001 .

$Figure 37 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 4 at some values of k k at time T = 1 T=1 . (a) Lax–Wendroff when k = 0.001 k=0.001 , (b) Lax–Wendroff when k = 0.0001 k=0.0001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , and (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 37

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 4 at some values of k at time T = 1 . (a) Lax–Wendroff when k = 0.001 , (b) Lax–Wendroff when k = 0.0001 , (c) Du Fort–Frankel when k = 0.01 , and (d) Du Fort–Frankel when k = 0.001 .

7.5 Results for scenario 5

From the stability analysis, we find that Lax–Wendroff and Du Fort–Frankel are stable when 0 < k ≤ 0.00284 and 0 < k ≤ 0.376892 , respectively. We present the profiles at T = 0.1 , using Lax–Wendroff when k = 0.001 and k = 0.0001 and using Du Fort–Frankel when k = 0.01 and k = 0.001 . NSFD is positivity preserving when k = 1.31350 × 1 0 − 3 . The contour plots of the numerical profiles are shown in Figures 38, 39, 40. There are massive dispersive oscillation when Du Fort–Frankel is used at k = 0.01 but reasonable profiles are obtained at k = 0.001 .

$Figure 38 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 and k = 1.3135 × 1 0 − 3 k=1.3135\times 1{0}^{-3} using NSFD schemes for scenario 5 at T = 0.1 T=0.1 and T = 1 T=1 . (a) NSFD when T = 0.1 T=0.1 and (b) NSFD when T = 1 T=1 .$

Figure 38

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 and k = 1.3135 × 1 0 − 3 using NSFD schemes for scenario 5 at T = 0.1 and T = 1 . (a) NSFD when T = 0.1 and (b) NSFD when T = 1 .

$Figure 39 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 5 at some values of k k at time T = 0.1 T=0.1 . (a) Lax–Wendroff when k = 0.001 k=0.001 , (b) Lax–Wendroff when k = 0.0001 k=0.0001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , and (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 39

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 5 at some values of k at time T = 0.1 . (a) Lax–Wendroff when k = 0.001 , (b) Lax–Wendroff when k = 0.0001 , (c) Du Fort–Frankel when k = 0.01 , and (d) Du Fort–Frankel when k = 0.001 .

$Figure 40 Contour plots of numerical solution vs x x vs y y with Δ x = Δ y = 0.05 \Delta x=\Delta y=0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 5 at some values of k k at time T = 1 T=1 . (a) Lax–Wendroff when k = 0.001 k=0.001 , (b) Lax–Wendroff when k = 0.0001 k=0.0001 , (c) Du Fort–Frankel when k = 0.01 k=0.01 , and (d) Du Fort–Frankel when k = 0.001 k=0.001 .$

Figure 40

Contour plots of numerical solution vs x vs y with Δ x = Δ y = 0.05 using Lax–Wendroff and Du Fort–Frankel schemes for scenario 5 at some values of k at time T = 1 . (a) Lax–Wendroff when k = 0.001 , (b) Lax–Wendroff when k = 0.0001 , (c) Du Fort–Frankel when k = 0.01 , and (d) Du Fort–Frankel when k = 0.001 .

8 Conclusion

This work is a major extension of the work by Appadu and Gidey [24] where three numerical schemes, namely, Lax–Wendroff, Du Fort–Frankel and NSFD are used to discretise a 2D nonconstant coefficient advection diffusion equation. A range of values of time-step size of Δ x = Δ y = 0.05 for the three methods at a fixed spatial step size for five scenarios is obtained for which the methods are stable. Numerical results from the three methods are presented for each of the five scenarios where an irregular domain is considered. The problem considered has no known exact solution.

The Du Fort–Frankel scheme has a much wider range of stability than the Lax–Wendroff scheme. However, it causes considerable dispersive oscillations at values of k close to its maximum k for stability and much smaller k must be used to obtain reasonable results. When the values D 1 and D 2 are decreased by factors of 100, 1,000, we observe that the ranges of values of k for stability also decrease by the same factors, in the case of Lax–Wendroff. We are able to construct useful nonstandard finite difference schemes with positivity preserving properties for scenarios 2–5 when the functional relationship is ϕ ( k ) [ ψ ( Δ x ) ] 2 = ϕ ( k ) [ ψ ( Δ y ) ] 2 = 0.5 .

Acknowledgments

The authors are grateful to the two 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 Nelson Mandela University (NMU) for allowing him to use his publication funds to pay for open access fees.
Author contributions: Both authors have accepted responsibility for the entire content of this manuscript and approved its submission. The plan of the work was provided by AR. Derivation was done by AR. Coding, typing of work was done by HG. Both authors were involved in writing up the paper. Work was supervised by AR. 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: All data generated or analysed during this study are included in this published article.

References

[1] Rahaman MM, Takia H, Hasan MK, Hossain MB, Mia S, Hossen K. Application of advection diffusion equation for determination of contaminants in aqueous solution: A mathematical analysis. Appl Math. 2022;10(1):24–31. Search in Google Scholar

[2] Szymkiewicz R. Numerical modeling in open channel hydraulics. vol. 83. Springer Dordrecht; 2010. p. 370.10.1007/978-90-481-3674-2Search in Google Scholar

[3] Kumar N. Unsteady flow against dispersion in finite porous media. J Hydrol. 1983;63(3–4):345–58. 10.1016/0022-1694(83)90050-1Search in Google Scholar

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

[5] Guvanasen V, Volker RE. Numerical solutions for solute transport in unconfined aquifers. Int J Numer Meth Fluids. 1983;3(2):103–23. 10.1002/fld.1650030203Search in Google Scholar

[6] Isenberg J, Gutfinger C. Heat transfer to a draining film. Int J Heat Mass Transfer. 1973;16(2):505–12. 10.1016/0017-9310(73)90075-6Search in Google Scholar

[7] Dehghan M. On the numerical solution of the one-dimensional convection-diffusion equation. Math Problems Eng. 2005;2005:61–74. 10.1155/MPE.2005.61Search in Google Scholar

[8] Dehghan M. Time-splitting procedures for the solution of the two-dimensional transport equation. Kybernetes. 2007;36(5/6):791–805. 10.1108/03684920710749857Search in Google Scholar

[9] Shukla A, Singh AK, Singh P. A recent development of numerical methods for solving convection-diffusion problems. Appl Math. 2011;1(1):1–12. 10.5923/j.am.20110101.01Search in Google Scholar

[10] 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

[11] Appadu AR, Gidey HH. Time-splitting procedures for the numerical solution of the 2D advection-diffusion equation. Math Problems Eng. 2013;2013(1):634657. 10.1155/2013/634657.Search in Google Scholar

[12] 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

[13] 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

[14] Jejeniwa OA, Gidey HH, Appadu AR. Numerical modeling of pollutant transport: Results and optimal parameters. Symmetry. 2022;14(12):2616. 10.3390/sym14122616Search in Google Scholar

[15] Appadu AR, Djoko JK, Gidey HH. Comparative study of some numerical methods to solve a 3D advection-diffusion equation. In: International Conference of Numerical Analysis and Applied Mathematics (ICNAAM 2016). vol. 1863; 2017. p. 030009. 10.1063/1.4992162Search in Google Scholar

[16] Verma AK, Kayenat S. An efficient Mickens’ type NSFD scheme for the generalized Burgers–Huxley equation. J Differ Equ Appl. 2020;26(9–10):1213–46. 10.1080/10236198.2020.1812594Search in Google Scholar

[17] Kayenat S, Verma AK. On the convergence of NSFD schemes for a new class of advection-diffusion-reaction equations. J Differ Equ Appl. 2022;28(7):946–70. 10.1080/10236198.2022.2102425Search in Google Scholar

[18] Kayenat S, Verma AK. Some novel exact and non-standard finite difference schemes for a class of diffusion-advection-reaction equation. J Differ Equ Appl. 2024;30:1–18. 10.1080/10236198.2024.2359106Search in Google Scholar

[19] Kayenat S, Verma AK. NSFD schemes for a class of nonlinear generalised advection-diffusion-reaction equation. Pramana. 2022;96(1):14. 10.1007/s12043-021-02239-1Search in Google Scholar

[20] Strikwerda JC. Finite difference schemes and partial differential equations. Philadelphia: SIAM; 2004. 10.1137/1.9780898717938Search in Google Scholar

[21] El-Nabulsi RA. Neutrons diffusion variable coefficient advection in nuclear reactors. Int J Adv Nucl Reactor Design Tech. 2021;3:102–7. 10.1016/j.jandt.2021.06.005Search in Google Scholar

[22] 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. IOP Publishing; 2019. p. 012009. 10.1088/1742-6596/1180/1/012009Search in Google Scholar

[23] 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 Methods Fluids. 1984;4(9):853–97. 10.1002/fld.1650040905Search in Google Scholar

[24] Appadu AR, Gidey HH. Stability Analysis and Numerical Results for some schemes discretising 2D nonconstant coefficient advection diffusion equations. Open Phys. 2024;22(1):20230195. 10.1515/phys-2023-0195Search in Google Scholar

[25] Corem N, Ditkowski A. New analysis of the Du Fort–Frankel methods. J Scientif Comput. 2012;53:35–54. 10.1007/s10915-012-9627-2Search in Google Scholar

[26] Tannehill JC, Anderson D, Pletcher R. Computational fluid mechanics and heat transfer. 2nd ed. Washington, DC: Taylor & Francis; 1997. Search in Google Scholar

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

[28] Morton KW, Mayers DF. Numerical solution of partial differential equations: an introduction. Cambridge, UK: Cambridge University Press; 2005. 10.1017/CBO9780511812248Search in Google Scholar

Received: 2024-11-26

Revised: 2025-02-14

Accepted: 2025-03-05

Published Online: 2025-04-29

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

Articles in the same Issue

https://doi.org/10.1515/phys-2025-0137

Keywords for this article

advection diffusion; Lax–Wendroff; Du Fort–Frankel; nonstandard finite difference; nonconstant coefficient; stability; irregular domain

Creative Commons

BY 4.0