Some standard and nonstandard finite difference schemes for a reaction–diffusion–chemotaxis model

Gysbert Nicolaas de Waal; Appanah Rao Appadu; Christiaan Johannes Pretorius

doi:10.1515/phys-2022-0231

Article Open Access

Some standard and nonstandard finite difference schemes for a reaction–diffusion–chemotaxis model

Gysbert Nicolaas de Waal , Appanah Rao Appadu and Christiaan Johannes Pretorius

Published/Copyright: March 1, 2023

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 21 Issue 1

Abstract

Two standard and two nonstandard finite difference schemes are constructed to solve a basic reaction–diffusion–chemotaxis model, for which no exact solution is known. The continuous model involves a system of nonlinear coupled partial differential equations subject to some specified initial and boundary conditions. It is not possible to obtain theoretically the stability region of the two standard finite difference schemes. Through running some numerical experiments, we deduce heuristically that these classical methods give reasonable solutions when the temporal step size k is chosen such that k ≤ 0.25 with the spatial step size h fixed at h = 1.0 (first novelty of this work). We observe that the standard finite difference schemes are not always positivity preserving, and this is why we consider nonstandard finite difference schemes. Two nonstandard methods abbreviated as NSFD1 and NSFD2 from Chapwanya et al. are considered. NSFD1 was not used by Chapwanya et al. to generate results for the basic reaction–diffusion–chemotaxis model. We find that NSFD1 preserves positivity of the continuous model if some criteria are satisfied, namely, ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ ≤ 1 2 σ + β and β ≤ σ , and this is the second novelty of this work. Chapwanya et al. modified NSFD1 to obtain NSFD2, which is positivity preserving if R = ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ and 2 σ R ≤ 1 , that is σ ≤ γ , and they presented some results. For the third highlight of this work, we show that NSFD2 is not always consistent and prove that consistency can be achieved if β → 0 and k h 2 → 0 . Fourthly, we show numerically that the rate of convergence in time of the four methods for case 2 is approximately one.

Keywords: standard finite difference method; nonstandard finite difference method; consistency; positivity preserving; cross-diffusion

Nomenclature

conc.: concentration
pop.: population
PDE: partial differential equation
SFD: standard finite difference
NSFD: nonstandard finite difference
k: temporal step size
h: spatial step size
n: time level
m: space level
η ( x , t ): population of bacteria
a ( x , t ): food concentration
N: approximation for population of bacteria
A: approximation for food concentration
σ: diffusion coefficient of the cells
γ: diffusion coefficient of the attractants
β: chemotactic term
λ: spontaneous production coefficient of attractant
ω: decay coefficient of attractant activity
: R = ϕ ( k ) [ ψ ( h ) ] 2

1 Introduction

Diffusion equations are often in the form of reaction–diffusion and advection–diffusion, which have been well studied in the modeling of biological processes [1]. Advection–diffusion models are often used in the field of oceanography to study ocean tracers, such as large-scale ocean circulation [2]. Reaction–diffusion processes can be widely observed in nature, notably in the formation and spread of patterns such as spots and stripes over the surface of animals through the chemical interaction between cells [3]. Unlike reaction–diffusion and advection–diffusion, the mathematical analysis for cross-diffusion equations is a challenge that is largely underdeveloped [1]. Cross-diffusion arises when the concentration gradient of one species induces the flux of another species [4]. Cross-diffusion equations are fundamental in the modeling of several natural processes such as cancer growth [5], population dynamics via Volterra–Lotka cross-diffusion systems [6] and chemotaxis [7]. Theoretically, equations consisting of cross-diffusion terms are challenging largely due to being strongly coupled nonlinear parabolic systems that do not enjoy the maximum principle, and thus, deriving appropriate estimates and proving the existence of positive solutions are not easy [1]. However, in refs [8] and [9], some results on global and local existence of solutions as well as on their long-time behaviour have been established. Recently, some results on the existence and uniqueness of solutions (wellposedness) were obtained in refs [10] and [11]. Standard finite difference (SFD) methods lack the ability to preserve positivity; in contrast, the numerical solutions from nonstandard finite difference (NSFD) methods can preserve properties from the exact solution when solving a differential equation, and this makes NSFD methods appealing [12]. Ronald E. Mickens started to work on NSFD methods around 1990. Another advantage of NSFD over classical methods for reaction–diffusion equations is that classical methods can experience blow up at long propagation time [13,14]. Studying the stability is not easy for classical finite difference methods discretising nonlinear partial differential equations (PDEs). The freezing coefficient technique and von Neumann stability analysis can be used, but freezing coefficient is not a very accurate technique.

NSFD methods are very popular to solve reaction–diffusion PDEs. For reaction–diffusion equations, we can obtain conditions for which NSFD methods are positivity preserving and we can prove boundedness [15]. Positivity and boundedness will ensure stability [15]. It is also difficult to study properties such as stability due to the nature of cross-diffusion systems when SFD methods are used. NSFD schemes are constructed by using some fundamental principles [16]:

The denominator of the discrete derivative must be replaced by a more general function, for example, ∂ u ∂ t ≈ u m n + 1 − u m n ϕ ( k ) , where, for example, ϕ ( k ) = exp ( k ) − 1 .
In general, we use non-local representation of nonlinear terms [17], for example, ( u m ) 2 ≈ u m u m + 1 and ( u m ) 3 ≈ 2 ( u m ) 3 − ( u m ) 2 u m + 1 .
The difference equation should have the same order as the original equation. In general, when the order of the difference equation is larger than the order of the differential equation, spurious solutions might appear as discussed in ref. [18].
The discrete approximation should preserve some important properties of the corresponding differential equation. Properties such as boundedness and positivity should be preserved [19].

If at least one of the first two principles given above is satisfied, we call the scheme a NSFD method [17,20]. Some literature and definitions can be found in ref. [21], and some previous works on reaction–diffusion using NSFD methods can be found in refs [22,23, 24,25].

This article is presented as follows. In Section 2, we describe the basic reaction–diffusion–chemotaxis model to be solved, and we also give some discussion on how the numerical rate of convergence can be calculated in Section 3. Section 4 is dedicated to the two standard schemes SFD1 and SFD2 to solve the reaction–diffusion–chemotaxis model. In Sections 4.3 and 4.4, we present some numerical results. Sections 5 and 6 are dedicated to the two nonstandard schemes NSFD1 and NSFD2 to solve the reaction–diffusion–chemotaxis model. In Section 5.1, we describe the NSFD1 scheme and check the consistency. We present some numerical results in 5.2 for NSFD1. In Section 6.1, we describe NSFD2 scheme [1], and we also check the consistency. In Section 6.2, some numerical results are presented for NSFD2. In Section 7, we give some concluding remarks.

2 The basic reaction–diffusion–chemotaxis model

Chemotaxis refers to the chemically directed movement of a bacterial population η up a gradient in the food a that the bacteria consume.

The mathematical model for the basic reaction–diffusion–chemotaxis model proposed in ref. [7], is described as follows:

(1) ∂ η ∂ t = σ ∂ 2 η ∂ x 2 − β ∂ ∂ x η ∂ a ∂ x ,

(2) ∂ a ∂ t = λ η − ω a + γ ∂ 2 a ∂ x 2 ,

where η ( x , t ) and a ( x , t ) describe the population of bacteria and food concentration, respectively.

The parameters λ , γ , β , ω , and σ are positive constants and γ > σ [7].

Since Eq. (1) contains ∂ η ∂ t on the left-hand side and ∂ ∂ x η ∂ a ∂ x on the right-hand side, the system considered is a cross-diffusion system, the concentration gradient of the bacterial population induces a flux on the concentration gradient of food (attractant).

We consider the following initial conditions given in [1]:

η ( x , 0 ) = exp ( − x 2 ) and a ( x , 0 ) = exp ( − x 2 ) .

We will work with the spatial domain x ∈ [ − 10 , 10 ] and the time domain t ∈ [ 0 , 40 ] . The zero-flux boundary conditions are considered at both the right boundary and the left boundary in this work [26].

The parameters are chosen as follows:

Case 1: λ = γ = β = ω = 1 and σ = 0.5 [1].

Case 2: λ = γ = ω = 1 , σ = 0.5 and β = 0.025 .

We have considered case 2, as one of the conditions for positivity of NSFD1 is that β ≤ σ . Also, to make NSFD2 consistent, we require β → 0 . (More explanation is provided on pages 12, 13, 14, 16 and 17.)

3 Analysis of convergence

As we do not have an analytical solution, we obtain the numerical rate of convergence in time by using

(3) R T = ln E k E k 2 ln ( 2 ) ,

where E k = ∣ ∣ N k − N 2 k ∣ ∣ and E k 2 = ∣ ∣ N k 2 − N k ∣ ∣ are discrete maximum norm errors [22,23].

All numerical simulations are done in MATLAB using an Intel Core i5-10600k with 16GB RAM.

4 Classical methods to solve the reaction–diffusion–chemotaxis model

4.1 SFD1 to solve the continuous model

We construct SFD1 by using forward difference approximation for ∂ η ∂ t , ∂ a ∂ t , ∂ η ∂ x , and ∂ a ∂ x , and centred difference approximation for ∂ 2 a ∂ x 2 and ∂ 2 η ∂ x 2 . By substituting these into Eqs. (1) and (2), we obtain the following system of equations:

(4) N m n + 1 − N m n k = σ N m + 1 n − 2 N m n + N m − 1 n h 2 − β N m + 1 n − N m n h A m + 1 n − A m n h − β N m n A m + 1 n − 2 A m n + A m − 1 n h 2 ,

and

(5) A m n + 1 − A m n k = λ N m n − ω A m n + γ A m + 1 n − 2 A m n + A m − 1 n h 2 ,

where h denotes the spatial step size and k is the temporal step size.

Rewriting Eqs. (4) and (5) yields

(6) N m n + 1 = k h 2 [ N m + 1 n ( σ − β A m + 1 n + β A m n ) ] + N m n − k h 2 [ N m n ( 2 σ − β A m n + β A m − 1 n ) − σ N m − 1 n ]

and

(7) A m n + 1 = λ k N m n − ω k A m n + A m n + γ k h 2 [ A m + 1 n − 2 A m n + A m − 1 n ] .

To check for consistency and to obtain the order of accuracy of SFD1 scheme, we obtain the Taylor series expansion of Eqs. (6) and (7) about the point ( t n , x m ) . By dividing throughout by k and after some rearrangement, we have

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = − k 2 ∂ 2 N ∂ t 2 − k 2 6 ∂ 3 N ∂ t 3 − k 3 24 ∂ 4 N ∂ t 4 + h 2 σ 12 ∂ 4 N ∂ x 4 − h β 2 ∂ N ∂ x ∂ 2 A ∂ x 2 + h 3 ∂ N ∂ x ∂ 3 A ∂ x 3 + h 2 12 ∂ N ∂ x ∂ 4 A ∂ x 4 − h 2 β 12 N ∂ 4 A ∂ x 4 + …

and

∂ A ∂ t − λ N + ω A − γ ∂ 2 A ∂ x 2 = − k 2 ∂ 2 A ∂ t 2 − k 2 6 ∂ 3 A ∂ t 3 − k 3 24 ∂ 4 A ∂ t 4 + h 2 γ 12 ∂ 4 A ∂ x 4 + …

As k , h → 0 , we have

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = 0 ,

and

∂ A ∂ t − λ N + ω A − γ ∂ 2 A ∂ x 2 = 0 .

We therefore conclude that SFD1 scheme is consistent with the PDEs given by Eqs. (1) and (2). The scheme is first-order accurate in time.

4.2 SFD2 to solve the continuous model

We construct SFD2 by using forward difference approximation for ∂ η ∂ t and centred difference approximation for ∂ η ∂ x , ∂ a ∂ x , ∂ 2 a ∂ x 2 , and ∂ 2 η ∂ x 2 . By substituting these into Eq. (1), we obtain SFD2:

(8) N m n + 1 − N m n k = σ N m + 1 n − 2 N m n + N m − 1 n h 2 − β N m + 1 n − N m − 1 n 2 h A m + 1 n − A m − 1 n 2 h − β N m n A m + 1 n − 2 A m n + A m − 1 n h 2 .

Eq. (2) is discretised in a similar way as in SFD1. We therefore have

(9) N m n + 1 = σ k h 2 [ N m + 1 n − 2 N m n + N m − 1 n ] − β k 4 h 2 [ N m + 1 n ( A m + 1 − A m − 1 n ) ] + β k 4 h 2 [ N m − 1 n ( A m + 1 − A m − 1 n ) ] − β k h 2 [ N m n ( A m + 1 − 2 A m n + A m − 1 n ) ] + N m n

and

A m n + 1 = λ k N m n − ω k A m n + A m n + γ k h 2 [ A m + 1 n − 2 A m n + A m − 1 n ] .

We now check for consistency and obtain the order of accuracy of SFD2 scheme. From the Taylor series expansion of Eqs. (7) and (9) about the point ( t n , x m ) , dividing throughout by k and after some rearrangement, we have

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = − k 2 ∂ 2 N ∂ t 2 − k 2 6 ∂ 3 N ∂ t 3 − k 3 24 ∂ 4 N ∂ t 4 + h 2 σ 12 ∂ 4 N ∂ x 4 − h 2 β 12 2 ∂ N ∂ x ∂ 3 A ∂ x 3 + 2 ∂ 3 N ∂ x 3 ∂ A ∂ x + N ∂ 4 A ∂ x 4 − h 4 β 36 ∂ 3 N ∂ x 3 ∂ 3 A ∂ x 3 + …

and

∂ A ∂ t − λ N + ω A − γ ∂ 2 A ∂ x 2 = − k 2 ∂ 2 A ∂ t 2 − k 2 6 ∂ 3 A ∂ t 3 − k 3 24 ∂ 4 A ∂ t 4 + h 2 γ 12 ∂ 4 A ∂ x 4 + …

As k , h → 0 , we have

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = 0 ,

and

∂ A ∂ t − λ N + ω A − γ ∂ 2 A ∂ x 2 = 0 .

We conclude that SFD2 scheme is consistent with the PDEs given by Eqs. (1) and (2). The scheme is first-order accurate in time and second order accurate in space.

4.3 Numerical results for SFD1 and SFD2: Case 1

Stability analysis using the Von Neumann condition cannot be used for the two standard schemes as the equations are coupled. We note that x ∈ [ − 10 , 10 ] and t ∈ [ 0 , 40 ] . We fix h = 1.0 , and observe that the numerical solutions are unbounded when k = 40 80 = 0.5 , as shown in Figure 1. We then keep h = 1.0 fixed and decrease k k = 40 81 , k = 40 82 , k = 40 83 , etc . until we obtain reasonable numerical solutions. We note that reasonable solutions here refer to solutions that are bounded with possibly some dispersive oscillations at some values of x and t . Figure 2 displays numerical solutions using SFD1 with h = 1.0 and k = 0.3008 . We see that the choice of h = 1.0 , k = 40 160 = 0.25 give reasonable numerical solutions as displayed in Figure 3. Figure 4 displays numerical solutions vs x at some values of time t using SFD1 with h = 1.0 , k = 0.25 .

The same approach is used to obtain some values of k and h with reasonable numerical solutions for the SFD2 scheme for case 1. Here also, the combination of h = 1.0 , k = 0.5 gives numerical solutions that are unbounded as seen in Figure 5. We have reasonable numerical solutions using h = 1.0 , k = 0.25 , as depicted in Figure 6. Plots of numerical solutions vs x at some values of time using SFD2 with h = 1.0 , k = 0.25 are displayed in Figure 7.

In Tables 1 and 2, the numerical rate of convergence in time is obtained using SFD1 and SFD2 for case 1, and we conclude that the order of convergence in time is one.

Table 1

E k errors and the numerical rate of convergence in time using SFD1 for case 1 at time t = 1.0 with h = 1.0 and x ∈ [ − 10 , 10 ]

k	E k N	E k A	R T N	R T A
2.5 × 1 0 − 2	8.1638 × 1 0 − 4	1.7850 × 1 0 − 3
1.25 × 1 0 − 2	3.9681 × 1 0 − 4	9.0491 × 1 0 − 4	1.0408	0.9801
6.25 × 1 0 − 3	1.9559 × 1 0 − 4	4.5532 × 1 0 − 4	1.0206	0.9909
3.125 × 1 0 − 3	9.7100 × 1 0 − 5	2.2834 × 1 0 − 4	1.0103	0.9957

Table 2

E k errors and the numerical rate of convergence in time using SFD2 for case 1 at time t = 1.0 with h = 1.0 and x ∈ [ − 10 , 10 ]

k	E k N	E k A	R T N	R T A
2.5 × 1 0 − 2	3.4792 × 1 0 − 4	1.5917 × 1 0 − 3
1.25 × 1 0 − 2	1.6935 × 1 0 − 4	7.9860 × 1 0 − 4	1.0388	0.9950
6.25 × 1 0 − 3	1.9319 × 1 0 − 5	4.0000 × 1 0 − 4	0.8141	0.9975
3.125 × 1 0 − 3	5.1047 × 1 0 − 5	2.0017 × 1 0 − 4	0.9160	0.9988

4.3.1 SFD1 (case 1)

$Figure 1 Results for the chemotaxis model using SFD1 for case 1 with k = 0.5 k=0.5 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , and (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 .$

Figure 1

Results for the chemotaxis model using SFD1 for case 1 with k = 0.5 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , and (c) plot of numerical solutions vs t at x ( 11 ) = 0 .

$Figure 2 Results for the chemotaxis model using SFD1 for case 1 with k = 0.3008 k=0.3008 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 2

Results for the chemotaxis model using SFD1 for case 1 with k = 0.3008 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

$Figure 3 Results for the chemotaxis model using SFD1 for case 1 with k = 0.25 k=0.25 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 3

Results for the chemotaxis model using SFD1 for case 1 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

Figure 4

Plot of numerical solutions vs x at some values of t using SFD1 for case 1 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution vs x at t = 0 , (b) plot of numerical solution vs x at t = 2 , (c) plot of numerical solutions vs t at t = 10 , and (d) plot of numerical solutions vs t at t = 20 .

4.3.2 SFD2 (case 1)

$Figure 5 Results for the chemotaxis model using SFD2 for case 1 with k = 0.5 k=0.5 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , and (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 .$

Figure 5

Results for the chemotaxis model using SFD2 for case 1 with k = 0.5 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , and (c) plot of numerical solutions vs t at x ( 11 ) = 0 .

$Figure 6 Results for the chemotaxis model using SFD2 for case 1 with k = 0.25 k=0.25 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 6

Results for the chemotaxis model using SFD2 for case 1 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

Figure 7

Plot of numerical solutions vs x at some values of t using SFD2 for case 1 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution vs x at t = 0 , (b) plot of numerical solution vs x at t = 2 , (c) plot of numerical solutions vs t at t = 10 , and (d) plot of numerical solutions vs t at t = 20 .

4.3.3 Numerical rate of convergence in time for SFD1 and SFD2: Case 1

4.4 Numerical results for SFD1 and SFD2: case 2

The same approach as discussed in Section 4.3 is used to display numerical solutions using SFD1 and SFD2 for case 2.

We present results in Figures 8, 9 and 10 using SFD1 and SFD2 and we observe that the numerical solutions are bounded and free of dispersive oscillations. We observe no overshoot, unlike those of case 1 (see Figures 3 and 6).

Figures 11 and 12 show the plot of the numerical solutions vs x at four values of t ( t = 0 , 2 , 10 , 20 ) using SFD1 and SFD2 with k = 0.25 and h = 1.0 .

We obtain the numerical rate of convergence in time using SFD1 and SFD2 for case 2 in Tables 3 and 4. We conclude that the order of convergence in time is one for both standard schemes.

Table 3

E k errors and the numerical rate of convergence in time using SFD1 for case 2 at time t = 1.0 with h = 1.0 and x ∈ [ − 10 , 10 ]

k	E k N	E k A	R T N	R T A
2.5 × 1 0 − 2	2.4958 × 1 0 − 3	2.6322 × 1 0 − 3
1.25 × 1 0 − 2	1.2297 × 1 0 − 3	1.3204 × 1 0 − 3	1.0212	0.9953
6.25 × 1 0 − 3	6.1041 × 1 0 − 4	6.6100 × 1 0 − 4	1.0105	0.9983

Table 4

E k errors and the numerical rate of convergence in time using SFD2 for case 2 at time t = 1.0 with h = 1.0 and x ∈ [ − 10 , 10 ]

k	E k N	E k A	R T N	R T A
2.5 × 1 0 − 2	2.4629 × 1 0 − 3	2.6134 × 1 0 − 3
1.25 × 1 0 − 2	1.2137 × 1 0 − 3	1.3113 × 1 0 − 3	1.0210	0.9949
6.25 × 1 0 − 3	6.0246 × 1 0 − 4	6.5653 × 1 0 − 4	1.0104	0.9981

4.4.1 SFD1 (case 2)

$Figure 8 Results for the chemotaxis model using SFD1 for case 2 with k = 0.3008 k=0.3008 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 8

Results for the chemotaxis model using SFD1 for case 2 with k = 0.3008 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

$Figure 9 Results for the chemotaxis model using SFD1 for case 2 with k = 0.25 k=0.25 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 9

Results for the chemotaxis model using SFD1 for case 2 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

$Figure 10 Results for the chemotaxis model using SFD2 for case 2 with k = 0.25 k=0.25 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 10

Results for the chemotaxis model using SFD2 for case 2 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

4.4.2 SFD2 (Case 2)

Figure 11

Plot of numerical solutions vs x at some values of t using SFD1 for case 2 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution vs x at t = 0 , (b) plot of numerical solution vs x at t = 2 , (c) plot of numerical solutions vs t at t = 10 , and (d) plot of numerical solutions vs t at t = 20 .

Figure 12

Plot of numerical solutions vs x at some values of t using SFD2 for case 2 with k = 0.25 and h = 1.0 . (a) Plot of numerical solution vs x at t = 0 , (b) plot of numerical solution vs x at t = 2 , (c) plot of numerical solutions vs t at t = 10 , and (d) plot of numerical solutions vs t at t = 20 .

4.4.3 Numerical rate of convergence in time for SFD1 and SFD2: case 2

5 NSFD method to solve the reaction–diffusion–chemotaxis model

The logistic equation is the space independent case of

(10) ∂ u ∂ t = u ( 1 − u ) − ∂ ∂ x u ∂ c ∂ x ,

which is

d u d t = u ( 1 − u ) .

An exact scheme constructed in [20] to discretise d u d t = u ( 1 − u ) is:

U n + 1 − U n ϕ ( k ) = U n ( 1 − U n + 1 ) ,

where ϕ ( k ) = exp ( k ) − 1 with k → 0 .

5.1 NSFD1

Anguelov et al. [27] discretised the Fisher-Kolmogorov-Petrovsky-Piskunov equation i.e. ∂ u ∂ t = ∂ 2 u ∂ x 2 + λ u ( 1 − u ) using a NSFD scheme for λ = 0.25 , and the scheme they propose is

U m n + 1 − U m n ϕ ( k ) = U m + 1 n − 2 U m n + U m − 1 n [ ψ ( h ) ] 2 + 25 ( 1 − U m n + 1 ) U m − 1 n + U m n + U m + 1 n 3 ,

where ϕ ( k ) = 1 − exp ( − 25 k ) 25 and ψ ( h ) = h .

We use the same discretisation for ∂ u ∂ t and ∂ 2 u ∂ x 2 as given in the study by Anguelov et al. [27]. The same discretisation for ∂ u ∂ x as Chapwanya et al. [1] is used.

This gives the following scheme referred to as NSFD1 for the continuous chemotaxis model:

(11) N m n + 1 − N m n ϕ ( k ) = σ N m + 1 n − 2 N m n + N m − 1 n [ ψ ( h ) ] 2 − β N m n − N m − 1 n ψ ( h ) A m n − A m − 1 n ψ ( h ) − β N m n A m + 1 n − 2 A m n + A m − 1 n [ ψ ( h ) ] 2 ,

and

(12) A m n + 1 − A m n ϕ ( k ) = λ N m n − ω A m n + 1 + γ A m + 1 n − 2 A m n + A m − 1 n [ ψ ( h ) ] 2 ,

where ϕ ( k ) = exp ( k ) − 1 and ψ ( h ) = h .

Eq. (11) can be rewritten as follows:

N m n + 1 − N m n = σ ϕ ( k ) [ ψ ( h ) ] 2 ( N m + 1 n − 2 N m n + N m − 1 n ) − β ϕ ( k ) [ ψ ( h ) ] 2 ( N m n ( A m n − N m n ) − N m − 1 n ( A m n − N m − 1 n ) ) − β ϕ ( k ) [ ψ ( h ) ] 2 N m n ( A m + 1 n − 2 A m n + A m − 1 n ) ,

(13) N m n + 1 = N m + 1 n σ ϕ ( k ) [ ψ ( h ) ] 2 + N m n 1 − 2 σ ϕ ( k ) [ ψ ( h ) ] 2 + β ϕ ( k ) [ ψ ( h ) ] 2 ( A m n − A m + 1 n ) + N m − 1 n σ ϕ ( k ) [ ψ ( h ) ] 2 + β ϕ ( k ) [ ψ ( h ) ] 2 ( A m n − A m − 1 n ) .

We note that NSFD1 was derived by Chapwanya et al. [1] for the chemotaxis model, but they did not use the scheme to present results and also did not study the positivity and boundedness of the method as their aim was to modify NSFD1 to obtain a positivity preserving method.

For simplification, the following relation is used, namely, ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ .

Hence, Eqs. (12) and (13) give

(14) N m n + 1 = σ 2 γ N m + 1 n + 1 − σ γ + β 2 γ ( A m n − A m + 1 n ) N m n + σ 2 γ + β 2 γ ( A m n − A m − 1 n ) N m − 1 n ,

and

(15) A m n + 1 = ϕ ( k ) λ N m n + ( A m + 1 n + A m − 1 n ) ∕ 2 1 + ω ϕ ( k ) .

We note that A m − 1 0 , A m 0 and A m + 1 0 lie from 0 to 1. The scheme given by Eq. (14) preserves positivity of the continuous model if

1 − σ γ − β 2 γ ≥ 0 ⇒ 2 γ ≥ β + 2 σ ,

and

σ 2 γ − β 2 γ ≥ 0 ⇒ σ ≥ β ,

and the initial conditions must be non-negative.

The scheme given by Eq. (15) preserves positivity of the continuous model if the initial conditions are non-negative, i.e., A m 0 ≥ 0 , N m 0 ≥ 0 , and there is no other condition needed.

To check for consistency and to obtain the order of accuracy of NSFD1 scheme, we consider Eqs. (14) and (15) and the Taylor series expansions about the point ( t n , x m ) . After some re-arrangement using k ϕ ( k ) ≈ 1 and ψ ( h ) = h , we obtain

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = − k 2 ∂ 2 N ∂ t 2 − k 2 6 ∂ 3 N ∂ t 3 − k 3 24 ∂ 4 N ∂ t 4 + h β 2 ∂ N ∂ x ∂ 2 A ∂ x 2 + ∂ 2 N ∂ x 2 ∂ A ∂ x + h 2 β 12 σ ∂ 4 N ∂ x 4 − h 2 β 12 N ∂ 4 A ∂ x 4 + 2 ∂ N ∂ x ∂ 3 A ∂ x 3 + 3 ∂ 2 N ∂ x 2 ∂ 2 A ∂ x 2 + 2 ∂ 3 N ∂ x 3 ∂ A ∂ x − h 3 β 24 ∂ N ∂ x ∂ 4 A ∂ x 4 + 2 ∂ 2 N ∂ x 2 ∂ 3 A ∂ x 3 + 3 ∂ 3 N ∂ x 3 ∂ 2 A ∂ x 2 + 2 ∂ 4 N ∂ x 4 ∂ A ∂ x + …

and

∂ A ∂ t − λ N + ω A − γ ∂ 2 A ∂ x 2 = − k 2 ∂ 2 A ∂ t 2 + 2 ω ∂ A ∂ t − k 2 6 ∂ 3 A ∂ t 3 + 3 ω ∂ 2 A ∂ t 2 − k 3 6 ω ∂ 3 A ∂ t 3 − k 4 24 ω ∂ 4 A ∂ t 4 + h 2 γ 12 ∂ 4 A ∂ x 4 + …

As k , h → 0 , we have

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = 0

and

∂ A ∂ t − λ N + ω A − γ ∂ 2 A ∂ x 2 = 0 .

We, therefore, conclude that the NSFD1 scheme is consistent with the PDEs given by Eqs. (1) and (2). The scheme is first-order accurate in time.

5.2 Numerical results for NSFD1

5.2.1 Case 1

In the construction of NSFD1, the relationship ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ is used and since γ = 1 , we have ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 . For NSFD1 to preserve positivity of the continuous model, we showed in Section 5.1 that we need 2 γ ≥ β + 2 σ and σ ≥ β . For case 1, 2 γ ≥ β + 2 σ is satisfied, but σ ≥ β is not satisfied as β = 1 and σ = 0.5 . This means that all the conditions for positivity cannot be met when NSFD1 is used to solve case 1.

Since ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 , we can choose the following combinations of k and h :

h = 1.0 , k = 0.4055 .
h = 0.5 , k = 0.1178 .

We present numerical results in Figures 13 and 14 using NSFD1 for case 1 and observe some overshoot and non-smooth profiles.

$Figure 13 Results for the chemotaxis model using NSFD1 with k = 0.4055 k=0.4055 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 13

Results for the chemotaxis model using NSFD1 with k = 0.4055 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

$Figure 14 Results for the chemotaxis model using NSFD1 with k = 0.1178 k=0.1178 and h = 0.5 h=0.5 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 21 ) = 0 x\left(21)=0 , and (d) zoomed area of subfigure (c).$

Figure 14

Results for the chemotaxis model using NSFD1 with k = 0.1178 and h = 0.5 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 21 ) = 0 , and (d) zoomed area of subfigure (c).

5.2.2 Case 2

For case 2, λ = γ = ω = 1 , σ = 0.5 , and β = 0.025 . Since ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ and γ = 1 , we have ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 . Also, the two conditions for positivity are met, i.e., 2 γ ≥ β + 2 σ and σ ≥ β . This means that all the conditions for positivity are met when NSFD1 is used to solve case 2.

Since ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 , we can choose the following combinations of k and h :

h = 1.0 , k = 0.4055 .
h = 0.5 , k = 0.1178

We present results in Figures 15 and 16 and observe quite smooth profiles as well as all numerical solutions are bounded between 0 and 1 for x ∈ [ − 10 , 10 ] and t ∈ [ 0 , 40 ] , for case 2.

$Figure 15 Results for the chemotaxis model using NSFD1 with k = 0.4055 k=0.4055 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 15

$Figure 16 Results for the chemotaxis model using NSFD1 with k = 0.1178 k=0.1178 and h = 0.5 h=0.5 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 21 ) = 0 x\left(21)=0 , and (d) zoomed area of subfigure (c).$

Figure 16

We verify numerically that the order of convergence in time is one as displayed from calculations in Table 5.

Table 5

E k errors and the numerical rate of convergence in time using NSFD1 for case 2 at time t = 1.0 with h = 2 ( exp ( k ) − 1 ) and x ∈ [ − 10 , 10 ]

k	E k N	E k A	R T N	R T A
1.0137 × 1 0 − 1	2.1687 × 1 0 − 1	1.9108 × 1 0 − 1
5.0683 × 1 0 − 2	1.0718 × 1 0 − 1	8.0873 × 1 0 − 2	1.0168	1.2404
2.5342 × 1 0 − 2	5.5393 × 1 0 − 2	3.9676 × 1 0 − 2	0.9523	1.0274

6 Modification of NSFD1 scheme

6.1 NSFD2

Chapwanya et al. [1] derived NSFD2 by modifying NSFD1.

In NSFD2, N m n + 1 is calculated differently to Eq. (14), but A m n + 1 is obtained similarly as in Eq. (15).

We start with Eq. (11), i.e.,

N m n + 1 − N m n ϕ ( k ) = σ N m + 1 n − 2 N m n + N m − 1 n [ ψ ( h ) ] 2 − β N m n − N m − 1 n ψ ( h ) A m n − A m − 1 n ψ ( h ) − β N m n A m + 1 n − 2 A m n + A m − 1 n [ ψ ( h ) ] 2 ,

which can be rewritten as follows:

N m n + 1 − N m n ϕ ( k ) = σ N m + 1 n − 2 N m n + N m − 1 n [ ψ ( h ) ] 2 − β ( − N m n A m n − N m − 1 n A m n + N m − 1 n A m − 1 n + N m n A m + 1 n ) [ ψ ( h ) ] 2 ,

(16) N m n + 1 − N m n ϕ ( k ) = σ N m + 1 n − 2 N m n + N m − 1 n [ ψ ( h ) ] 2 + β [ ψ ( h ) ] 2 A m n ( N m n + N m − 1 n ) − β [ ψ ( h ) ] 2 ( N m n A m + 1 n + N m − 1 n A m − 1 n ) .

Chapwanya et al. [1] multiplied the negative term of the RHS of Eq. (16) by N m n + 1 N m n .

Eq. (16) is modified to

(17) N m n + 1 − N m n ϕ ( k ) = σ N m + 1 n − 2 N m n + N m − 1 n [ ψ ( h ) ] 2 + β [ ψ ( h ) ] 2 A m n ( N m n + N m − 1 n ) − β [ ψ ( h ) ] 2 ( N m n A m + 1 n + N m − 1 n A m − 1 n ) × N m n + 1 N m n .

The first equation for NSFD2 is therefore

(18) N m n + 1 = N m n ( 1 − 2 σ R ) 1 + β R ( N m n A m + 1 n + N m − 1 n A m − 1 n ) ∕ N m n + σ R ( N m + 1 n + N m − 1 n ) + β R ( N m n + N m − 1 n ) A m n 1 + β R ( N m n A m + 1 n + N m − 1 n A m − 1 n ) ∕ N m n ,

where R = ϕ ( k ) [ ψ ( h ) ] 2 .

From Eq. (12), we have

A m n + 1 − A m n ϕ ( k ) = λ N m n − ω A m n + 1 + γ A m + 1 n − 2 A m n + A m − 1 n [ ψ ( h ) ] 2 ,

and using the functional relation ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ , we obtain

A m n + 1 = ϕ ( k ) λ N m n + ( A m + 1 n + A m − 1 n ) ∕ 2 1 + ω ϕ ( k ) .

Theorem 1

NSFD2 preserves positivity of the continuous model if 2 σ R ≤ 1 :

N m n , A m n ≥ 0 ⇒ N m n + 1 , A m n + 1 ≥ 0 .

Since R = ϕ ( k ) [ ψ ( h ) ] 2 = 1 2 γ , we have that 2 σ R ≤ 1 , which gives 2 σ 2 γ ≤ 1 , i.e., σ ≤ γ .

To check for consistency and to obtain the order of accuracy of NSFD2 scheme, we consider Eq. (18) and the Taylor series expansions about the point ( t n , x m ) . After some rearrangement using k ϕ ( k ) ≈ 1 and R = ϕ ( k ) [ ψ ( h ) ] 2 , we obtain

∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = − k 2 ∂ 2 N ∂ t 2 − k 2 6 ∂ 3 N ∂ t 3 − k 3 24 ∂ 4 N ∂ t 4 + h 2 β ∂ N ∂ x ∂ 2 A ∂ x 2 − h 2 12 β N ∂ 4 A ∂ x 4 + h 3 24 β ∂ N ∂ x ∂ 4 A ∂ x 4 − k h 2 β N ∂ N ∂ t + k 2 ∂ 2 N ∂ t 2 + k 2 6 ∂ 3 N ∂ t 3 + k 3 24 ∂ 4 N ∂ t 4 × 2 N A − h ∂ N ∂ x A + h 2 N ∂ 2 A ∂ x 2 + ∂ N ∂ x ∂ A ∂ x − h 3 2 ∂ N ∂ x ∂ 2 A ∂ x 2 + …

As k , h → 0 , we have

(19) ∂ N ∂ t − σ ∂ 2 N ∂ x 2 + β ∂ N ∂ x ∂ A ∂ x + β N ∂ 2 A ∂ x 2 = − 2 k β h 2 A ∂ N ∂ t .

We therefore conclude that NSFD2 scheme is not consistent with the PDE given by Eq. (1). We would like to mention that the consistency of NSFD2 was not checked in [1]. This is a novel result obtained.

From Eq. (19), we see that NSFD2 can be made consistent if β → 0 and k h 2 → 0 . For this reason, we considered case 2 with β = 0.025 .

Choosing different combinations of k and h result in the following modified equations for case 2:

For k = 0.4055 , h = 1.0 , we have

(20) N t − 0.5 N x x + ( 0.025 ) N x A x + ( 0.025 ) N A x x = − 0.4055 2 N t t − ( 0.4055 ) 2 6 N t t t − ( 0.4055 ) 3 24 N t t t t + 1 80 N x A x x − 1 480 N A x x x x + 1 960 N x A x x x x − 811 8000 1 N × N t + 0.4055 2 N t t + 0.405 5 2 6 N t t t + 0.405 5 3 24 N t t t t × 2 N A − N x ( A − A x + 1 2 A x x ) + N A x x + …

For k = 0.1178 , h = 0.5 , we have

(21) N t − 0.5 N x x + ( 0.025 ) N x A x + ( 0.025 ) N A x x = − 0.1178 2 N t t − 0.117 8 2 6 N t t t − 0.117 8 3 24 N t t t t + 1 160 N x A x x − 1 1920 N A x x x x + 1 7680 N x A x x x x − 589 50000 1 N × N t + 0.1178 2 N t t + 0.117 8 2 6 N t t t + 0.117 8 3 24 N t t t t × 2 N A − N x 1 2 A − 1 4 A x + 1 128 A x x + 1 4 N A x x + …

From Eqs. (20) and (21), we can deduce that the choices k = 0.4055 , h = 1.0 , and k = 0.1178 , h = 0.5 seem reasonable which makes NSFD2 consistent.

6.2 Numerical results for NSFD2

We now present results for the two cases using NSFD2. The profiles from case 1 are not very informative as NSFD2 is not consistent for case 1.

We can choose the following combinations of k and h :

h = 1.0 , k = 0.4055 .
h = 0.5 , k = 0.1178 .

We display plots of numerical solutions using NSFD2 in Figures 17 and 18 for case 1 and Figures 19 and 20 for case 2.

From Figures 17 and 18, we observe some overshoot and uneven profiles when case 1 ( β > σ ) is considered, similar to that of NSFD1. We observe from Figures 19 and 20 that case 2 ( β ≤ σ ) results in numerical solutions that are bounded between 0 and 1 and fairly smooth profiles are obtained.

Table 6 gives the maximum norm errors ( E k ) and the numerical rate of convergence in time ( R T ) for NSFD2. We conclude that NSFD2 is first-order accurate in time for case 2.

Table 6

E k errors and the numerical rate of convergence in time using NSFD2 for case 2 at time t = 1.0 with h = 2 ( exp ( k ) − 1 ) and x ∈ [ − 10 , 10 ]

k	E k N	E k A	R T N	R T A
1.0137 × 1 0 − 1	2.1897 × 1 0 − 1	1.9171 × 1 0 − 1
5.0683 × 1 0 − 2	1.0757 × 1 0 − 1	8.0954 × 1 0 − 2	1.0255	1.2438
2.5342 × 1 0 − 2	5.5645 × 1 0 − 2	3.9689 × 1 0 − 2	0.9509	1.0284

6.2.1 Case 1

$Figure 17 Results for the chemotaxis model using NSFD2 with k = 0.4055 k=0.4055 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 17

Results for the chemotaxis model using NSFD2 with k = 0.4055 and h = 1.0 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 11 ) = 0 , and (d) zoomed area of subfigure (c).

$Figure 18 Results for the chemotaxis model using NSFD2 with k = 0.1178 k=0.1178 and h = 0.5 h=0.5 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 21 ) = 0 x\left(21)=0 , and (d) zoomed area of subfigure (c).$

Figure 18

Results for the chemotaxis model using NSFD2 with k = 0.1178 and h = 0.5 . (a) Plot of numerical solution for food conc. vs x vs t , (b) plot of numerical solution for pop. of bacteria vs x vs t , (c) plot of numerical solutions vs t at x ( 21 ) = 0 , and (d) zoomed area of subfigure (c).

6.2.2 Case 2

$Figure 19 Results for the chemotaxis model using NSFD2 with k = 0.4055 k=0.4055 and h = 1.0 h=1.0 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 11 ) = 0 x\left(11)=0 , and (d) zoomed area of subfigure (c).$

Figure 19

$Figure 20 Results for the chemotaxis model using NSFD2 with k = 0.1178 k=0.1178 and h = 0.5 h=0.5 . (a) Plot of numerical solution for food conc. vs x x vs t t , (b) plot of numerical solution for pop. of bacteria vs x x vs t t , (c) plot of numerical solutions vs t t at x ( 21 ) = 0 x\left(21)=0 , and (d) zoomed area of subfigure (c).$

Figure 20

7 Conclusion

In this work, two SFD methods and two NSFD methods were used to solve a basic reaction–diffusion–chemotaxis model consisting of a cross-diffusion term and a system of nonlinear coupled PDEs requiring positivity preserving numerical solutions. The standard methods SFD1 and SFD2 resulted in unreasonable numerical solutions at some combinations of values of k and h , and it was not possible to study the stability theoretically.

NSFD1 was not always positivity preserving depending on the parameters used. NSFD1 can be used when conditions for positivity hold, that is, when 2 γ ≥ 2 σ + β and σ ≥ β , such as case 2. The classical schemes (SFD1 and SFD2) and the nonstandard schemes (NSFD1 and NSFD2) resulted in more favourable numerical solutions when case 2 ( β ≤ σ ) was considered. NSFD2 was obtained by appropriate modification of NSFD1 and preserved positivity when 2 σ R ≤ 1 . However, we established that NSFD2 is in general not consistent with the PDE describing the population of bacteria given by Eq. (1). We found that a suitable choice of β can render NSFD2 consistent.

Acknowledgments

The authors would like to thank the anonymous reviewers and editor for comments and suggestions that helped to improve the article considerably.

Funding information: G.N. de Waal would like to thank the department of Mathematics and Applied Mathematics of Nelson Mandela University for providing some funding toward his MSc studies. A.R. Appadu is grateful to Nelson Mandela University and the National Research Foundation of South Africa (Grant Number: 136161), which allowed this work to be carried out. C.J. Pretorius is grateful to the National Research Foundation of South Africa (Grant Number: 132544).
Author contributions: The plan of the work was provided by Appadu. Most derivation and analysis was done by de Waal under supervision of Appadu and Pretorius. 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.

References

[1] Chapwanya M, Lubuma JM-S, Mickens RE. Positivity-preserving non-standard finite difference schemes for cross-diffusion equations in biosciences. Comput Math Appl. 2014;68(9):1071–82. 10.1016/j.camwa.2014.04.021. Search in Google Scholar

[2] Roy-Barman M, Jeandel C. Marine geochemistry: ocean circulation, carbon cycle and climate change. Online Edition. Oxford: Oxford Academic; 2016. p. 978-0191829604. 10.1093/acprof:oso/9780198787495.001.0001Search in Google Scholar

[3] Volpert V, Petrovskii S. Reaction-diffusion waves in biology. Phys Life Rev. 2009;6(4):267–310. 10.1016/j.plrev.2009.10.002. Search in Google Scholar PubMed

[4] Vanag, VK, Epstein, IV. Cross-diffusion and pattern formation in reaction–diffusion systems. Phys Chemistry Chem Phys 2009;11(6):897–912. 10.1039/B813825G. Search in Google Scholar PubMed

[5] Marchant BP, Norbury J, Perumpanani AJ. Traveling shock waves arising in a model of malignant invasion. SIAM J Appl Math. 2000;60(2):463–76. 10.1137/S0036139998328034. Search in Google Scholar

[6] Murray JD. Mathematical biology II: spatial models and biomedical applications. 3rd ed. New York: Springer; 2003. p. 978-038795228, 814. 10.1007/b98869Search in Google Scholar

[7] Murray JD. Mathematical biology I: an introduction. 3rd Ed. New York: Springer; 2002. 978-0387952239. p. 551. 10.1007/b98868Search in Google Scholar

[8] Chen L, Jüngel A. Analysis of a parabolic cross-diffusion population model without self-diffusion. J Differ Equ. 2006;60(1):39–59. 10.1016/j.jde.2005.08.002. Search in Google Scholar

[9] Le D. Cross-diffusion equations on n spatial dimensional domains. In: Fifth Mississippi State Conference on Differential Equations and Computational Simulations, Electron. J. Diff. Equ. Conference. vol. 10, 2003, Mississippi. p. 193–210. Search in Google Scholar

[10] Seis, S, Winkler, D. A well-posedness result for a system of cross-diffusion equations. J Evolut Equ. 2020;21:2471–89. 10.1007/S00028-021-00690-6. Search in Google Scholar

[11] Chen, L, Daus, ES, Jüngel, A. Global existence analysis of cross-diffusion population systems for multiple species. Archive Rational Mech Anal. 2018;227(2018):715–47. 10.1007/S00205-017-1172-6. Search in Google Scholar

[12] Mickens RE. Advances in the application of non-standard finite difference schemes. Singapore: World Scientific; 2005. p. 978-9812564047. 664. 10.1142/5884Search in Google Scholar

[13] Sun GF, Liu GR, Li M. An efficient explicit finite-difference scheme for simulating coupled biomass growth on nutritive substrates. Math Probl Eng. 2015;2015:1–17. 10.1155/2015/708497. Search in Google Scholar

[14] Sun GF, Liu GR, Li M. A novel explicit positivity preserving finite-difference scheme for simulating bounded growth of biological films. Int J Comput Meth. 2016;13(2):1640013. 10.1142/S0219876216400132. Search in Google Scholar

[15] Yu Y, Deng W, Wu Y. Positivity and boundedness preserving schemes from space-time fractional predator-prey reaction–diffusion model. Comp Math Appl. 2015;69(8):743–59. 10.1016/j.camwa.2015.02.024. Search in Google Scholar

[16] Mickens RE. Exact solutions to a finite-difference model of a nonlinear reaction-advection equation: implications for numerical analysis. Numer Meth Partial Differ Equ. 1989;5(4):313–25. 10.1002/num.1690050404. Search in Google Scholar

[17] Anguelov R, Lubuma JM-S. Contributions to the mathematics of the non-standard finite difference method and applications. Numer Methods Partial Differ Equ. 2001;17(5):518–43. 10.1002/num.1025. Search in Google Scholar

[18] Hildebrand FB. Finite-difference equations and simulations. New Jersey: Prentice-Hall; 1968. p. 978-0133172300, 338. Search in Google Scholar

[19] Mickens RE. Nonstandard finite difference models of differential equations. Singapore: World Scientific; 1994. p. 978-9810214586, 264. 10.1142/2081Search in Google Scholar

[20] Anguelov R, Lubuma JM-S. Nonstandard finite difference method by non-local approximation. Math Comput Simul. 2003;61(2003):465–75. 10.1016/S0378-4754(02)00106-4. Search in Google Scholar

[21] Mickens RE. Nonstandard finite difference models of differential equations. Singapore: World Scientific; 2000. p. 9978-981-4493-98-7, 664. Search in Google Scholar

[22] Bhatt HP, Khaliq AQM. Fourth-order compact schemes for the numerical simulation of coupled Burgers’ equation. Comput Phys Commun. 2015;200:117–38. 10.1016/j.cpc.2015.11.007Search in Google Scholar

[23] Tijani YO, Appadu AR, Aderogba AA. Some finite difference methods to model biofilm growth and decay: classical and non-standard. Comput. 2021;9(11):123. 10.3390/computation9110123. Search in Google Scholar

[24] Appadu AR, Tijani YO, Aderogba AA. On the performance of some NSFD methods for a 2-D generalized Burgers-Huxley equation. J Differ Equ Appl. 2021;27(11):1537–73. 10.1080/10236198.2021.1999433. Search in Google Scholar

[25] Tijani YO, Appadu AR. Unconditionally positive NSFD and classical finite difference schemes for biofilm formation on medical implant using Allen-Cahn equation. Demonstr Math. 2022;55(1):40–60. 10.1515/dema-2022-0006. Search in Google Scholar

[26] Kovel M, Zubik-Kowal B. Numerical solutions for a model of tissue invasion and migration of tumour cells. Comput Math Meth Med. 2010;2011:1–16. 10.1155/2011/452320. Search in Google Scholar PubMed PubMed Central

[27] Anguelov R, Kama P, Lubuma JM-S. On non-standard finite difference models of reaction–diffusion equations. J Comput Appl Math. 2005;175(1):11–29. 10.1016/j.cam.2004.06.002. Search in Google Scholar

Received: 2022-11-20

Revised: 2023-01-20

Accepted: 2023-01-31

Published Online: 2023-03-01

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

Articles in the same Issue

https://doi.org/10.1515/phys-2022-0231

Keywords for this article

standard finite difference method; nonstandard finite difference method; consistency; positivity preserving; cross-diffusion

Creative Commons

BY 4.0