Next Article in Journal
Operator Smith Algorithm for Coupled Stein Equations from Jump Control Systems
Previous Article in Journal
Computational Approach to Third-Order Nonlinear Boundary Value Problems via Efficient Decomposition Shooting Method
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Developing Higher-Order Unconditionally Positive Finite Difference Methods for the Advection Diffusion Reaction Equations

by
Ndivhuwo Ndou
1,2,*,†,
Phumlani Dlamini
1,† and
Byron Alexander Jacobs
1,†
1
Department of Mathematics and Applied Mathematics, University of Johannesburg, P.O. Box 524, Auckland Park, Johannesburg 2006, South Africa
2
Department of Mathematical and Computational Sciences, University of Venda, Private Bag X5050, Thohoyandou 0950, South Africa
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Axioms 2024, 13(4), 247; https://doi.org/10.3390/axioms13040247
Submission received: 1 February 2024 / Revised: 21 March 2024 / Accepted: 1 April 2024 / Published: 9 April 2024
(This article belongs to the Section Mathematical Analysis)

Abstract

:
This study introduces the higher-order unconditionally positive finite difference (HUPFD) methods to solve the linear, nonlinear, and system of advection–diffusion–reaction (ADR) equations. The stability and consistency of the developed methods are analyzed, which are necessary and sufficient for the numerical approach to converge to the exact solution. The problem under consideration is of the Cauchy type, and hence, Von Neumann stability analysis is used to analyze the stability of the proposed schemes. The HUPFD’s efficacy and efficiency are investigated by calculating the error, convergence rate, and computing time. For validation purposes, the higher-order unconditionally positive finite difference solutions are compared to analytical calculations. The numerical results demonstrate that the proposed methods produce accurate solutions to solve the advection diffusion reaction equations. The results also show that increasing the order of the unconditionally positive finite difference leads an implicit scheme that is conditionally stable and has a higher order of accuracy with respect to time and space.

1. Introduction

The advection–diffusion–reaction (ADR) equation is important in the application of various physical and chemical processes, such as absorption of pollutants in soil, exponential travelling waves, modeling of semiconductors, and modeling of biological processes, among others [1,2,3]. The ADR equation is a partial differential equation generally represented by the equation below:
u t + U a u x D 2 u x 2 = f ( t , x , u ) , ( x , t ) [ a , b ] × [ 0 , T ] t > 0 , u ( 0 , x ) = u 0 ( x ) 0 , u ( a , t ) = u a ( t ) u ( b , t ) = u b ( t ) .
where D is the diffusion coefficient, U a represent the advection coefficient, u x ( t , x ) is the advection term, u x x ( t , x ) is the diffusion term, u t ( t , x ) is the rate at which the concentration of substances changes over time, and f ( t , x , u ) is the reaction term or the source [4]. The advection–diffusion–reaction equation governs the process of advection and diffusion simultaneously [5,6,7,8].
In this paper, we develop the higher-order unconditionally positive finite difference methods (HUPFD) and apply them to solve the ADR equations. These methods extends to the unconditionally positive finite difference method (UPFD) developed by Chen-Charpentier and Kojouharov [1] and it has subsequently been utilized by several other researchers, such as [9,10,11]. The primary advantage of the UPFD is that it guarantees the positivity of solutions regardless of time step or mesh size. Numerical techniques that preserve solution positivity are crucial in physical applications, in quantities such as chemical species concentration, population sizes, and neutron numbers [1] wherein only positive solution have meaning. The ADR equations are often difficult to solve exactly, and hence, the need to develop numerical methods that are good in solving the problem as accurately as possible with less computational time. Benito M et al. studied an unconditionally positivity-preserving scheme for solving the ADR equation. In this study, they found the UPFD schme produces good results when solving the ADR equation as compared to Crank–Nicolson method, and also that the method works for a larger time steps [1]. Appadu A.R. studied the unconditionally positive schemes for biofilm formation on medical implant using the Allen–Cahn equation. In this study, he found that the proposed method produced good results, and that it is stable at all times for all step sizes h , k > 0 [9]. Mahamond et al. tested and improved a nonconventional UPFD method, in which they found that the UPFD scheme is convergent time integrator with order 1, which necessitates investigation of the higher-order UPFD scheme [12]. The study of higher-order numerical methods is of keen interest in the field of fluid dynamics, mathematical physics, engineering, etc, because they are known to produce highly accurate solutions [13,14]. Much research effort is put into developing efficient numerical methods that converge to the solution quickly and with less computational time [15,16,17].
The novel HUPFD schemes we developed are validated on the linear, nonlinear, and system of ADR equations, and its findings are compared with those obtained using the exact solution, UPFD, Crank–Nicolson, and nonstandard finite difference methods. The efficiency of the developed approaches is assessed by comparing the errors, convergence rates, and computing times. In recent years, there has been significant interest in higher-order finite difference approaches. Kovac et al. (2021) introduced a novel approach for solving the diffusion or heat equation using explicit numerical methods. These approaches incorporate Huxley and Fisher’s reaction terms, which proved analytically that the convergence of the approaches for linear ordinary differential systems occurs at a fourth order in the time step size. Additionally, it was found that the provided methods are applicable to solve nonlinear equations in stiff cases. Kovacs et al. (2024) developed a novel two-stage explicit approach for the solution of partial differential equations that involve a diffusion and two reaction terms. A large system with random characteristics and discontinuous initial condition was utilized in this investigation. The researchers demonstrated that in the linear cases, the accuracy of the approaches follows a second order of accuracy, and it is unconditionally stable. Gurarslan et al. (2013) introduced a sixth-order compact difference scheme and a fourth-order Runge–Kutta scheme to generate numerical solutions for the one-dimensional advection–diffusion equation. In this study, the solutions were found to be highly accurate in solving the concurrent transport equation for the one-dimensional advection–diffusion equation. Other examples can be found in [18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33]. The need for the development of HUPFD approaches arise from the limitations observed in numerical methods. In order to achieve accuracy in solutions for solving majority of the partial differential equations, a significant amount of computational time is required.
This study primarily focuses on the development and evaluation of the HUPFD, showing that augmenting the order of the UPFD schemes to higher-order results in an implicit scheme that is unconditionally stable, and progressively improves its accuracy with respect to time and space.
The subsequent sections of the paper are structured sequentially. Section 2 provides an introduction to the UPFD, Section 3 focuses on the development of the HUPFD, and Section 4 and Section 5 include the numerical analysis of the linear advection–diffusion–reaction equation using the HUPFD, which include stability analysis and consistency of the methods. The efficiency and effectiveness of the HUPFD approach in solving the linear advection–diffusion–reaction equations has been shown by numerical examples in Section 6. These examples are presented in comparison to the analytical results. The conclusion of the investigation is presented in Section 7.

2. The Unconditionally Positive Finite Difference Scheme

In this section, we present details of the UPFD schemes [1] as given below:
u x { u i n u i 1 n + 1 Δ x (2a) u i n + 1 u i 1 n Δ x (2b)
u t u i n + 1 u i n Δ t
2 u x 2 u i + 1 n 2 u i n + 1 + u i 1 n ( Δ x ) 2
where the spatial step size is Δ x = b a m 1 and the grid points are given as x i = a + Δ x ( i 1 ) , 1 i m . The time step size Δ t = T N 1 , with the corresponding grid points given as t n = Δ t ( n 1 ) , 1 n N .
Note that when approximating the derivatives of the first and second space, the terms in the finite difference schemes are evaluated at different time levels n and n + 1 . This is required to preserve the positivity of the solution. Equation (2a) is used when the coefficient of the first space derivative is negative, while Equation (2b) is used when the derivative is positive. Conversely, Equation (4) is employed regardless of the positive or negative nature of the coefficient of the second space derivative. However, when the negative coefficient to the diffusion term is used, the solution is expected to blow up, which will undoubtedly impact stability.

3. Development of the Higher-Order Unconditionally Positive Finite Difference Schemes

The development of the higher-order unconditionally positive finite difference (HUPFD) involves the study of three distinct cases, resulting in the formulation of three different algorithms, denoted as HUPFD 1, HUPFD 2, and HUPFD 3. The objective is to develop a robust approach, based on the principles of the UPFD, which preserves the positivity of solutions. The process of developing finite difference schemes for space derivatives in practical situations involves a careful consideration of the tradeoffs between accuracy, stability, consistency, and computing costs. Engineers and scientists employ varying methods in accordance with the specific demands of their simulations and the limitations imposed by their computational resources, with a primary focus on attaining accuracy. This study focuses on the construction of HUPFD schemes utilizing various orders within finite difference schemes for space derivative schemes. The objective is to explore combinations that yield improved results. The unconditionally positive finite difference methods are a category of methods used to approximate the solution of partial differential equations. These methods are applicable to various real-life problems, including heat transfer, stress/strain mechanics, fluid dynamics, and electromagnetics.

3.1. HUPFD 1

In this case, the higher-order method is achieved by considering the first-order formula for the first space derivative and the fourth-order approximation for the second space derivative, as shown below:
u x { u i n u i 1 n + 1 Δ x (5a) u i n + 1 u i 1 n Δ x (5b)
u t u i n + 1 u i n Δ t
2 u x 2 u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2
where the spatial step size is Δ x = b a m 1 and the grid points are given as x i = a + Δ x ( i 1 ) , 2 i m 1 .

3.2. HUPFD 2

In this case, the HUPFD scheme is obtained by considering the second-order formula for the first space derivative and the fourth-order formula for the second space derivative, as shown below:
u x { u i + 1 n + 1 u i 1 n 2 Δ x (8a) u i + 1 n u i 1 n + 1 2 Δ x (8b)
u t u i n + 1 u i n Δ t ,
2 u x 2 u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2 .
where the spatial step size is Δ x = b a m 1 and the grid points are given as x i = a + Δ x ( i 1 ) , 2 i m 1 .

3.3. HUPFD 3

In this case, HUPFD schemes are obtained by considering the fourth-order formula for the derivatives of the first and second space, as shown below:
u x { u i + 2 n + 8 u i + 1 n + 1 8 u i 1 n + u i 2 n + 1 12 Δ x (11a) u i + 2 n + 1 + 8 u i + 1 n 8 u i 1 n + 1 + u i 2 n 12 Δ x (11b)
u t u i n + 1 u i n Δ t
2 u x 2 u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2
where the spatial step size is Δ x = b a m 1 and the grid points are given as x i 1 = a + Δ x ( i 2 ) , 2 i m 1 .

4. Implementation of the Higher-Order Unconditionally Positive Finite Difference Schemes on the Linear ADR Equation

In this section, we set f ( t , x , u ) = u , U a = 1 , and D = 1 in Equation (1) to obtain the equation given below:
u t + u x = 2 u x 2 u ,
subject to the following initial and boundary conditions:
u ( x , 0 ) = e x , 0 x 10 u ( 0 , t ) = e t , 0 t 0.85 u x ( 10 , t ) = u ( 10 , t ) , 0 t 0.85
with the exact solution given by u ( x , t ) = e ( t x ) [1].

4.1. Solving the Linear ADR Equation Using the HUPFD 1 Scheme

Since the coefficient of u x in Equation (14) is positive, we use Equations (5b)–(7) to find the UPFD scheme. Therefore, the HUPFD discretization of Equation (14) is given by the equation below:
u i n + 1 u i n Δ t + u i n + 1 u i 1 n Δ x = u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2 u i n + 1 ,
which simplifies to the following:
1 12 Δ x 2 u i 2 n + 1 + 1 Δ t + 1 Δ x + 30 12 Δ x 2 + 1 u i n + 1 + 1 12 Δ x 2 u i + 2 n + 1 = 1 Δ t u i n + 1 Δ x + 16 12 D x 2 u i 1 2 + 16 12 Δ x 2 u i + 1 n .
Let p 1 = 1 Δ t + 1 Δ x + 30 12 Δ x 2 + 1 , p 2 = 1 12 Δ x 2 , p 3 = 1 2 Δ t , p 4 = 1 Δ x + 16 12 Δ x 2 , thus
p 1 u i n + 1 + p 2 u i + 2 n + 1 + p 2 u i 2 n + 1 = p 3 u i n + p 4 u i 1 n + 16 p 2 u i + 1 n .
Equation (18) is simplified to the following:
A u n + 1 = B u n + d n ,
where:
A u n + 1 = p 1 0 p 2 0 0 0 0 p 1 0 p 2 0 p 2 0 p 1 0 p 2 0 p 2 0 p 1 0 0 0 0 p 2 0 p 1 p 2 0 0 p 2 0 p 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1
B u n = p 3 16 p 2 0 0 0 0 p 4 p 3 16 p 2 0 0 0 p 4 p 3 16 p 2 0 0 0 p 4 p 3 16 p 2 0 0 0 0 p 4 p 3 0 16 p 2 0 0 p 4 p 3 u 1 n u 2 n u 3 n u m 1 n u m n
and:
d n = p 3 u 0 n 0 0 · · · · · · · · u m n p 2 u 0 n + 1 0 0 · · · · · · · · u m + 1 n + 1
where u is the unknown vector of the order m × 1 at any time level t n + 1 and A and B are the coefficient square matrix of order m × m with a triangular structure.

4.1.1. Consistency

We investigate the consistency by using the Taylor expansion on the point ( i , n ) . The Taylor expansions of u i n + 1 , u i n 1 , u i 1 n , u i + 1 n , u i 2 n + 1 , u i + 2 n + 1 are as follows:
u i n + 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + u i n 1 u i n Δ t u t + ( Δ t ) 2 2 2 u t 2 ( Δ t ) 3 6 3 u t 3 + u i 1 n u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + u i + 1 n u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + u i 2 n + 1 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + u i + 2 n + 1 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 +
By substituting Equation (23) into (18), we obtain the following:
p 2 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + + p 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + + p 2 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + = p 3 u i n + p 4 u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + + 16 p 2 u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + .
Simplifying Equation (24) leads to the following equation:
p 1 14 p 2 p 3 p 4 u i n + p 1 Δ t + 2 p 2 Δ t u t + ( p 4 Δ x 16 p 2 Δ x ) Δ x u x + 4 ( Δ x ) 2 p 2 ( Δ x ) 2 2 p 4 2 u x 2 + ( Δ t ) 2 p 2 + ( Δ t ) 2 2 p 1 2 u t 2 = 0 ,
Further simplification of Equation (25) leads to the following equation:
u i n + Δ t Δ x + 32 Δ t 12 ( Δ x ) 2 + Δ t + 1 u t + u x 1 + Δ x 2 2 u x 2 + 31 12 Δ t Δ x 2 + 1 2 ( Δ t ) 2 ( Δ x ) + ( Δ t ) 2 2 + ( Δ t ) 2 2 2 u t 2 = 0 ,
When Δ t 0 and Δ x 0 , Equation (26) is not consistent, and hence, we set Δ t = ( Δ x ) 3 to obtain the following equation:
u i n + ( Δ x ) 2 + 32 12 Δ x + ( Δ x ) 3 + 1 u t + u x 1 + Δ x 2 2 u x 2 + ( Δ t ) 3 2 + ( Δ x ) 5 2 + 30 ( Δ x ) 4 12 + ( Δ x ) 6 2 + ( Δ x ) 2 12 2 u t 2 = 0 ,
When Δ x 0 , Equation (27) leads to the follow equation:
u t + u x = 2 u x 2 u .
Therefore, the original Equation (14) is derived, indicating that the scheme exhibits conditional consistency when Δ t = ( Δ x ) 3 .

4.1.2. Stability

The Von Neumann stability analysis is used to investigate the stability region for the finite difference schemes. First, we define the following equation:
z i n u ˜ i n u i n = ξ n e j Δ t i Δ x ,
where u ˜ i n is exact solution, u i n is numerical solution, ξ is the amplification factor, and ω = Δ t Δ x is a phase angle, where the phase angle ranges from [ π , π ] . The following results are obtained by applying Fourier series analysis to the terms in Equation (18):
u i n = ξ n e j Δ t i Δ x u i 2 n + 1 = ξ n + 1 e j Δ t ( i 2 ) Δ x u i n + 1 = ξ n + 1 e j Δ t ( i ) Δ x u i + 2 n + 1 = ξ n + 1 e j Δ t ( i + 2 ) Δ x u i 1 n = ξ n e j Δ t ( i 1 ) Δ x u i 1 n 1 = ξ n 1 e j Δ t ( i 1 ) Δ x u i + 1 n = ξ n e j Δ t ( i + 1 ) Δ x u i n 1 = ξ n 1 e j Δ t ( i ) Δ x ,
By substituting Equation (30) into (18), we obtain the following:
p 1 ξ n + 1 e j Δ t ( i 2 ) Δ x + p 2 ξ n + 1 e j Δ t ( i ) Δ x + p 1 ξ n + 1 e j Δ t ( i + 2 ) Δ x = ( 16 p 1 + p 4 ) ξ n e j Δ t ( i 1 ) Δ x + p 3 ξ n 1 e j Δ t ( i ) Δ x + 16 p 1 ξ n e j Δ t ( i + 1 ) Δ x ,
By simplifying Equation (31), we obtain the following:
ξ n + 1 p 1 e 2 j Δ Δ x + p 2 + p 1 e 2 j Δ t Δ x = ξ n ( 16 p 1 + p 4 ) e j Δ t Δ x + 16 p 1 e j Δ t Δ x + p 3 ξ n 1 ,
Further simplification to Equation (32) leads to the following equation
ξ 2 p 1 + s 2 e 2 j Δ t Δ x + p 2 e 2 j Δ t Δ x ξ p 3 + p 4 e j Δ t Δ x + 16 p 2 e j Δ t Δ x = 0 ,
Since | ξ | 2 1 from the Von Neuman stability analysis, we have the following:
4 ( 3 Δ x + 8 ) 3 Δ x 2 + 6 Δ x + 16 2 Δ t
The HUPFD scheme is stable, which implies the scheme is convergence since is stable and conditionally consistent when Δ t = ( Δ x ) 3 .

4.2. Solving the Linear ADR Equation Using the HUPFD 2 Scheme

Since the coefficient of u x in Equation (14) is positive, we use Equations (8b)–(10) to find the UPFD scheme. Therefore, the HUPFD discretized Equation (14) is given by the scheme below:
u i n + 1 u i n Δ t + u i + 1 n + 1 u i 1 n 2 Δ x = u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2 u i n + 1 ,
which simplifies to the following:
1 12 Δ x 2 u i 2 n + 1 + 1 Δ t + 30 12 Δ x 2 + 1 u i n + 1 + 1 12 Δ x 2 u i + 2 n + 1 = 1 Δ t u i n + 16 12 Δ x 2 u i + 1 n + u i 1 n + 1 2 Δ x u i 1 n .
Let s 1 = 1 Δ t + 30 12 Δ x 2 + 1 , s 2 = 1 2 Δ x , s 3 = 1 12 Δ x 2 , s 4 = 1 Δ t , and thus:
s 1 u i n + 1 + s 2 u i + 1 n + 1 + s 3 u i 2 n + 1 + s 3 u i + 2 n + 1 = s 3 ( 16 u i + 1 n + 16 u i 1 n ) + s 4 u i n + 16 s 2 u i 1 n .
Equation (36) is simplified to the following equation:
A u n + 1 = B u n + d n ,
where:
A u n + 1 = s 1 s 2 s 3 0 0 0 0 s 1 s 2 s 3 0 s 3 0 s 1 s 2 s 3 0 s 3 0 s 1 0 0 0 0 s 3 0 s 1 s 3 s 2 0 s 3 0 s 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1 ,
B u n = s 4 16 s 3 0 0 0 0 16 s 3 + s 2 s 4 16 s 3 0 0 0 16 s 3 + s 2 s 4 16 s 3 0 0 0 16 s 3 + s 2 s 4 16 s 3 0 0 0 0 16 s 3 + s 2 s 4 0 16 s 3 0 0 16 s 3 + s 2 s 4 u 1 n u 2 n u 3 n u m 1 n u m n ,
and:
d n = s 3 u 0 n 0 0 · · · · · · · · u m n s 2 u 0 n + 1 0 0 · · · · · · · · u m + 1 n + 1
where u is the unknown vector of the order m × 1 at any time level t n + 1 and A and B are the coefficient square matrix of order m × m with a triangular structure.

4.2.1. Consistency

The Taylor expansion of u i n + 1 , u i n 1 , u i 1 n , u i + 1 n , u i 2 n + 1 , u i + 2 n + 1 is obtained from (23). By substituting Equation (23) into (18), we obtain the following:
p 2 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + + p 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + + p 2 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + = p 3 u i n + p 4 u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + + 16 p 2 u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + ,
Simplification of Equation (41) leads to the following equation:
p 1 14 p 2 p 3 p 4 u i n + p 1 Δ t + 2 p 2 Δ t u t + ( p 4 Δ x 16 p 2 Δ x ) Δ x u x + 4 ( Δ x ) 2 p 2 ( Δ x ) 2 2 p 4 2 u x 2 + ( Δ t ) 2 p 2 + ( Δ t ) 2 2 p 1 2 u t 2 = 0 ,
Futher simplification of Equation (42) leads to the following equation:
u i n + Δ t Δ x + 32 Δ t 12 ( Δ x ) 2 + Δ t + 1 u t + u x 1 + Δ x 2 2 u x 2 + 31 12 Δ t Δ x 2 + 1 2 ( Δ t ) 2 ( Δ x ) + ( Δ t ) 2 2 + ( Δ t ) 2 2 2 u t 2 = 0 ,
When Δ t 0 and Δ x 0 , Equation (43) is not consistent, and hence, we set Δ t = ( Δ x ) 3 to obtain the following equation:
u i n + ( Δ x ) 2 + 32 12 Δ x + ( Δ x ) 3 + 1 u t + u x 1 + Δ x 2 2 u x 2 + ( Δ ) 3 2 + ( Δ x ) 5 2 + 30 ( Δ x ) 4 12 + ( Δ x ) 6 2 + ( Δ x ) 2 12 2 u t 2 = 0 ,
When Δ x 0 , Equation (44) leads to the follow equation:
u t + u x = 2 u x 2 u .
Thus, the original Equation (14) is obtained, which implies the scheme is consistent.

4.2.2. Stability

The Von Neuman stability analyses is applied to terms in Equation (37). By substituting Equation (30) into (37), we obtain the following:
ξ n + 1 s 1 + s 2 e j Δ t Δ x + s 3 e 2 j Δ Δ x + s 3 e 2 j Δ t Δ x = ξ n ( 16 s 3 + s 2 ) e j Δ t Δ x + 16 s 3 e j Δ t Δ x + s 3 ,
Since | ξ | 2 1 from the Von Neuman stability analysis, we have the following
1 + 1 3 ( Δ x ) 2 + 1 Δ t
The HUPFD scheme is stable, which implies the scheme is convergence, since it is stable and conditionally consistent when Δ t = ( Δ x ) 3 .

4.3. Solving the Linear ADR Equation Using HUPFD 3 Scheme

In this section, the linear ADR equation is solved by the fourth order UPFD scheme.
By substituting Equations (11b)–(13) into (14), the following equation is obtained:
u i n + 1 u i n Δ t + u i + 2 n + 8 u i + 1 n + 1 8 u i 1 n + u i 2 n + 1 12 Δ x = u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n + u i 2 n + 1 12 Δ x 2 u i n + 1 ,
By simplifying Equation (82), the following is obtained:
1 12 Δ x + 1 12 Δ x 2 u i 2 n + 1 + 1 Δ t + 30 12 Δ x 2 + 1 u i n + 1 + 8 12 Δ x u i + 1 n + 1 + 1 12 Δ x 2 u i + 2 n + 1 = 16 12 Δ x 2 + 8 12 Δ x u i 1 n + 1 Δ t u i n + 16 12 Δ x 2 u i + 1 n + 1 12 Δ x u i + 2 n ,
Let s 1 = 1 12 Δ x + 1 12 Δ x 2 , s 2 = 1 Δ t + 30 12 Δ x 2 + 1 , s 3 = 8 12 Δ x , s 4 = 1 12 Δ x 2 , s 5 = 16 12 Δ x 2 + 8 12 Δ x , s 6 = 16 12 Δ x 2 , s 7 = 1 12 Δ x , s 8 = 1 Δ t .
Equation (83) is simplified to the following equation:
s 1 u i 2 n + 1 + s 2 u i n + 1 + s 3 u i + 1 n + 1 + s 4 u i + 2 n + 1 = s 5 u i 1 n + s 6 u i + 1 n + s 7 u i + 2 n + s 8 u i n ,
The fourth-order UPFD scheme to solve linear ADR equation is as follows:
A u n + 1 = B u n + d
where:
A u n + 1 = s 2 s 3 s 4 0 0 0 0 s 2 s 3 s 4 0 s 1 0 s 2 s 3 s 4 0 s 1 0 s 2 s 3 0 0 0 s 1 0 s 2 s 4 s 3 0 s 1 0 s 2 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1
B u n = s 8 s 6 s 7 0 0 0 s 5 s 8 s 6 s 7 0 0 s 5 s 8 s 6 s 7 0 0 s 5 s 8 s 6 0 0 0 0 s 5 s 8 s 7 s 6 0 0 s 5 s 8 u 1 n u 2 n u 3 n u m 1 n u m n
and:
d n = s 1 u 0 n 0 0 · · · · · · s 4 u m n s 5 u 0 n + 1 0 0 · · · · · · s 7 u m + 1 n + 1
where u is the unknown vector of the order m × 1 at any time level t n + 1 and A and B are the coefficient square matrix of order m × m with a triangular structure.

4.3.1. Consistency

The Taylor expansion of u i n + 1 , u i n 1 , u i 1 n , u i + 1 n , u i 2 n + 1 , u i + 2 n + 1 is obtained from (23). By substituting Equation (23) into (84), we obtain the following:
8 3 Δ x 2 + 1 u i n + 9 Δ t 12 Δ x + 8 Δ t 3 Δ x 2 + Δ t + 1 u t + 7 6 1 6 Δ x u x 2 u x 2 + 3 ( Δ t ) 2 3 Δ x + 21 ( Δ t ) 2 12 ( Δ x ) 2 + ( Δ t ) 2 2 2 u t 2 + Δ t 4 Δ t 12 Δ x + Δ t 12 Δ x 2 2 u t x = 0 .
when Δ t 0 , Δ x 0 , the original reaction-diffusion equation is not produced, which implies the fourth order of unconditionally positive on advection diffusion reaction equation is not consistent.

4.3.2. Stability

The Von Neuman stability analysis is applied to terms in Equation (84). By substituting Equation (30) into (84), we obtain the following:
s 1 ξ n + 1 e j Δ t ( i 2 ) Δ x + s 2 ξ n + 1 e j Δ t ( i ) Δ x + s 3 ξ n + 1 e j Δ t ( i + 1 ) Δ x + s 4 ξ n + 1 e j Δ t ( i + 2 ) Δ x = s 5 ξ n e j Δ t ( i 1 ) Δ x + s 6 ξ n + 1 e j Δ t ( i + 1 ) Δ x + s 7 ξ n e j Δ t ( i + 2 ) Δ x + s 8 ξ n 1 e j Δ t ( i ) Δ x ,
By simplifying Equation (91), the following is obtained:
ξ 2 s 1 e 2 j Δ t Δ x + s 2 + s 3 e j Δ t Δ x + s 4 e 2 j Δ t Δ x ξ s 5 e j Δ t Δ x + s 6 e j Δ t Δ x + s 7 e 2 j Δ t Δ x s 8 = 0 ,
Since | ξ | 2 1 from the Von Neumann stability analysis, the following is obtained:
Δ t 7 8 Δ x 44 3 ( Δ x ) 2 4 ± 7 8 Δ x 44 3 ( Δ x ) 2 4 2 + 3 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 2 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 ,
and:
Δ t 77 24 Δ x + 6 ( Δ x ) 2 + 1 2 ± 77 24 Δ x + 6 ( Δ x ) 2 + 1 2 2 + 3 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 2 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 .
The scheme is not consistent but stable, implying that the scheme is not convergent.

5. Implementation of the Higher-Order Unconditionally Positive Finite Difference Schemes on the Nonlinear ADR Equation

In this section, we set f ( t , x , u ) = r u ( 1 u ) , and U a = 0 in Equation (58), which yields the nonlinear reaction diffusion equation given below:
u t D 2 u x 2 = r u ( 1 u )
with initial and boundary conditions given by:
u ( x , 0 ) = 0.7 e σ ( x x 0 ) 2 , | x x 0 | 0.06 0 , otherwise
u ( 0 , t ) = 0 , 0 t T u ( 1 , t ) = 0 ,
where σ = 80 , x 0 = 0.5 [1].
With the exact solution given as:
u ( t , x ) = 1 1 + 1 u ( x , 0 ) u ( x , 0 ) exp 1 6 D 6 D r x exp 5 6 r t + 1 6 D 6 D r x 2
where D = 0.0002 , r = 0.05 .
The discretization of the nonlinear ADR Equation (58) using Equation (10) is given below:
u i n + 1 u i n 1 2 Δ t D u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2 = r u i n u m a x u i n + 1 ,
By simplifying Equation (61), the following is obtained:
1 2 Δ t + 30 12 ( Δ x ) 2 + r u m a x u i n + 1 + D 12 ( Δ x ) 2 u i + 2 n + 1 + u i 2 n + 1 = D 12 ( Δ x ) 2 16 u i + 1 n + 16 u i 1 n + r u i n + 1 2 Δ t u i n 1 ,
Let s 1 = 1 2 Δ t + 30 12 Δ x + r u m a x , s 2 = D 12 ( Δ x ) 2 , s 3 = 1 2 Δ t , and thus:
s 2 u i 2 n + 1 + s 1 u i n + 1 + s 2 u i + 2 n + 1 = 16 s 2 u i 1 n + r u i n + 16 s 2 u i + 1 n + s 3 u i n 1 ,
Equation (63) is simplified to the following equation:
A u n + 1 = B u n + d n
where:
A u n + 1 = s 1 0 s 2 0 0 0 0 s 1 0 s 2 0 s 2 0 s 1 0 s 2 0 s 2 0 s 1 0 0 0 0 s 2 0 s 1 s 2 0 0 s 2 s 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1 ,
B u n = r 16 s 2 0 0 0 0 16 s 2 r 16 s 2 0 0 0 16 s 2 r 16 s 2 0 0 0 16 s 2 r 16 s 2 0 0 0 0 16 s 2 r 0 16 s 2 0 0 16 s 2 r u 1 n u 2 n u 3 n u m 1 n u m n ,
and:
d = s 3 u 1 n u 2 n u 3 n · · · · · · · u m 1 n u m n .

5.1. Consistency

Consistency of the scheme is investigated by using the Taylor expansion on the point ( m , n ) . By substituting Equation (63) into Equation (23), the following equation is obtained:
s 2 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + + s 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + + s 2 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + = 16 s 2 u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + + r u i n + 16 s 2 u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + + s 3 u i n Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + ,
By simplifying Equation (68), we obtain the following:
s 1 30 s 2 s 3 u i n + 2 Δ t s 2 + Δ t s 1 + Δ t s 3 u t + 12 ( Δ x ) 2 s 2 2 u x 2 + ( Δ t ) 2 s 2 + ( Δ t ) 2 2 s 1 ( Δ t ) 2 s 3 2 u t 2 = 0 ,
Further simplification of Equation (69) leads to the following outcome:
r u m a x u i n + 1 + 2 Δ t D 12 ( Δ x ) 2 + 30 Δ t 23 ( Δ x ) + r Δ t u m a x u t D 2 u x 2 + ( Δ t ) 2 D 12 ( Δ x ) 2 + 30 ( Δ t ) 2 24 ( Δ x ) + r ( Δ t ) 2 u m a x 2 2 u t 2 = 0 ,
When Δ t 0 and Δ x 0 , the scheme becomes:
u t D 2 u x 2 = r u r u m a x u .
Thus, the original PDE is obtained, which implies that the scheme is consistent.

5.2. Stability

The Von Neumann stability analysis will be the one to focus on to obtain the stability region for the finite difference schemes. Firstly, we define:
z m n u ˜ m n u m n = ξ n e j Δ t i Δ x ,
where u ˜ m n is the exact solution, u m n is the numerical solution, ξ is the amplification factor. and ω = Δ t Δ x is a phase angle, where the phase angle ranges from [ π , π ] .
By substituting Equation (72) into (58), the following is obtained:
ξ n + 1 s 2 e 2 j Δ t Δ x + s 1 + s 2 e 2 j Δ t Δ x = ξ n 16 s 2 e j Δ t Δ x + r + 16 s 2 e j Δ t Δ x + ξ n 1 s 3 ,
Equation (73) is simplified to the equation below:
2 s 2 cos ( 2 Δ t Δ x ) + s 1 ξ 2 32 s 2 cos ( Δ t Δ x ) + r ξ s 3 = 0 ,
By using the quadratic Equation ξ = b ± b 2 4 a c 2 a , we obtain the following:
ξ = 32 s 2 cos ( Δ t Δ x ) + r ± 32 s 2 cos ( Δ t Δ x ) + r 2 + 4 s 3 2 s 2 cos ( 2 Δ t Δ x ) + s 1 2 2 s 2 cos ( 2 Δ t Δ x ) + s 1 ,
By using the Von Neumman condition ξ 2 1 to Equation (75), where ξ = ( ( ξ ) ) 2 + ( ( ξ ) ) 2 , the following is obtained:
ξ = 32 s 2 cos ( Δ t Δ x ) + r ± 32 s 2 cos ( Δ t Δ x ) + r 2 + 4 s 3 2 s 2 cos ( 2 Δ t Δ x ) + s 1 2 2 s 2 cos ( 2 Δ t Δ x ) + s 1 ,
By simplifying Equation (76), the following equation is obtained:
ξ 2 = 32 s 2 cos ( Δ t Δ x ) + r ± 32 s 2 cos ( Δ t Δ x ) + r 2 + 4 s 3 2 s 2 cos ( 2 Δ t Δ x ) + s 1 2 2 s 2 cos ( 2 Δ t Δ x ) + s 1 ,
for ω = Δ t Δ x :
ξ 2 = 32 s 2 cos ( ω ) + r ± 32 s 2 cos ( ω ) + r 2 + 4 s 3 2 s 2 cos ( 2 ω ) + s 1 2 2 s 2 cos ( 2 ω ) + s 1 ,
To find the region of stability, we considered an approach followed by Hindmarsh, and Sousa, where two cases are considered; this the case corresponding with ω = π .
Δ t < 32 D 2 ( Δ x ) 2 + 6 r u m a x + 2 r 8 D 3 ( Δ x ) 2 2 r u m a x r 2 r 32 D 12 ( Δ x ) 2 2 ,
and:
Δ t > 55 D 3 ( Δ x ) 2 + 6 r u m a x 8 D ( Δ x ) 2 + 2 r u u m a x 2 r 32 D 12 ( Δ ) 2 2 ,
thus:
55 D 3 ( Δ x ) 2 + 6 r u m a x 8 D ( Δ x ) 2 + 2 r u u m a x 2 r 32 D 12 ( Δ ) 2 2 < Δ t < 32 D 2 ( Δ x ) 2 + 6 r u m a x + 2 r 8 D 3 ( Δ x ) 2 2 r u m a x r 2 r 32 D 12 ( Δ x ) 2 2 .

6. Solving the Nonlinear ADR Equation by Using HUPFD 2

In this section, we solve the linear ADR equation by the HUPFD 2 scheme.
By substituting Equation (10) into (58), the following equation is obtained:
u i n + 1 u i n 1 2 Δ t + u i + 1 n + 8 u i + 1 n + 1 8 u i 1 n 2 + u i 2 n + 1 12 Δ x = u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n + u i 2 n + 1 12 Δ x 2 u i n + 1 ,
By simplifying Equation (82), the following is obtained:
1 12 Δ x + 1 12 Δ x 2 u i 2 n + 1 + 1 2 Δ t + 30 12 Δ x 2 + 1 u i n + 1 + 8 12 Δ x u i + 1 n + 1 + 1 12 Δ x 2 u i + 2 n + 1 = 16 12 Δ x 2 + 8 12 Δ x u i 1 n + 1 2 Δ t u i n 1 + 16 12 Δ x 2 u i + 1 n + 1 12 Δ x u i + 2 n ,
Let s 1 = 1 12 Δ x + 1 12 Δ x 2 , s 2 = 1 2 Δ t + 30 12 Δ x 2 + 1 , s 3 = 8 12 Δ x , s 4 = 1 12 Δ x 2 , s 5 = 16 12 Δ x 2 + 8 12 Δ x , s 6 = 16 12 Δ x 2 , s 7 = 1 12 Δ x , s 8 = 1 2 Δ t .
Equation (83) is simplified to the following equation:
s 1 u i 2 n + 1 + s 2 u i n + 1 + s 3 u i + 1 n + 1 + s 4 u i + 2 n + 1 = s 5 u i 1 n + s 6 u i + 1 n + s 7 u i + 2 n + s 8 u i n 1 ,
The HUPFD 2 scheme to solve nonlinear ADR equation is as follows:
A u n + 1 = B u n + d
where:
A u n + 1 = s 2 s 3 s 4 0 0 0 0 s 2 s 3 s 4 0 s 1 0 s 2 s 3 s 4 0 s 1 0 s 2 s 3 0 0 0 s 1 0 s 2 s 4 s 3 0 s 1 s 2 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1
B u n = 0 s 6 s 7 0 0 0 s 5 0 s 6 s 7 0 0 s 5 0 s 6 s 7 0 0 s 5 0 s 6 0 0 0 0 s 5 0 s 7 s 6 0 0 s 5 0 u 1 n u 2 n u 3 n u m 1 n u m n
and:
d = s 8 u 1 n u 2 n u 3 n · · · · · · u m n .

6.1. Consistency

By using the Taylor expansion given below, the following system of Equations are obtained:
u i n + 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + u i n 1 u i n Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + u i 1 n u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + u i + 1 n u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + u i 2 n + 1 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + u i + 2 n + 1 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + u i + 1 n + 1 u i n + Δ t u t + Δ x u x + ( Δ x ) 2 2 2 u x 2 + 1 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + u i + 2 n u i n + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 8 6 ( Δ x ) 3 3 u x 3 +
By substituting Equation (88) into (84), we obtain the following:
8 3 Δ x 2 + 1 u i n + 9 Δ t 12 Δ x + 8 Δ t 3 Δ x 2 + Δ t + 1 u t + 7 6 1 6 Δ x u x 2 u x 2 + 3 ( Δ t ) 2 3 Δ x + 21 ( Δ t ) 2 12 ( Δ x ) 2 + ( Δ t ) 2 2 2 u t 2 + Δ t 4 Δ t 12 Δ x + Δ t 12 Δ x 2 2 u t x = 0 .
when Δ t 0 , Δ x 0 , the original reaction diffusion equation is not produced, which implies the fourth order of unconditionally positive on advection diffusion reaction equation is not consistent.

6.2. Stability

By applying the Fourier series analyses to the scheme (84), we then obtain the following:
u i n = ξ n e j Δ t i Δ x u i 2 n + 1 = ξ n + 1 e j Δ t ( i 2 ) Δ x u i n + 1 = ξ n + 1 e j Δ t ( i ) Δ x u i + 2 n + 1 = ξ n + 1 e j Δ t ( i + 2 ) Δ x u i 1 n = ξ n e j Δ t ( i 1 ) Δ x u i + 1 n = ξ n + 1 e j Δ t ( i + 1 ) Δ x u i + 2 n = ξ n e j Δ t ( i + 2 ) Δ x u i 1 n 1 = ξ n 1 e j Δ t ( i 1 ) Δ x u i + 1 n = ξ n e j Δ t ( i + 1 ) Δ x u i n 1 = ξ n 1 e j Δ t ( i ) Δ x
By substituting Equation (90) into (84), we obtain the following:
s 1 ξ n + 1 e j Δ t ( i 2 ) Δ x + s 2 ξ n + 1 e j Δ t ( i ) Δ x + s 3 ξ n + 1 e j Δ t ( i + 1 ) Δ x + s 4 ξ n + 1 e j Δ t ( i + 2 ) Δ x = s 5 ξ n e j Δ t ( i 1 ) Δ x + s 6 ξ n + 1 e j Δ t ( i + 1 ) Δ x + s 7 ξ n e j Δ t ( i + 2 ) Δ x + s 8 ξ n 1 e j Δ t ( i ) Δ x ,
By simplifying Equation (91), the following is obtained:
ξ 2 s 1 e 2 j Δ t Δ x + s 2 + s 3 e j Δ t Δ x + s 4 e 2 j Δ t Δ x ξ s 5 e j Δ t Δ x + s 6 e j Δ t Δ x + s 7 e 2 j Δ t Δ x s 8 = 0 ,
Since | ξ | 2 < 1 from the Von Neumann stability analysis, the following is obtained
Δ t < 7 8 Δ x 44 3 ( Δ x ) 2 4 ± 7 8 Δ x 44 3 ( Δ x ) 2 4 2 + 3 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 2 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 ,
and:
Δ t > 77 24 Δ x + 6 ( Δ x ) 2 + 1 2 ± 77 24 Δ x + 6 ( Δ x ) 2 + 1 2 2 + 3 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 2 7 12 Δ x 2 7 4 Δ x + 7 3 ( Δ x ) 2 2 .
The scheme is not consistent, but is stable, and that implies that the scheme is not convergent.

7. Solving Nonlinear ADR by Using HUPFD 3

In this section, we consider solving nonlinear ADR equation by using the HUPFD 3 scheme (82). By substituting (82) into (58), the following is obtained:
u i n + 1 u i n 1 2 Δ t D u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 Δ x 2 = r u i n u m a x u i n + 1 ,
where by using the freezing coefficient technique, u 2 is replace by u u m a x and u is frozen to a constant u m a x , which represents the maximum value of the solution u. By simplifying Equation (95), the following is obtained:
1 2 Δ t + 30 12 ( Δ x ) 2 + r u m a x u i n + 1 + D 12 ( Δ x ) 2 u i + 2 n + 1 + u i 2 n + 1 = D 12 ( Δ x ) 2 16 u i + 1 n + 16 u i 1 n + r u i n + 1 2 Δ t u i n 1 ,
Let s 1 = 1 2 Δ t + 30 12 Δ x + r u m a x , s 2 = D 12 ( Δ x ) 2 , s 3 = 1 2 Δ t , and thus:
s 2 u i 2 n + 1 + s 1 u i n + 1 + s 2 u i + 2 n + 1 = 16 s 2 u i 1 n + r u i n + 16 s 2 u i + 1 n + s 3 u i n 1 ,
The HUPFD 3 scheme to solve the nonlinear ADR equation is as follows:
A u n + 1 = B u n + d n
where:
A u n + 1 = s 1 0 s 2 0 0 0 0 s 1 0 s 2 0 s 2 0 s 1 0 s 2 0 s 2 0 s 1 0 0 0 0 s 2 0 s 1 s 2 0 0 s 2 s 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1 ,
B u n = r 16 s 2 0 0 0 0 16 s 2 r 16 s 2 0 0 0 16 s 2 r 16 s 2 0 0 0 16 s 2 r 16 s 2 0 0 0 0 16 s 2 r 0 16 s 2 0 0 16 s 2 r u 1 n u 2 n u 3 n u m 1 n u m n ,
and:
d = s 3 u 1 n u 2 n u 3 n · · · · · · u m n .

7.1. Consistency

Consistency of the scheme is investigated by using the Taylor expansion on the point ( m , n ) .
u i n + 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + u i n 1 u i n Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + u i 1 n u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + u i + 1 n u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + u i 2 n + 1 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + u i + 2 n + 1 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 +
By substituting Equation (101) into (97), we obtain the following:
s 2 u i n + Δ t u t 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + + s 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + + s 2 u i n + Δ t u t + 2 Δ x u x + 2 ( Δ x ) 2 2 u x 2 + 2 ( Δ t Δ x ) 2 u t x + ( Δ t ) 2 2 2 u t 2 + = 16 s 2 u i n Δ x u x + ( Δ x ) 2 2 2 u x 2 ( Δ x ) 3 6 3 u x 3 + + r u i n + 16 s 2 u i n + Δ x u x + ( Δ x ) 2 2 2 u x 2 + ( Δ x ) 3 6 3 u x 3 + + s 3 u i n Δ t u t + ( Δ t ) 2 2 2 u t 2 + ( Δ t ) 3 6 3 u t 3 + ,
By simplifying Equation (102), the following is obtained:
s 1 30 s 2 s 3 u i n + 2 Δ t s 2 + Δ t s 1 + Δ t s 3 u t + 12 ( Δ x ) 2 s 2 2 u x 2 + ( Δ t ) 2 s 2 + ( Δ t ) 2 2 s 1 ( Δ t ) 2 s 3 2 u t 2 = 0 ,
Further simplification to Equation (103) gives the following outcome:
r u m a x u i n + 1 + 2 Δ t D 12 ( Δ x ) 2 + 30 Δ t 23 ( Δ x ) + r Δ t u m a x u t D 2 u x 2 + ( Δ t ) 2 D 12 ( Δ x ) 2 + 30 ( Δ t ) 2 24 ( Δ x ) + r ( Δ t ) 2 u m a x 2 2 u t 2 = 0 ,
When Δ t 0 and Δ x 0 , the scheme becomes:
u t D 2 u x 2 = r u r u m a x u
Thus, the original PDE equation is obtained, which implies that the scheme is consistent.

7.2. Stability

Von Neumann stability analysis is the focus on to obtain the stability region for the finite difference schemes. Firstly, we define:
z m n u ˜ m n u m n = ξ n e j Δ t i Δ x ,
where u ˜ m n is exact solution, u m n is numerical solution, ξ is the amplification factor, and ω = Δ t Δ x is a phase angle, where the phase angle ranges from [ π , π ]
u i n = ξ n e j Δ t i Δ x u i 2 n + 1 = ξ n + 1 e j Δ t ( i 2 ) Δ x u i n + 1 = ξ n + 1 e j Δ t ( i ) Δ x u i + 2 n + 1 = ξ n + 1 e j Δ t ( i + 2 ) Δ x u i 1 n = ξ n e j Δ t ( i 1 ) Δ x u i 1 n 1 = ξ n 1 e j Δ t ( i 1 ) Δ x u i + 1 n = ξ n e j Δ t ( i + 1 ) Δ x u i n 1 = ξ n 1 e j Δ t ( i ) Δ x
By substituting Equations (107) into (97), the following is obtained:
ξ n + 1 s 2 e 2 j Δ t Δ x + s 1 + s 2 e 2 j Δ t Δ x = ξ n 16 s 2 e j Δ t Δ x + r + 16 s 2 e j Δ t Δ x + ξ n 1 s 3 ,
By simplifying Equation, the following equation is obtained:
ξ n + 1 2 s 2 cos ( 2 Δ t Δ x ) + s 1 = ξ n 32 cos ( Δ t Δ x ) + r + ξ n 1 s 3 ,
By simplifying Equation (109), the following is obtained:
2 s 2 cos ( 2 Δ t Δ x ) + s 1 ξ 2 32 s 2 cos ( Δ t Δ x ) + r ξ s 3 = 0 ,
By using the quadratic Equation ξ = b ± b 2 4 a c 2 a , the following is obtained:
ξ = 32 s 2 cos ( Δ t Δ x ) + r ± 32 s 2 cos ( Δ t Δ x ) + r 2 + 4 s 3 2 s 2 cos ( 2 Δ t Δ x ) + s 1 2 2 s 2 cos ( 2 Δ t Δ x ) + s 1 ,
By using the ξ 2 1 , where ξ = ( ( ξ ) ) 2 + ( ( ξ ) ) 2 , thus:
ξ = 32 s 2 cos ( Δ t Δ x ) + r ± 32 s 2 cos ( Δ t Δ x ) + r 2 + 4 s 3 2 s 2 cos ( 2 Δ t Δ x ) + s 1 2 2 s 2 cos ( 2 Δ t Δ x ) + s 1 ,
Further simplification of (112) leads to the following:
ξ 2 = 32 s 2 cos ( Δ t Δ x ) + r ± 32 s 2 cos ( Δ t Δ x ) + r 2 + 4 s 3 2 s 2 cos ( 2 Δ t Δ x ) + s 1 2 2 s 2 cos ( 2 Δ t Δ x ) + s 1 ,
For ω = Δ t Δ x :
ξ 2 = 32 s 2 cos ( ω ) + r ± 32 s 2 cos ( ω ) + r 2 + 4 s 3 2 s 2 cos ( 2 ω ) + s 1 2 2 s 2 cos ( 2 ω ) + s 1 ,
To find the region of stability, we considered an approach followed by Hindmarsh and Sousa, where two cases are considered, in this the case corresponding with ω = π .
Δ t < 32 D 2 ( Δ x ) 2 + 6 r u m a x + 2 r 8 D 3 ( Δ x ) 2 2 r u m a x r 2 r 32 D 12 ( Δ x ) 2 2 ,
and:
Δ t > 55 D 3 ( Δ x ) 2 + 6 r u m a x 8 D ( Δ x ) 2 + 2 r u u m a x 2 r 32 D 12 ( Δ ) 2 2 ,
thus:
55 D 3 ( Δ x ) 2 + 6 r u m a x 8 D ( Δ x ) 2 + 2 r u u m a x 2 r 32 D 12 ( Δ ) 2 2 < Δ t < 32 D 2 ( Δ x ) 2 + 6 r u m a x + 2 r 8 D 3 ( Δ x ) 2 2 r u m a x r 2 r 32 D 12 ( Δ x ) 2 2 .
The scheme is not consistent, but is stable, which implies that the scheme is not convergent.

8. Implementation of the UPFD and the HUPFD Schemes on the Linear Schnakenberg Model

The linear Schnakenberg model is given by the equations below:
u t + c u x = ( q 1 1 ) u + q 2 v + μ 1 2 u x 2 ,
v t + c v x = q 1 u q 2 v + μ 1 2 v x 2 .
Subject to the following initial and boundary conditions:
u ( x , 0 ) = e x , 0 x 10 u ( 0 , t ) = e t , 0 t 0.85 u x ( 10 , t ) = u ( 10 , t ) , 0 t 0.85
and:
v ( x , 0 ) = e x , 0 x 10 v ( 0 , t ) = e t , 0 t 0.85 v x ( 10 , t ) = v ( 10 , t ) , 0 t 0.85
with constants s = 5 , c = 0 , q 1 = q 2 = 0 , μ 1 = 1 , μ 2 = 5 .

8.1. Solving the Linear Schnakenberg Model Using the UPFD Scheme

The linear Schnakenberg model is solved by using the UPFD. Since the coefficient of u x in Equation (14) is positive, we use Equations (5b)–(7) to find the UPFD scheme. Therefore, the UPFD discretization of Equation (14) is given by the equation below:
u i n + 1 u i n Δ t + c u i n + 1 u i n Δ x = ( q 1 1 ) u i n + 1 + ( q 2 + s ) v i n + μ 1 u i + 1 n 2 u i n + 1 + u i 1 n ( Δ x ) 2 ,
v i n + 1 v i n Δ t + c v i n + 1 v i n Δ x = q 1 u i n + 1 q 2 v i n + 1 + μ 2 v i + 1 n 2 v i n + 1 + v i 1 n ( Δ x ) 2 .
which simplifies to the following:
1 Δ t + c Δ x + ( q 1 1 ) + 2 μ 1 ( Δ x ) 2 u i n + 1 = ( q 2 + s ) v i n + μ 1 Δ x ( u i + 1 n + u i 1 n ) + c Δ x u i n + 1 Δ t u i n , 1 Δ t + c Δ x + q 2 + 2 μ 2 ( Δ x ) 2 u i n + 1 = q 1 u i n + 1 + 1 Δ t v i n + c Δ x v i n + μ 2 ( Δ x ) 2 ( v i + 1 n + v i 1 n ) .
Let s 1 = 1 Δ t + c Δ x + ( q 1 1 ) + 2 μ 1 ( Δ x ) 2 , s 2 = 1 Δ t + c Δ x + q 2 + 2 μ 2 ( Δ x ) 2 , and thus:
u i n + 1 = ( q 2 + s ) s 1 v i n + μ 1 Δ x s 1 ( u i + 1 n + u i 1 n ) + c Δ x s 1 u i n + 1 Δ t s 1 u i n ,
v i n + 1 = q 1 s 2 u i n + 1 + 1 Δ t s 2 v i n + c Δ x s 2 v i n + μ 2 ( Δ x ) 2 s 2 ( v i + 1 n + v i 1 n ) .

8.2. Solving the Schnakenberg Model Using the HUPFD1 Scheme

The Schnakenberg model is solved by using the HUPFD1. Since the coefficient of u x in Equation (14) is positive, we use Equations (5b)–(7) to find the HUPFD1 scheme. Therefore, the HUPFD1 discretization of Equation (14) is given by the equations below:
u i n + 1 u i n Δ t + c u i n + 1 u i 1 n Δ x = ( q 1 1 ) u i n + 1 + ( q 2 + s ) v i n ,
+ μ 1 u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 ( Δ x ) 2 , v i n + 1 v i n Δ t + c v i n + 1 v i 1 n Δ x = q 1 u i n + 1 q 2 v i n + 1 ,
+ μ 2 v i + 2 n + 1 + 16 v i + 1 n 30 v i n + 1 + 16 v i 1 n v i 2 n + 1 12 ( Δ x ) 2 .
which are simplified to the following:
d 1 u i n + 1 + d 2 u i + 2 n + 1 + d 2 u i 2 n + 1 = d 4 u i n + ( q 2 + s ) v i n + 16 d 2 u i + 1 n + d 3 u i 1 n ,
p 1 v i n + 1 + p 2 v i + 2 n + 1 + p 2 v i 2 n + 1 + q 1 u i n + 1 = d 4 v i n + 1 + p 3 v i 1 n + 16 p 2 v i + 1 n .
where a 1 = 1 Δ t + c Δ x + 30 μ 1 12 ( Δ ) 2 , a 2 = μ 1 12 ( Δ x ) 2 , a 3 = c Δ x + 16 μ 1 12 ( Δ x ) 2 , a 4 = 1 Δ t , p 1 = 1 Δ t + c Δ x + q 2 + 30 μ 2 12 ( Δ x ) 2 , p 2 = μ 2 12 ( Δ x ) 2 , p 3 = c Δ x + 16 μ 2 12 ( Δ x ) 2 .
Equations (128) and (129), respectively, are simplified to the following:
A 1 u n + 1 = B 1 u n + d 1 n and A 2 u n + 1 = B 2 u n + d 2 n ,
where:
A 1 u n + 1 = a 1 0 a 2 0 0 0 0 a 1 0 a 2 0 a 2 0 a 1 0 a 2 0 a 2 0 a 1 0 0 0 0 a 2 0 a 1 a 2 0 0 a 2 0 a 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1
B 1 u n = a 4 16 a 2 0 0 0 0 a 3 a 4 16 a 2 0 0 0 a 3 a 4 16 a 2 0 0 0 a 3 a 4 16 a 2 0 0 0 0 d 3 a 4 0 16 a 2 0 0 a 3 a 4 u 1 n u 2 n u 3 n u m 1 n u m n
and:
d 1 n = a 2 u 0 n 0 0 · · · · · · · · u M n a 3 u 0 n + 1 0 0 · · · · · · · · u m + 1 n + 1
and:
A 2 u n + 1 = p 1 0 p 2 0 0 0 0 p 1 0 p 2 0 p 2 0 p 1 0 p 2 0 p 2 0 p 1 0 0 0 0 p 2 0 p 1 p 2 0 0 p 2 0 p 1 v 1 n + 1 v 2 n + 1 v 3 n + 1 v m 1 n + 1 v m n + 1
B 2 u n = d 4 16 p 2 0 0 0 0 p 3 d 4 16 p 2 0 0 0 p 3 d 4 16 p 2 0 0 0 p 3 d 4 16 p 2 0 0 0 0 p 3 d 4 0 16 p 2 0 0 p 3 d 4 v 1 n v 2 n v 3 n v m 1 n v m n
and:
d 2 n = a 2 v 0 n 0 0 · · · · · · · · v m n a 3 u 0 n + 1 0 0 · · · · · · · · v m + 1 n + 1
respectively, where u is the unknown vector of the order m × 1 at any time level t n + 1 and A and B are the coefficient square matrix of order m × m with a triangular structure.

8.3. Solving the Schnakenberg Model Using the HUPFD2 Scheme

The Schnakenberg model is solved by using the HUPFD2. Since the coefficient of u x in Equation (14) is positive, we use Equations (5b)–(7) to find the HUPFD2 scheme. Therefore, the HUPFD2 discretization of Equation (14) is given by the equation below:
u i n + 1 u i n Δ t + c u i + 1 n + 1 u i 1 n Δ x = ( q 1 1 ) u i n + 1 + ( q 2 + s ) v i n
+ μ 1 u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 ( Δ x ) 2 v i n + 1 v i n Δ t + c v i + 1 n + 1 v i 1 n Δ x = q 1 u i n + 1 q 2 v i n + 1
+ μ 2 v i + 2 n + 1 + 16 v i + 1 n 30 v i n + 1 + 16 v i 1 n v i 2 n + 1 12 ( Δ x ) 2
which are simplified to the following:
b 1 u i n + 1 + b 2 u i + 1 n + 1 + b 3 u i + 2 n + 1 + b 3 u i 2 n + 1 = b 4 u i n + b 5 u i 1 n + 16 b 3 u i + 1 n ,
c 1 v i n + 1 + c 2 v i + 1 n + 1 + c 3 v i + 2 n + 1 + c 3 v i 2 n + 1 = c 4 v i n + 1 + c 4 v i 1 n + 16 c 3 v i + 1 n
where b 1 = 1 Δ t ( q 1 1 ) + 30 μ 1 12 ( Δ x ) 2 , b 2 = c 2 ( Δ x ) 2 , b 3 = μ 1 12 ( Δ x ) 2 , b 4 = 1 Δ t , b 5 = c 2 Δ x + 16 μ 1 12 ( Δ x ) 2 , c 1 = 1 Δ t + q 2 + 30 μ 2 12 ( Δ x ) 2 , c 2 = c 2 Δ x , c 3 = μ 2 12 ( Δ x ) 2 , c 4 = c 2 Δ x + 16 μ 2 12 ( Δ x ) 2 .
Equations (139) and (140) are simplified to the following:
A 1 u n + 1 = B 1 u n + d 1 n ,
where:
A 1 u n + 1 = b 1 b 2 b 3 0 0 0 0 b 1 b 2 b 3 0 b 3 0 b 1 b 2 b 3 0 b 3 0 b 1 b 2 0 0 0 b 3 0 b 1 b 3 b 2 0 b 3 0 b 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1
B 1 u n = b 4 16 b 3 0 0 0 0 b 5 b 4 16 b 3 0 0 0 b 5 b 4 16 b 3 0 0 0 b 5 b 4 16 b 3 0 0 0 0 b 5 b 4 0 16 b 3 0 0 b 5 b 4 u 1 n u 2 n u 3 n u m 1 n u m n
and:
d 1 n = b 3 u 0 n 0 0 · · · · · · · · u m n b 3 u 0 n + 1 0 0 · · · · · · · · u m + 1 n + 1
and:
A 2 u n + 1 = c 1 c 2 c 3 0 0 0 0 c 1 c 2 c 3 0 c 3 0 c 1 c 2 c 3 0 c 3 0 c 1 c 2 0 0 0 c 3 0 c 1 c 3 c 2 0 c 3 0 c 1 v 1 n + 1 v 2 n + 1 v 3 n + 1 v m 1 n + 1 v m n + 1
B 2 u n = d 4 16 c 3 0 0 0 0 c 4 d 4 16 c 3 0 0 0 c 4 d 4 16 c 3 0 0 0 c 4 d 4 16 c 3 0 0 0 0 c 4 d 4 0 16 c 3 0 0 c 4 d 4 v 1 n v 2 n v 3 n v m 1 n v m n
and:
d 2 n = c 3 v 0 n 0 0 · · · · · · · · v m n c 3 u 0 n + 1 0 0 · · · · · · · · v m + 1 n + 1
respectively, where u is the unknown vector of the order m × 1 at any time level t n + 1 and A and B are the coefficient square matrix of order m × m with a triangular structure.

8.4. Solving the Schnakenberg Model Using the HUPFD3 Scheme

The Schnakenberg model is solved by using the HUPFD3. Since the coefficient of u x in Equation (14) is positive, we use Equations (5b)–(7) to find the HUPFD3 scheme. Therefore, the HUPFD3 discretization of Equation (14) is given by the equation below:
u i n + 1 u i n Δ t + c u i + 1 n + 8 u i + 1 n + 1 8 u i 1 n + u i 2 n + 1 12 Δ x = ( q 1 1 ) u i n + 1 + ( q 2 + s ) v i n
+ μ 1 u i + 2 n + 1 + 16 u i + 1 n 30 u i n + 1 + 16 u i 1 n u i 2 n + 1 12 ( Δ x ) 2 ,
v i n + 1 v i n Δ t + c v i + 1 n + 8 v i + 1 n + 1 8 v i 1 n + v i 2 n + 1 12 Δ x = q 1 u i n + 1 q 2 v i n + 1
+ μ 2 v i + 2 n + 1 + 16 v i + 1 n 30 v i n + 1 + 16 v i 1 n v i 2 n + 1 12 ( Δ x ) 2 ,
which are simplified to the following:
e 1 u i n + 1 + e 2 u i + 1 n + 1 + e 3 u i 2 n + 1 + e 4 u i + 2 n + 1 = e 5 u i n + q 5 v i n + e 5 u i + 2 n + ( 8 e 3 + 16 e 4 ) u i 1 n + 16 e 4 u i + 1 n ,
f 1 v i n + 1 + q 1 u i n + 1 + 8 f 2 v i + 1 n + 1 + f 3 v i + 2 n + 1 + f 4 v i 2 n + 1 = e 5 v i n + 16 f 3 v i + 1 n + f 2 v i + 2 n + f 5 v i 1 n ,
where e 1 = 1 Δ t ( q 1 1 ) + 30 μ 1 12 ( Δ x ) 2 , e 2 = 8 c 12 Δ x , e 3 = c 12 Δ + μ 1 12 ( Δ x ) 2 , e 4 = μ 1 12 ( Δ x ) 2 , e 5 = 1 Δ t   f 1 = 1 Δ t + q 2 + 30 μ 2 12 ( Δ x ) 2 , f 2 = c 12 Δ x , f 3 = μ 2 12 ( Δ x ) 2 , f 4 = c 12 Δ x + μ 2 12 ( Δ x ) 2 , f 5 = 8 c 12 Δ x + 16 μ 2 12 ( Δ x ) 2 .
Equations (150) and (151) are simplified to the following:
A 1 u n + 1 = B 1 u n + d 1 n ,
where:
A 1 u n + 1 = e 1 e 2 e 4 0 0 0 0 e 1 e 2 e 4 0 e 3 0 e 1 e 2 e 4 0 e 3 0 e 1 e 4 0 0 0 e 3 0 e 1 e 4 e 2 0 e 3 0 e 1 u 1 n + 1 u 2 n + 1 u 3 n + 1 u m 1 n + 1 u m n + 1
B 1 u n = e 5 16 e 4 0 0 0 0 8 e 3 + 16 e 4 e 5 16 e 4 0 0 0 8 e 3 + 16 e 4 e 5 16 e 4 0 0 0 8 e 3 + 16 e 4 e 5 16 e 4 0 0 0 0 8 e 3 + 16 e 4 e 5 0 16 e 4 0 0 8 e 3 + 16 e 4 e 5 u 1 n u 2 n u 3 n u m 1 n u m n
and:
d 1 n = e 3 u 0 n 0 0 · · · · · · · · u m n e 3 u 0 n + 1 0 0 · · · · · · · · u m + 1 n + 1
and:
A 2 u n + 1 = f 1 8 f 2 f 3 0 0 0 f 4 f 1 8 f 2 f 3 0 0 f 4 f 1 8 f 2 f 3 0 0 f 4 f 1 8 f 2 0 0 0 0 f 4 f 1 f 3 8 f 2 0 0 f 4 f 1 v 1 n + 1 v 2 n + 1 v 3 n + 1 v m 1 n + 1 v m n + 1
B 2 u n = e 5 16 f 3 f 2 0 0 0 f 5 e 5 16 f 3 f 2 0 0 f 5 e 5 16 f 3 f 2 0 0 f 5 e 5 16 f 3 0 0 0 0 f 5 e 5 f 2 16 f 3 0 0 f 5 e 5 v 1 n v 2 n v 3 n v m 1 n v m n
and:
d 2 n = f 3 v 0 n 0 0 · · · · · · · · v m n f 3 u 0 n + 1 0 0 · · · · · · · · v m + 1 n + 1
respectively, where u is the unknown vector of the order m × 1 at any time level t n + 1 and A and B are the coefficient square matrix of order m × m with a triangular structure.
Theorem 1 
(Solution positivity). Let u 0 ( x ) 0 , 0 g 0 ( x ) 1 , v 0 ( x ) 0 for the numerical scheme (18), (37), (124), (125), (128), (129), (139), (140), (150), and (151); we have u i n + 1 0 , 0 g i n + 1 , v i n + 1 0 , independent of the mesh step size [11], where g is the reaction term.
Proof. 
Let us assume that u i n + 1 0 , 0 g i n + 1 1 , v i n 0 , for i = 0 , 1 , , m , n = 0 , m .
Equations (19), (38), (51), (130), (141), and (152) is composed of a coefficient matrix with strictly diagonally dominant tridiagonal matrix, with positive values on diagonal and subdiagonals, which implies nonnegative inverse. Thus, the solution to the algebraic system is nonnegative u i n + 1 0 and v i n + 1 0 , i = 0 , 1 , m .
Furthermore, it is observed that g i n + 1 0 in Equations (19), (38), (51), (130), (141), and (152). By applying the same consideration at each time level, we complete the proof. □

9. Numerical Results

In this section, we present solutions obtained by the positivity-preserving higher-order methods for solving the linear, nonlinear ADR equations and Schnakenberg system of equations. The validity of the method’s outcomes is verified by a comparative analysis with the exact, Crank Nicolson, and NSFD solutions. The performance of the approaches is evaluated by calculating the convergence rates, absolute errors, and computing time. The numerical solution is represented by u ¯ i , j , while the analytical solution is denoted by u i , j . The error is represented by the symbol e i , j , and its magnitude is quantified by the L -norm, which is determined by the following formula:
| | e i , n | | = | | u i , n u ¯ i , n | | = M a x | u i , n u ¯ i , n |
The convergence rate is calculated by using the formula below:
q = l o g 2 e k e k + 1 l o g 2 Δ t Δ t k + 1

9.1. Numerical Results for the Linear Advection–Diffusion–Reaction Equation

Table 1 presents the convergence rates of the UPFD, Crank–Nicolson, NSFD, and HUPFD solutions with respect to the time and spatial steps. The results show that the HUPFD1 approach shows a faster convergence to the solution compared to the UPFD, HUPFD2, and HUPFD3 approaches. In Table 2, solutions obtained by the exact, UPFD, Crank–Nicolson, NSFD, and HUPFD 1, HUPFD 2, and HUPFD 3 solutions at ( t , x ) = 1.33334 , 0.0002 , 2.6668 , 0.0004 , and 4.0002 , 0.0006 are compared. Table 3 presents a comparison of the infinite norm results obtained from the HUPDF methods developed in this study, with respect to the Crank Nicolson and NSFD methods, along with their respective computational times. The results show that the UPFD approach exhibits faster convergence to the solution compared to the HUPFD techniques. This is primarily attributed to the fact that the HUPFD methods require an increased number of nodal points for approximating solutions to the problem.
The infinity norms of the proposed HUPFD methods are compared with the UPFD, Crank–Nicolson, and NSFD methods in solving the linear ADR equation, as shown in Table 3. In this example, the results demonstrate that HUPFD approaches preserve the accuracy of the solution. The surface plots of the UPFD, HUPFD1, and exact solutions of linear advection–diffusion–reaction equation are shown in Figure 1a to Figure 2a.

9.2. Numerical Results for the Nonlinear Advection Diffusion Equation

Table 4 presents the convergence rates of the UPFD and the HUPFD with respect to the time and spatial steps. The convergence rate of HUPFD is slightly higher compared to that of UPFD. The comparison of the exact, UPFD, and HUPFD solutions at ( x , t ) = 0.3 , 0.06 , 0.4 , 0.008 , and 0.5 , 0.1 is presented in Table 5. Additionally, the computational time for the UPFD and HUPFD methods are shown. The accuracy of the HUPFD approach is slightly greater in comparison to the UPFD method. However, the UPFD method exhibits low convergence rate as compared to the HUPFD method when solving the linear ADR equation.
The infinity norms of the proposed HUPFD are compared with the UPFD and NSFD methods in solving the linear ADR equation, as shown in Table 6. The findings indicate that the HUPFD approach achieved slightly high accuracy (Figure 3 and Figure 4).
Figure 5 show the absolute errors between the exact solution and HUPFD solution for the nonlinear ADR against time.

9.3. Numerical Results for the Schnakenberg Model

Table 7 presents the convergence rates of the UPFD, HUPFD 1, HUPFD 2, and HUPFD 3 in relation to the time and spatial steps. The findings show that the HUPFD methods provide a slightly higher convergence rate compared to the UPFD, HUPFD2, and HUPFD 3 approaches while solving the Schnakenberg model. The solutions of the UPFD, HUPFD 1, HUPFD 2, and HUPFD 3 solutions at ( x , t ) = 1.33334 , 0.0002 , 2.6668 , 0.0004 , and 4.0002 , 0.0006 are compared in Table 8 and Table 9. Additionally, the computing times are also provided. The results show that the HUPFD 1 approach has less computational time in comparison to the UPFD, HUPFD 2, and HUPFD 3 solutions when solving the problem. The surface plots of the UPFD, HUPFD1, HUPFD2, and HUPFD3 of the Schnakenberg model are shown in Figure 6a to Figure 7b. The figure demonstrates that the implemented approaches achieve positivity of solution when solving the problem.

10. Conclusions

The HUPFD scheme was developed and used in this study to solve the linear, nonlinear, and system of ADR equations. The goal of this study was to investigate the performance of the HUPFD in various aspects, including stability analysis, consistency, convergence rate, error, computational time, and accuracy. The findings show that the HUPFD performs well as an approximation method to solve the ADR equations, thereby confirming its efficiency and effectiveness. However, it is important to note that the consistency of the HUPFD is conditional. Additionally, we noted that augmenting the order of the UPFD scheme results in an implicit scheme that exhibits conditional stability and enhanced accuracy in terms of time and space. Nevertheless, these schemes also possess certain limitations, including requiring more computational resources than explicit methods. This study has the potential to be extended to a two-dimensional case in order to replicate real-world challenges. It is possible to build higher-order schemes for the time and space fractional advection diffusion equations, which could result in a numerical method that is highly efficient.

Author Contributions

Formal analysis, N.N.; Investigation, N.N.; Supervision, P.D. and B.A.J.; Writing—original draft, N.N.; Writing—review editing, P.D. and B.A.J. All authors have read and agreed to the published version of the manuscript.

Funding

The authors wish to acknowledge the financial support from the University of Johannesburg and the University of Venda.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare that there are no conflicts of interests regarding the publication of this paper.

Abbreviations

HUPFDHigher-order unconditionally positive finite difference methods
ADRAdvection–diffusion–reaction
PDEsPartial differential equations
CPU timeComputational time

Nomenclature

xSpatial variable
tTime variable
Δ x Change in spatial variable
Δ t Change in time variable
u t Rate at which concentration of substances changes over time
u x Advection term
u x x Diffusion term
UAdvection coefficient
DDiffusion coefficient

References

  1. Chen-Charpentier, B.; Kojouharov, H. An unconditionally positivity preserving scheme for advection—Diffusion reaction equations. Math. Comput. Model. 2013, 57, 2177–2185. [Google Scholar] [CrossRef]
  2. Bohn, H. Soil absorption of air pollutants. J. Environ. Qual. 1972, 1, 372–377. [Google Scholar] [CrossRef]
  3. Guo, Y. Exponential stability analysis of travelling waves solutions for nonlinear delayed cellular neural networks. Dyn. Syst. 2017, 32, 490–503. [Google Scholar] [CrossRef]
  4. Britton, N. Others Reaction-Diffusion Equations and Their Applications to Biology; Academic Press: Cambridge, MA, USA, 1986. [Google Scholar]
  5. Hao, X.; Liu, L.; Wu, Y. Iterative solution for nonlinear impulsive advection-reaction-diffusion equations. J. Nonlinear Sci. Appl. 2016, 9, 4070–4077. [Google Scholar] [CrossRef]
  6. Liu, B.; Allen, M.; Kojouharov, H.; Chen, B. Others Finite-element solution of reaction-diffusion equations with advection. Comput. Methods Water Resour. XI 1996, 1, 3–12. [Google Scholar]
  7. Chapwanya, M.; Lubuma, J.; Mickens, R. Nonstandard finite difference schemes for Michaelis-Menten type reaction-diffusion equations. Numer. Methods Partial. Differ. Equ. 2013, 29, 337–360. [Google Scholar] [CrossRef]
  8. Khan, H.; Gómez-Aguilar, J.; Khan, A.; Khan, T. Stability analysis for fractional order advection–reaction diffusion system. Phys. A Stat. Mech. Its Appl. 2019, 521, 737–751. [Google Scholar] [CrossRef]
  9. Appadu, A. Analysis of the unconditionally positive finite difference scheme for advection-diffusion-reaction equations with different regimes. AIP Conf. Proc. 2016, 1738, 030005. [Google Scholar]
  10. Tijani, Y.; Appadu, A. Unconditionally positive NSFD and classical finite difference schemes for biofilm formation on medical implant using Allen-Cahn equation. Demonstr. Math. 2022, 55, 40–60. [Google Scholar] [CrossRef]
  11. Ndou, N.; Dlamini, P.; Jacobs, B. Enhanced Unconditionally Positive Finite Difference Method for Advection-Diffusion-Reaction Equations. Mathematics 2022, 10, 2639. [Google Scholar] [CrossRef]
  12. Saleh, M.; Kovács, E.; Pszota, G. Testing and improving a non-conventional unconditionally positive finite difference method. Multidiszcip. TudomáNyok Miskolci Egy. KöZleméNye 2020, 10, 206–213. [Google Scholar] [CrossRef]
  13. Deville, M.; Fischer, P.; Mund, E. High-Order Methods for Incompressible Fluid Flow; Cambridge University Press: Cambridge, UK, 2002. [Google Scholar]
  14. Gassner, G.; Beck, A. On the accuracy of high-order discretizations for underresolved turbulence simulations. Theor. Comput. Fluid Dyn. 2013, 27, 221–237. [Google Scholar] [CrossRef]
  15. Bini, D.; Latouche, G.; Meini, B. A family of fast fixed point iterations for M/G/1-type Markov chains. IMA J. Numer. Anal. 2022, 42, 1454–1477. [Google Scholar] [CrossRef]
  16. Southworth, B.; Krzysik, O.; Pazner, W.; Sterck, H. Fast Solution of Fully Implicit Runge–Kutta and Discontinuous Galerkin in Time for Numerical PDEs, Part I: The Linear Setting. SIAM J. Sci. Comput. 2022, 44, A416–A443. [Google Scholar] [CrossRef]
  17. Chen, J.; Huang, Z. On a fast deterministic block Kaczmarz method for solving large-scale linear systems. Numer. Algorithms 2022, 89, 1007–1029. [Google Scholar] [CrossRef]
  18. Gurarslan, G.; Karahan, H.; Alkaya, D.; Sari, M.; Yasar, M. Numerical solution of advection-diffusion equation using a sixth-order compact finite difference method. Math. Probl. Eng. 2013, 2013, 672936. [Google Scholar] [CrossRef]
  19. Cui, M. Compact finite difference method for the fractional diffusion equation. J. Comput. Phys. 2009, 228, 7792–7804. [Google Scholar] [CrossRef]
  20. Tian, Z.; Ge, Y. A fourth-order compact ADI method for solving two-dimensional unsteady convection—Diffusion problems. J. Comput. Appl. Math. 2007, 198, 268–286. [Google Scholar] [CrossRef]
  21. Li, Y.; Kim, J. An efficient and stable compact fourth-order finite difference scheme for the phase field crystal equation. Comput. Methods Appl. Mech. Eng. 2017, 319, 194–216. [Google Scholar] [CrossRef]
  22. Li, L.; Jiang, Z.; Yin, Z. Fourth-order compact finite difference method for solving two-dimensional convection—Diffusion equation. Adv. Differ. Equ. 2018, 2018, 234. [Google Scholar] [CrossRef]
  23. Visbal, M.; Gaitonde, D. On the use of higher-order finite-difference schemes on curvilinear and deforming meshes. J. Comput. Phys. 2002, 181, 155–185. [Google Scholar] [CrossRef]
  24. Carpenter, M.; Gottlieb, D.; Abarbanel, S. The stability of numerical boundary treatments for compact high-order finite-difference schemes. J. Comput. Phys. 1993, 108, 272–295. [Google Scholar] [CrossRef]
  25. Vasilyev, O. High order finite difference schemes on non-uniform meshes with good conservation properties. J. Comput. Phys. 2000, 157, 746–761. [Google Scholar] [CrossRef]
  26. Chung, Y.; Tucker, P. Accuracy of higher-order finite difference schemes on nonuniform grids. AIAA J. 2003, 41, 1609–1611. [Google Scholar] [CrossRef]
  27. Ndou, N.; Dlamini, P.; Jacobs, B. Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method. Mathematics 2024, 12, 1009. [Google Scholar]
  28. Calvo, M.; De Frutos, J.; Novo, J. Linearly implicit Runge–Kutta methods for advection–reaction–diffusion equations. Appl. Numer. Math. 2001, 37, 535–549. [Google Scholar] [CrossRef]
  29. Zhao, G.; Ruan, S. Time periodic traveling wave solutions for periodic advection–reaction–diffusion systems. J. Differ. Equations 2014, 257, 1078–1147. [Google Scholar] [CrossRef]
  30. Dwivedi, K.; Das, S.; Baleanu, D. Numerical solution of nonlinear space–time fractional-order advection–reaction–diffusion equation. J. Comput. Nonlinear Dyn. 2020, 15, 061005. [Google Scholar] [CrossRef]
  31. Bennequin, D.; Gander, M.; Gouarin, L.; Halpern, L. Optimized Schwarz waveform relaxation for advection reaction diffusion equations in two dimensions. Numer. Math. 2016, 134, 513–567. [Google Scholar] [CrossRef]
  32. Sari, M.; Gürarslan, G.; Zeytinoğlu, A. High-order finite difference schemes for solving the advection-diffusion equation. Math. Comput. Appl. 2010, 15, 449–460. [Google Scholar] [CrossRef]
  33. Kovács, E.; Majár, J.; Saleh, M. Unconditionally Positive, Explicit, Fourth Order Method for the Diffusion-and Nagumo-Type Diffusion–Reaction Equations. J. Sci. Comput. 2024, 98, 39. [Google Scholar] [CrossRef]
Figure 1. (a) UPFD solution of the nonlinear ADR equation. (b) HUPFD solution of the linear ADR equation.
Figure 1. (a) UPFD solution of the nonlinear ADR equation. (b) HUPFD solution of the linear ADR equation.
Axioms 13 00247 g001
Figure 2. (a) Exact solution of the linear ADR equation. (b) Absolute error obtained by the HUPFD method.
Figure 2. (a) Exact solution of the linear ADR equation. (b) Absolute error obtained by the HUPFD method.
Axioms 13 00247 g002
Figure 3. (a) Analytical solution of the nonlinear ADR equation. (b) HUPFD solution of the nonlinear ADR equation.
Figure 3. (a) Analytical solution of the nonlinear ADR equation. (b) HUPFD solution of the nonlinear ADR equation.
Axioms 13 00247 g003
Figure 4. (a) UPFD solution of the nonlinear ADR equation. (b) Comparison of the analytical and numerical solutions of the nonlinear ADR equation.
Figure 4. (a) UPFD solution of the nonlinear ADR equation. (b) Comparison of the analytical and numerical solutions of the nonlinear ADR equation.
Axioms 13 00247 g004
Figure 5. Absolute difference between analytical and HUPFD of the nonlinear ADR equation.
Figure 5. Absolute difference between analytical and HUPFD of the nonlinear ADR equation.
Axioms 13 00247 g005
Figure 6. (a) UPFD solution of the Schnakenberg system equations. (b) HUPFD1 solution of the Schnakenberg system equations.
Figure 6. (a) UPFD solution of the Schnakenberg system equations. (b) HUPFD1 solution of the Schnakenberg system equations.
Axioms 13 00247 g006
Figure 7. (a) HUPFD2 solution of the Schnakenberg system equations. (b) HUPFD3 solution of the Schnakenberg system equations.
Figure 7. (a) HUPFD2 solution of the Schnakenberg system equations. (b) HUPFD3 solution of the Schnakenberg system equations.
Axioms 13 00247 g007
Table 1. Convergence rates for the linear ADR equation by the UPFD, NSFD, Crank–Nicolson, and HUPFD with respect to varying Δ x and Δ t .
Table 1. Convergence rates for the linear ADR equation by the UPFD, NSFD, Crank–Nicolson, and HUPFD with respect to varying Δ x and Δ t .
q ( Δ t k )
Δ t k Δ x k UPFDHUPFD 1HUPFD 2HUPFD 3Crank–NicolsonNSFD
0.00005 0.6667 2.1675 2.5233 1.6978 1.0006 4.5606 2.0579
0.0001 0.6667 0.9973 1.0054 1.1996 1.0212 1.0051 1.9427
0.0002 0.6667 1.8201 1.9990 1.5640 1.0023 2.7372 1.8871
0.0004 0.6667
Table 2. Numerical results of the linear ADR equation at a fixed Δ x = 0.6667 and Δ t = 0.0001 .
Table 2. Numerical results of the linear ADR equation at a fixed Δ x = 0.6667 and Δ t = 0.0001 .
( x , t ) Exact SolutionUPFD SolutionAbsolute Error (Exact and UPFD)HUPFD 1Absolute Error (Exact and HUPFD 1)HUPFD 2Absolute Error (Exact and HUPFD 2)HUPFD 3Absolute Error (Exact and HUPFD 3)
( 1.3334 , 0.0002 ) 0.0048 0.0048 3.0635 × 10 7 0.0049 5.3037 × 10 5 0.0048 5.1776 × 10 6 0.0049 2.9297 × 10 5
( 2.6668 , 0.0004 ) 0.0185 0.0186 9.2639 × 10 5 0.0186 8.1664 × 10 5 0.0184 2.2392 × 10 5 0.0188 3.9185 × 10 4
( 4.0002 , 0.0006 ) 0.0709 0.0715 5.9143 × 10 4 0.0714 5.1625 × 10 4 0.0700 1.4136 × 10 4 0.0722 2.4000 × 10 3
Table 3. Infinity norm results of the linear ADR equation by using the exact, UPFD, Crank–Nicolson, NSFD, and HUPFD solutions.
Table 3. Infinity norm results of the linear ADR equation by using the exact, UPFD, Crank–Nicolson, NSFD, and HUPFD solutions.
( x , t ) ExactUPFDErrorHUPFD 1ErrorHUPFD 2ErrorHUPFD 3ErrorCrank NicolsonErrorNSFDError
u ( : , 0.0002 ) 1.1654 1.1654 0.0025 1.1654 0.0025 1.1655 0.0025 1.1653 0.0025 1.1654 0.0024 1.1654 0.0025
u ( : , 0.0004 ) 1.1657 1.1657 0.0026 1.1657 0.0025 1.1657 0.0025 1.1652 0.0027 1.1657 0.0025 1.1657 0.0024
u ( : , 0.0006 ) 1.1659 1.1660 0.0026 1.1660 0.0026 1.1660 0.0028 1.1652 0.0030 1.1660 0.0025 1.1660 0.0026
Time 0.001230 0.002748 0.002717 0.003089 0.002689 0.001033
Table 4. Convergence rates for the nonlinear ADR reaction equation by the UPFD and HUPFD with respect to varying Δ x and Δ t .
Table 4. Convergence rates for the nonlinear ADR reaction equation by the UPFD and HUPFD with respect to varying Δ x and Δ t .
q ( Δ t k )
Δ t k Δ x k UPFDHUPFDCrank–Nicolson
0.00005 0.6667 1.0002 1.0029 1.9634
0.0001 0.6667 1.0004 1.0058 1.9055
0.0002 0.6667 1.0007 1.0116 2.2109
0.0004 0.6667
Table 5. Numerical results of the nonlinear ADR equation at a fixed Δ x = 0.01 , Δ t = 0.05 , r = 0.05 , and D = 0.0002 .
Table 5. Numerical results of the nonlinear ADR equation at a fixed Δ x = 0.01 , Δ t = 0.05 , r = 0.05 , and D = 0.0002 .
( x , t ) Exact SolutionUPFD SolutionHUPFD SolutionAbsolute Error (Exact and UPFD)Absolute Error (Exact and HUPFD)
( 0.3 , 0.06 ) 1.327 × 10 7 3.9120 × 10 8 1.314 × 10 7 9.346 × 10 8 1.201 × 10 9
( 0.4 , 0.08 ) 5.27 × 10 7 1.76 × 10 7 5.204 × 10 7 3.509 × 10 7 6.681 × 10 9
( 0.5 , 0.1 ) 1.964 × 10 6 7.322 × 10 7 1.933 × 10 6 1.232 × 10 6 3.195 × 10 8
Time 0.0060.086
Table 6. Infinity norm results of the nonlinear ADR equation by using the exact, UPFD, NSFD, and HUPFD solutions.
Table 6. Infinity norm results of the nonlinear ADR equation by using the exact, UPFD, NSFD, and HUPFD solutions.
( x , t ) ExactUPFDErrorHUPFDErrorNSFDError
u ( : , 0.1 ) 0.71010.69950.01060.70000.01010.69940.0107
u ( : , 0.2 ) 0.71010.69850.01160.70000.01160.69830.0118
u ( : , 0.4 ) 0.71010.69790.01220.70000.01220.69720.0129
Table 7. Convergence rates for the Schnakenberg model by the UPFD, HUPFD 1, HUPFD 2, and HUPFD 3 with respect to varying Δ x and Δ t .
Table 7. Convergence rates for the Schnakenberg model by the UPFD, HUPFD 1, HUPFD 2, and HUPFD 3 with respect to varying Δ x and Δ t .
q ( Δ t k )
Δ t k Δ x k UPFDHUPFD 1HUPFD 2HUPFD 3
0.00005 0.6667 1.0002 1.0029 1.0027 1.0027
0.0001 0.6667 1.0004 1.0058 1.0055 1.0055
0.0002 0.6667 1.0007 1.0116 1.0109 1.0109
0.0004 0.6667
Table 8. Comparison of solutions u of the Schnakenberg model by using the UPFD, HUPFD1, HUPFD2, and HUPFD3 solutions.
Table 8. Comparison of solutions u of the Schnakenberg model by using the UPFD, HUPFD1, HUPFD2, and HUPFD3 solutions.
( x , t ) UPFDErrorHUPFD 1ErrorHUPFD 2ErrorHUPFD 3Error
u ( : , 0.0002 ) 0.0048 0.0013 0.0050 0.0148 0.5254 0.0142 0.5254 0.0142
u ( : , 0.0004 ) 0.0182 0.0013 0.0198 0.0151 0.1456 0.0146 0.1456 0.0146
u ( : , 0.0006 ) 0.0690 0.0013 0.0788 0.0154 0.0403 0.0149 0.0403 0.0149
Time 0.048534 0.046653 0.057943 0.048868
Table 9. Comparison of solutions v of the Schnakenberg model by using the UPFD, HUPFD1, HUPFD2, and HUPFD3 solutions.
Table 9. Comparison of solutions v of the Schnakenberg model by using the UPFD, HUPFD1, HUPFD2, and HUPFD3 solutions.
( x , t ) UPFDErrorHUPFD 1ErrorHUPFD 2ErrorHUPFD 3Error
u ( : , 0.0002 ) 0.0048 0.0030 0.5310 0.0186 0.5026 0.0108 0.5271 0.0161
u ( : , 0.0004 ) 0.0186 0.0029 0.1429 0.0187 0.1350 0.0108 0.1460 0.0163
u ( : , 0.0006 ) 0.0711 0.0029 0.0390 0.0188 0.0355 0.0108 0.0405 0.0164
Time 0.048534 0.046653 0.057943 0.048868
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Ndou, N.; Dlamini, P.; Jacobs, B.A. Developing Higher-Order Unconditionally Positive Finite Difference Methods for the Advection Diffusion Reaction Equations. Axioms 2024, 13, 247. https://doi.org/10.3390/axioms13040247

AMA Style

Ndou N, Dlamini P, Jacobs BA. Developing Higher-Order Unconditionally Positive Finite Difference Methods for the Advection Diffusion Reaction Equations. Axioms. 2024; 13(4):247. https://doi.org/10.3390/axioms13040247

Chicago/Turabian Style

Ndou, Ndivhuwo, Phumlani Dlamini, and Byron Alexander Jacobs. 2024. "Developing Higher-Order Unconditionally Positive Finite Difference Methods for the Advection Diffusion Reaction Equations" Axioms 13, no. 4: 247. https://doi.org/10.3390/axioms13040247

APA Style

Ndou, N., Dlamini, P., & Jacobs, B. A. (2024). Developing Higher-Order Unconditionally Positive Finite Difference Methods for the Advection Diffusion Reaction Equations. Axioms, 13(4), 247. https://doi.org/10.3390/axioms13040247

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop