Next Article in Journal
Polynomial Intermediate Checksum for Integrity under Releasing Unverified Plaintext and Its Application to COPA
Previous Article in Journal
On Semi-Infinite Optimization Problems with Vanishing Constraints Involving Interval-Valued Functions
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method

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 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.
Mathematics 2024, 12(7), 1009; https://doi.org/10.3390/math12071009
Submission received: 20 February 2024 / Revised: 18 March 2024 / Accepted: 22 March 2024 / Published: 28 March 2024
(This article belongs to the Section Computational and Applied Mathematics)

Abstract

:
In this paper, the enhanced higher-order unconditionally positive finite difference method is developed to solve the linear, non-linear and system advection diffusion reaction equations. Investigation into the effectiveness and efficiency of the proposed method is carried out by calculating the convergence rate, error and computational time. A comparison of the solutions obtained by the enhanced higher-order unconditionally positive finite difference and exact solution is conducted for validation purposes. The numerical results show that the developed method reduced the time taken to solve the linear and non-linear advection diffusion reaction equations as compared to the results obtained by the higher-order unconditionally positive finite difference method.

1. Introduction

Numerical methods typically require large mesh sizes to give better approximations to solutions of differential equations. Consequently, their memory requirements are huge and, therefore, computational efficiency is affected. In this work, we aim to improve the computational efficiency of the higher-order unconditionally positive finite difference method (HUPFD). The HUPFD is based on using higher-order finite difference schemes to solve differential equations, with its main attribute being the preservation of the positivity of solutions. The positivity of solutions is essential in several physical settings—for example, quantities like population sizes, chemical species concentration and neutron numbers require positive solutions in order to be meaningful.
One of the techniques researchers employ to reduce the computational span of computational methods is the proper orthogonal decomposition method (POD) [1]. This has previously been used to reduce the dimensions of numerical methods such as non-standard finite difference and Crank–Nicolson to solve problems such as the Sobolev, hyperbolic, Burgers and Navier–Stokes equations [2,3,4]. In this paper, we couple the POD and HUPFD to obtain a new enhanced hybrid method that we name the enhanced higher-order unconditionally positive finite difference method (EHUPFD). By utilizing the POD, the EHUPFD extracts a set of basis functions from the HUPFD solution called the snapshot matrix and then uses a small subset of leading basis functions to construct state variable approximations.
We test the applicability of the EHUPFD on the advection diffusion reaction equation below
u t + U a u x D 2 u x 2 = f ( t , x , u ) , a x b t > 0 , x Ω 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 is the advection coefficient, u x ( t , x ) is the advection term, u x x ( t , x ) is the diffusion term, u t ( t , x ) represent the rate of change of concentration of substances and f ( t , x , u ) is the source or reaction term [5]. The ADR equation governs the process of diffusion and advection simultaneously [6]. The ADR equation is used to model the exponential traveling waves semiconductor modeling, biological system modeling, pollutant absorption in soil and neutron diffusion [6,7,8,9,10]. The UPFD method has previously been used to solve the ADR equation, amongst other known numerical methods reported in the literature, such as the standard finite difference, Crank–Nicolson finite difference, and non-standard finite difference methods and acts. The UPFD method has been proven to be a stable numerical method, accurate and consistent when solving the ADR equation in rectangular domains. The method is also known to preserve the positivity of solution regardless of step size [11]. For example, Appadu, A and Rao used UPFD to approximate the solution to the linear ADR equation, which models exponential traveling waves. The authors studied various scenarios in which advection or diffusion is dominating, as well as when both are equal. In this study, they discovered that the UPFD scheme is successful in solving the ADR equation when advection or diffusion factors are dominating [12]. However, as always, there is a need to constantly refine the numerical methods in order to improve the method’s accuracy, computational time and convergence rate.
The HUPFD method was developed to solve the advection diffusion reaction equations to improve the solution’s accuracy [13]. The higher-order numerical methods have been used to solve mathematical problems such as [13,14,15,16,17,18,19,20]. The standard finite difference methods, such as the backward, forward and central difference schemes, struggle to approach the exact solution of most partial differential models with large step sizes. As a result, relatively small step sizes are chosen to improve the accuracy of the solution to the problem at the cost of increased computational time due to the increased number of grid points. With a smaller grid size, HUPFD calculates a more accurate solution to the ADR equation’s linear, non-linear and system versions compared to the UPFD results.
The highly accurate finite difference method has gained much attention in recent years [21,22,23,24,25]. Gemechis FD and Tesfaye AB studied the higher-order numerical method to solve the singularly perturbed delay reaction–diffusion equation [26]. In order to achieve the six higher-order schemes, they used the Richardson extrapolation technique [26]. They found that the developed method is stable and consistent, and produces a more accurate solution than some of the methods reported in the literature. Ranjana KM and Venn G investigated the fourth-order finite difference method based on spline in tension approximation for the solution of second-order quasi-linear hyperbolic equations in one space dimension [27]. They proposed a three-level implicit nine-point compact finite difference formulation of order two in the time and four in the space directions based on split tension approximation in the t-direction for the numerical solution of space dimensional second-order quasi-linear hyperbolic equations with a first-order space derivative term. The proposed numerical method for the wave equation in polar coordinates is conditionally stable [27]. In contrast, the method for the damped wave and telegraphic equations is unconditionally stable. Sari M et al. [28] proposed a tenth-order finite difference scheme for solving a one-dimensional diffusion equation. Gurarslan et al. [29] presented a sixth-order compact difference scheme and a fourth-order Runge–Kutta scheme in time to produce numerical solutions of the one-dimensional advection–diffusion equation, which was found to be accurate in solving the concurrent transport equation for the one-dimensional advection diffusion reaction equation.
The rest of the sections in the paper are organized as follows. Section 2 and Section 3 introduce the UPFD and HUPFD, respectively. Subsequently, Section 4 and Section 5 apply the schemes in Section 2 and Section 3 to solve the ADR equations. Section 6 introduces the POD technique used to enhance the HUPFD, which results in the development of EHUPFD. The results obtained by the EHUPFD are illustrated in Section 6 and a summary of the results is provided in Section 7.

2. The Unconditionally Positive Finite Difference Scheme

The UPFD schemes are presented below [30],
u x { (2a) u i n u n 1 n + 1 Δ x (2b) u i n + 1 u n 1 n Δ x
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 grid points are calculated by using the formula x i = a + Δ x ( i 1 ) for 1 i m and spatial step size Δ x = b a m 1 . The time step is calculated by using the formula Δ t = T N 1 with grid points given by t n = Δ t ( n 1 ) , for 1 n N . Equation (2a) is used when the coefficient of the first space derivative is negative, whereas Equation (2b) is used when the first space derivative term is positive. Equation (4) is used when the second space derivative is negative or positive.

3. The Higher-Order Unconditionally Positive Finite Difference Schemes for Coupling with POD

In this section, we present the HUPFD that will be coupled with the POD to develop the EHUPFD. The EHUPFD is then tested on the ADR equations.

HUPFD

In this case, we consider the first-order formula (2a) and (2b) for the first space derivative, (3) for the first time derivative, and the fourth order for the second space derivative approximations to obtain the higher-order schemes as given below,
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 Δ x = b a m 1 and the grid points are given as x i = a + Δ x ( i 1 ) and 2 i m + 1 . The time step size Δ t = T N 1 , with the corresponding grid points given as t n = Δ t ( n 1 ) , 1 n N .
In the development of the higher-order unconditionally positive finite difference (HUPFD), the goal is to produce a higher-order method, based on the principles of the UPFD, that will still preserve the positivity of solutions. We test the effectiveness of the scheme on the ADR equations. In practice, selecting different orders in finite difference schemes for space derivatives involves a tradeoff between accuracy, stability, consistency and computational cost. Engineers and scientists choose different orders of the schemes based on the specific requirements of their simulations and the constraints of their computational resources, which is often about achieving accuracy. In this study, the HUPFD schemes are constructed by using different orders in finite difference schemes for the space derivative scheme with the idea of investigating a combination that gives better results.

4. Implementation of the HUPFD 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 following equation
u t + u x = 2 u x 2 u ,
The boundary and initial conditions used in this problem are given below
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 ) [30].

4.1. Solving the Linear ADR Equation Using the HUPFD Scheme

Because the coefficient of u x in Equation (6) is positive, we find the UPFD scheme using Equations (2b), (3) and (5). As a result, the following equation gives the HUPFD discretization of Equation (6)
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 (10) 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.2. Solving the Non-Linear ADR Equation Using the HUPFD Scheme

In this section, the EHUPFD is used to solve the non-linear advection diffusion reaction equation given below, with f ( t , x , u ) = r u ( 1 u ) and U a = 0
u t D 2 u x 2 = r u ( 1 u ) ,
The boundary and initial conditions of (15) are given below
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 and x 0 = 0.5 [30].
The exact solution is given as follows
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 and r = 0.05 .
The discretization of the non-linear ADR using Equations (2b), (3) and (15) 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 ,
where u m a x is a frozen coefficient. By simplifying Equation (18), 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 and s 3 = 1 2 Δ t ; 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 (19) 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 n = s 3 u 1 n u 2 n u 3 n · · · · · · u M 1 n u M n .

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

The linear Schnakenberg model is given by the equation below
u t + c u x = ( q 1 1 ) u + ( q 2 + s ) v + μ 1 2 u x 2 ,
v t + c v x = q 1 u q 2 v + μ 2 2 v x 2 ,
subject to the boundary and initial conditions given below
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 and μ 2 = 5 .

4.3.1. Solving the Linear Schnakenberg Model Using the UPFD Scheme

The UPFD scheme is applied to the linear Schnakenberg model and the discretization of the scheme is given 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 and s 2 = 1 Δ t + c Δ x + q 2 + 2 μ 2 ( Δ x ) 2 ; 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 ,
u 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 ) ,

4.3.2. Prove the Stability of the System ADR Equation with the UPFD

By discretizing (54) and (55) by using the UPFD scheme, we obtain the following
1 Δ t + ( 1 q 1 ) + 2 μ 1 ( Δ x ) 2 u i n + 1 + c Δ x u i + 1 n + 1 = q 2 v i n + μ 1 ( Δ x ) 2 u i + 1 n + u i 1 n + 1 Δ t u i n + c Δ x u i n ,
1 Δ t + q 2 + 2 μ 2 ( Δ x ) 2 v i n + 1 + c Δ x v i + 1 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 ,
By using the parameters s = 5 , c = 0 , q 1 = q 2 = 0 , μ 1 = 1 and μ 2 = 5 , Equations (34) and (35) are simplified to
1 Δ t + 1 + 2 ( Δ x ) 2 u i n + 1 = 1 ( Δ x ) 2 u i + 1 n + u i 1 n + 1 Δ t u i n
1 Δ t + 10 ( Δ x ) 2 v i n + 1 = 1 Δ t v i n + 5 ( Δ x ) 2 v i + 1 n + v i 1 n
The von Neumann stability analysis, u i n = ξ n e j Δ t i Δ x , is used to check the stability region in Equations (36) and (37) to obtain the following
1 Δ t + 1 + 2 ( Δ x ) 2 ξ n + 1 e j Δ t i Δ x = 1 ( Δ x ) 2 ξ n e j Δ t ( i + 1 ) Δ x + ξ n e j Δ t ( i 1 ) Δ x + 1 Δ t ξ n e j Δ t i Δ x ,
1 Δ t + 10 ( Δ x ) 2 ξ n + 1 e j Δ t i Δ x = 1 Δ t ξ n e j Δ t i Δ x + 5 ( Δ x ) 2 ξ n e j Δ t ( i + 1 ) Δ x + ξ n e j Δ t ( i 1 ) Δ x ,
Equations (38) and (39) are simplified to obtain the following
ξ = 2 μ 1 ( Δ x ) 2 cos ( ω ) + 1 Δ t 1 Δ t + 1 + 2 ( Δ x ) 2 ,
and
ξ = 1 Δ t + 10 cos ( ω ) ( Δ x ) 2 1 Δ t + 10 ( Δ x ) 2 ,
where ω = Δ t Δ x . By using the inequality | ξ | 1 , when ω = π , we obtain the following, respectively
Δ t > 2 .

4.3.3. Proof of Consistency with the System ADR Equation by Using UPFD Scheme

By discretizing (24) and (25) by using the UPFD scheme, we obtain the following
u i n + 1 u i n Δ t + c u i + 1 n + 1 u i 2 Δ x = ( q 1 1 ) u i n + 1 + ( q 2 + s ) v i n + μ 1 ( Δ x ) 2 u i + 1 n 2 u i n + 1 + u i 1 n ,
v i n + 1 v i n Δ t + c v i + 1 n + 1 v i 2 Δ x = q 1 u i n + 1 q 2 v i n + μ 2 ( Δ x ) 2 v i + 1 n 2 v i n + 1 + v i 1 n .
From (43) and (44), we obtain the following Taylor’s expansions
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 applying Taylor’s expansions to Equations (43) and (44), we obtain the following simplified equations
( 1 q 1 ) u i n + 1 + ( 1 q 1 ) Δ t 2 μ 1 Δ x 2 Δ t + c Δ t Δ x u t + c u x + Δ t 2 + ( 1 q 1 ) ( Δ t ) 2 2 + 2 μ 1 ( Δ x ) 2 + 2 μ 1 ( Δ t ) 3 6 ( Δ x ) 2 + c ( Δ t ) Δ x 2 u t 2 + 2 c Δ t 2 u t x q 2 v i n μ 1 2 u x 2 = 0 ,
and
1 Δ t + q 2 + 2 μ 2 ( Δ x ) 2 + c Δ x v i n + 1 + q 2 Δ t + 2 μ 2 Δ t ( Δ x ) 2 + c Δ t Δ x v t + c v x + Δ t 2 + q 2 ( Δ t ) 2 2 + μ 2 ( Δ t ) 2 ( Δ x ) 2 + c ( Δ t ) 2 2 ( Δ x ) 2 v t 2 + 2 Δ t 2 v t x + c ( Δ x ) 2 2 v x 2 = 1 1 u i n + 1 Δ t + 2 μ 2 ( Δ x ) 2 v i n q 1 Δ t u t q 1 ( Δ t ) 2 2 2 u t 2 q 1 ( Δ t ) 3 6 3 u t 3 + μ 2 2 v x 2 c Δ x v i n ,
Let Δ t = ( Δ x ) 3 , which leads to Equations (51) and (52) being written as follows
( 1 q 1 ) u i n + 1 + ( 1 q 1 ) ( Δ x ) 3 + 2 μ 1 Δ x + c ( Δ x ) 2 u x + c u x + ( Δ x ) 3 2 + ( 1 q 1 ) ( Δ x ) 6 2 + μ 1 ( Δ x ) 4 3 + c ( Δ x ) 2 2 u t 2 + 2 c ( Δ x ) 3 2 u t x q 2 v i n μ 1 2 u x 2 = 0 ,
and
q 2 v i n + q 1 u i n + 1 + q 2 ( Δ x ) 3 + 2 μ 2 Δ x + c ( Δ x ) 2 + q 1 ( Δ x ) 3 v t + c v x     + ( Δ x ) 3 2 + q 2 ( Δ x ) 6 2 + μ 2 ( Δ x ) 4 + d ( Δ x ) 5 2 + q 1 ( Δ x ) 6 2 2 u t 2    + 2 Δ t 2 v t x + c ( Δ x ) 2 μ 2 = 0 ,
When Δ x 0 , we obtain the original system ADR equations
u t + c u x = ( q 1 1 ) u + ( q 2 + s ) v + μ 1 2 u x 2
v t + c v x = q 1 u q 2 v + μ 2 2 v x 2 ,

4.3.4. Solving the Schnakenberg Model Using the HUPFD Scheme

The HUPFD scheme is applied to the linear Schnakenberg model and the discretization of the scheme is given 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 and p 3 = c Δ x + 16 μ 2 12 ( Δ x ) 2
Equations (58) and (59) 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 , respectively
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 v 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 v 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 v 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.

4.3.5. Prove Stability When Solving System ADR Equation by the HUPFD

By discretizing (56) and (57) by using the UPFD scheme, we obtain the following
u i n + 1 u i n Δ t + c u i + 1 n + 1 u i 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 n Δ x = q 1 u i n + 1 q 2 v i n + 1 + μ 1 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 ,
By using parameters s = 5 , c = 0 , q 1 = q 2 = 0 , μ 1 = 1 and μ 2 = 5 , Equations (67) and (70) are simplified to
1 Δ t + 30 12 ( Δ x ) 2 v i n + 1 + 5 12 ( Δ x ) 2 v i + 2 n + 1 + v i 1 n + 1 = 1 Δ t v i n + 80 12 ( Δ x ) 2 v i + 1 n + v i 1 n ,
v i n + 1 v i n Δ t + c v i + 1 n + 1 v i n Δ x = q 1 u i n + 1 q 2 v i n + 1        + μ 1 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 ,
The von Neumann stability analysis, u i n = ξ n e j Δ t i Δ x , is used to check the stability region in Equations (67) and (70) to obtain the following
1 Δ t + 1 q 1 + 30 μ 1 12 ( Δ x ) 2 ξ n + 1 e j Δ t i Δ x + 1 12 ( Δ x ) 2 ξ n + 1 e j Δ t ( i + 1 ) Δ x      + 1 12 ( Δ x ) 2 ξ n + 1 e j Δ t ( i 1 ) Δ x = 16 12 ( Δ x ) 2 ξ n e j Δ t ( i 1 ) Δ x     + 1 Δ t ξ n e j Δ t i Δ x + 16 12 ( Δ x ) 2 ξ n e j Δ t ( i + 1 ) Δ x ,
1 Δ t + 30 12 ( Δ x ) 2 ξ n + 1 e j Δ t i Δ x + 5 12 ( Δ x ) 2 ξ n + 1 e j Δ t ( i + 1 ) Δ x + ξ n + 1 e j Δ t ( i 2 ) Δ x        = 1 Δ t ξ n e j Δ t i Δ x + 80 12 ( Δ x ) 2 ξ n e j Δ t ( i + 1 ) Δ x + ξ n e j Δ t ( i 1 ) Δ x ,
Equations (71) and (72) are simplified to obtain the following
ξ = 16 6 ( Δ x ) 2 cos ( ω ) + 1 Δ t 1 Δ t + 1 + 30 12 ( Δ x ) 2 + 1 6 ( Δ x ) 2 cos ( ω ) ,
and
ξ = 1 Δ t + 80 cos ( ω ) 6 ( Δ x ) 2 1 Δ t + 30 ( Δ x ) 2 + 5 6 ( Δ x ) 2 cos ( ω ) ,
where ω = Δ t Δ x . By using the inequality | ξ | 1 , when ω = π , we obtain the following, respectively
Δ t > 24 4 12 ( Δ x ) 2 ,
and
Δ t > 6 ( Δ x ) 2 35 ,

4.3.6. Consistency with the System of ADR by Using HUPFD

By discretizing (56) and (57) by using the HUPFD1 scheme, we obtain the following
v i n + 1 v i n Δ t + c v i + 1 n + 1 v i n Δ x = q 1 u i n + 1 q 2 v i n + 1 + μ 1 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 ,
From (77), we obtain the following Taylor’s expansions
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 +
and
v i n + 1 v i n + Δ t v t + ( Δ t ) 2 2 2 v t 2 + ( Δ t ) 3 6 3 v t 3 +
v i n 1 v i n Δ t v t + ( Δ t ) 2 2 2 v t 2 ( Δ t ) 3 6 3 v t 3 +
v i 1 n v i n Δ x v x + ( Δ x ) 2 2 2 v x 2 ( Δ x ) 3 6 3 v x 3 +
v i + 1 n v i n + Δ x v x + ( Δ x ) 2 2 2 v x 2 + ( Δ x ) 3 6 3 v x 3 +
v i 2 n + 1 v i n + Δ t v t 2 Δ x v x + 2 ( Δ x ) 2 2 v x 2 2 ( Δ t Δ x ) 2 v t x + ( Δ t ) 2 2 2 v t 2 +
v i + 2 n + 1 v i n + Δ t v t + 2 Δ x v x + 2 ( Δ x ) 2 2 v x 2 + 2 ( Δ t Δ x ) 2 v t x + ( Δ t ) 2 2 2 v t 2 +
By applying Taylor’s expansions to (77), respectively, we obtain the following simplified equations
( 1 q 1 ) u i n + 1 + Δ t q 1 Δ t + 32 μ 2 Δ t 12 ( Δ x ) 2 + c Δ t 12 u t + c u x + Δ t 2 + ( Δ t ) 2 2 q 1 ( Δ t ) 2 2 + 17 μ 1 ( Δ t ) 2 12 ( Δ x ) 2 + c ( Δ t ) 2 2 Δ x 2 u t 2 + c ( Δ x ) 2 μ 1 2 u x 2 + c Δ t 2 2 u t x q 2 v i n = 0 ,
and
q 2 v i n + 1 + q 2 Δ t + 32 Δ t μ 2 12 ( Δ x ) 2 + c ( Δ t ) 2 Δ x v t + c v x + Δ t + q 2 ( Δ t ) 2 + 32 μ 2 ( Δ t ) 2 12 ( Δ ) 2 + c ( Δ ) 2 2 ( Δ x ) 2 v t 2 + c ( Δ x ) 2 4 μ 2 12 12 μ 2 12 2 v x 2    + q 1 u i n + Δ t u t + ( Δ t ) 2 2 2 u t 2 = 0 ,
Let Δ t = ( Δ x ) 3 , which leads to Equations (90) and (91) being written as follows
( 1 q 1 ) u i n + 1 + ( Δ x ) 3 q 1 ( Δ x ) 3 + 32 μ 1 Δ x 12 + c ( Δ x ) 2 u t + c u x + ( Δ x ) 3 2 + ( Δ x ) 6 2 q 1 ( Δ x ) 6 2 + 17 μ 1 ( Δ x ) 4 12 + c ( Δ x ) 5 2 2 u t 2 + c ( Δ x ) 2 μ 1 2 u x 2 + c ( Δ x ) 3 2 2 u t x q 2 v i n = 0 ,
and
q 2 v i n + 1 + q 2 ( Δ x ) 3 + 32 μ 2 Δ x 12 + c ( Δ x ) 5 v t + c v x + ( Δ x ) 3 + q 2 ( Δ x ) 6 + 32 μ 2 ( Δ x ) 4 12 + c ( Δ x ) 5 2 2 v t 2 + c ( Δ x ) 2 μ 2 2 v x 2 + q 1 u i n + ( Δ x ) 3 u t + ( Δ x ) 6 2 2 u t 2 = 0 ,
When Δ x 0 , we obtain the original system ADR equations
u t + c u x = ( q 1 1 ) u + ( q 2 + s ) v + μ 1 2 u x 2 ,
v t + c v x = q 1 u q 2 v + μ 2 2 v x 2 ,
Theorem 1
(Solution positivity). Let u 0 ( x ) 0 , 0 f 0 ( x ) 1 and v 0 ( x ) 0 for the numerical scheme (11), (20) and (60); we have u i n + 1 0 , 0 f i n + 1 and v i n + 1 0 , independent of the mesh step size [31].
Proof. 
Let us assume that u i n + 1 0 , 0 f i n + 1 1 and v i n 0 , for i = 0 , 1 , , N , n = 0 , N .
Equations (11), (20) and (60) comprise a coefficient matrix with a strictly diagonally dominant tridiagonal matrix, with positive values on the diagonal and sub-diagonals, which implies Equations (11), (20) and (60) yield non-negative inverses. Thus, the solution to the algebraic system is non-negative in the form of u i n + 1 0 and v i n + 1 0 , for i = 0 , 1 , N .
Furthermore, it is observed that f i n + 1 0 in equations (6), (15), (54) and (55). By applying the same consideration at each time level, we complete the proof. □

5. Formulation of the POD Technique

In this section, the POD basis function is constructed and used to reduce the degree of freedom to the UPFD and HUPFD scheme.
The POD technique is used for the HUPFD by extracting the snapshot matrix in order to produce the EHUPFD, which leads to the construction of the approximate state variable ψ i by developing a small subset of the leading basis function obtained from the SVD factorization. The POD technique seeks a low-dimensional basis function that can accurately approximate the state variable, which is a left singular vector from the SVD.
The first L sequence of solutions is extracted from the HUPFD solutions as given by u i n , v i n and s i n , which gives the formulated n × L snapshot matrix given below
A = u 1 1 u 1 2 · · · · · · u 1 L u 2 1 u 2 2 · · · · · · u 2 L · · · · · · · · · · · · · · · · · · · · · · · · · · · u m 1 u m 2 · · · · · · u m L ,
where m and L represent the size of the snapshot matrix, which is factorized by using the singular value decomposition (SVD) as given below
A = U Σ V T = U Σ l × l O l × ( L l ) O ( m 1 ) × l O ( m 1 ) × ( L l ) V T ,
where U = ( ψ 1 , ψ 2 , , ψ m ) is an m × m matrix representing an orthogonal matrix consisting of the eigenvectors of A A T , Σ l × l = diag ( σ 1 , σ 2 , , σ l ) represents the diagonal matrix from the SVD, given in decreasing order σ 1 σ 2 σ l > 0 , and V = ( ϕ 1 , ϕ 2 , , ϕ L ) is an L × L orthogonal matrix which represent the orthogonal eigenvector of A T A and l = r a n k A ; O is a zero matrix [32].
From identical eigenvalues of A T A and A A T , we obtain the eigenvalues of λ 1 λ 2 λ l > 0 ( l = r a n k A ) for matrices A A T and corresponding eigenvectors ψ j , which is given by the formula below
ψ j = A ϕ j / λ j , j = 1 , 2 , , l ,
This leads to the construction of the matrix A M given below
A M = U Σ V T = U Σ M × M O M × ( L M ) O ( m M ) × l O ( m M ) × ( L M ) V T ,
where the diagonal matrix given by Σ M × M = d i a g ( σ 1 , σ 1 , , σ m ) is constituted of the first M positive singular values of the diagonal matrix.

5.1. The Implementation of the Algorithm of the EHUPFD Scheme for the Advection Diffusion Reaction Equations

Below are the steps to implement the algorithm of the EHUPFD to solve the linear, non-linear and system ADR equations [32].
Step 1:
Snapshot matrix for HUPFD is extracted.
Step 2:
Formulation of the matrix A = ( u i n ) m × L and calculation of eigenvalues and eigenvectors of matrix A T A and A A T executed.
Step 3:
Choice of POD basis for the error tolerance of μ = O ( Δ t , Δ x 2 ) is executed.
Step 4:
The EHUPFD computes the reduced-order solution.
Step 5:
Accuracy of the solution is checked for the renewal of the POD basis.

5.2. The EHUPFD Error Estimate

The following equation is used to estimate the error of the EHUPFD when solving the ADR equations.
u n u * n 2 = u n u M n 2 = ( A ϕ ϕ T A ) ξ n 2 λ M + 1 ,
n = 1 , 2 , , L
Theorem 2
(Accuracy). The error estimate for the HUPFD u i n and u i * n from the EHUPFD is calculated by using the following equation
| u i n u i * n | ( 1 + δ ) n L λ M + 1 for ( L + 1 n N ) ,
where δ = Δ t Δ x 2 B 2 , 2 .
The absolute error, calculated between the actual and approximate solutions from the POD reduced-order model satisfy the following equation [32]
| u ( x i , t n ) u i * n | = O ( 1 + δ ) n L λ M + 1 , Δ t , Δ x 2 , 1 i M , 1 n N ,

6. Results and Discussion

In this section, we present solutions to the linear, non-linear and system ADR equations solved by using the EHUPFD, UPFD, NSFD and Crank–Nicolson methods. To check the performance of the methods, absolute errors, convergence rates and computational time are measured. The actual and numerical solutions are denoted by u i , j and u ¯ i , j , respectively. The error is represented by e i , j and its magnitude is calculated by using norm to infinity denoted by L -norm given by the formula below,
| | e i , n | | = | | u i , n u ¯ i , n | | = M a x | u i , n u ¯ i , n | ,
The formula to calculate the convergence rate is given below,
q ( Δ t k ) = l o g 2 e k e k + 1 l o g 2 Δ t Δ t k + 1 or q ( Δ x k ) = l o g 2 e k e k + 1 l o g 2 Δ x Δ x k + 1 ,
Finding the best snapshot size is paramount because it influences the EHUPFD solution accuracy. The best snapshot is determined by finding the solution to the problem when L = 1 and then an increment of 1 to L is evaluated until there is no significant change in the solution. The tolerance is calculated using the formula below
| | u e x a c t u L | | ϵ , L = 1 , 2 , 3 ,
where u L is a solution for a specific value of L and u e x a c t is the exact solution. The first value of L that satisfies Equation (106) is the optimal value of L. We set ϵ = 10 15 .

6.1. Results Obtained to Solve the Linear ADR Equation

The convergence rates of the UPFD, Crank–Nicolson, NSFD and HUPFD with respect to the time and spatial steps are shown in Table 1. The results show that the EHUPFD converges to the solution faster than the UPFD. Table 2 compares solutions obtained by the exact, UPFD, Crank–Nicolson, NSFD and EHUPFD solutions at ( t , x ) = 1.33334 , 0.0002 , 2.6668 , 0.0004 and 4.0002 , 0.0006 . Table 3 compares the infinity norm results of the developed EHUPFD method against the Crank Nicolson and NSFD, as well as their computational times. The results show that the EHUPFD takes less computational time to converge to the solution as compared to the HUPFD methods.
Table 3 compares the infinity norms of the proposed EHUPFD methods with the UPFD, Crank–Nicolson and NSFD methods in solving the linear ADR equation. The results show that EHUPFD methods preserve solution accuracy.
Figure 1 shows the contribution made by each POD basis mode. The first five modes contain the predominant information where the EHUPFD meets the theoretical accuracy requirement. We observed that five modes are significant in optimizing the linear ADR equation. Figure 2b compares the number of POD modes against the absolute error. We observe that only five POD modes are required in order to converge to the toleration of 10 14 .
Figure 3a to Figure 4a show the surface plots of the exact, UPFD and EUPFD solutions of the linear ADR. Figure 4b shows the increase in error when the POD modes increases. The figures show that the developed methods achieve positivity of the solutions and the results are comparable with well-known methods such as Crank–Nicolson and NSFD.

6.2. Results Obtained by Using the Non-Linear ADR Equation

The convergence rates of the UPFD, Crank–Nicolson, NSFD and the EHUPFD with respect to the time and spatial steps are shown in Table 4. Table 5 compares the exact, UPFD and EHUPFD solutions at ( t , x ) = 0.03338 , 0.002 , 0.0676 , 0.004 and 0.1014 , 0.006 . The HUPFD has a slightly higher convergence rate than UPFD, while the EHUPFD has significantly reduced the computational time taken to converge to the solution of the ADR equation.
In Table 6, we compare the infinity norms of the proposed HUPFD with the UPFD, Crank–Nicolson and NSFD methods in solving the linear ADR equation. For this example, the results show that the EHUPFD method preserves the accuracy and reduces the computational time taken to converge to the solution of the ADR equation.
Figure 5a shows the contribution made by each POD basis mode. The first five modes contain the predominant information where the EHUPFD meets the theoretical accuracy requirement. We observed five modes are significant in optimizing the linear ADR equation. Figure 5b shows the increase in computational time when the number of modes increases.
Figure 6a compares the results obtained by the EHUPFD, Crank–Nicolson and NSFD methods. The results show that the EHUPFD solution is comparable with the Crank–Nicolson and NSFD solutions.
Figure 6b compares the error against the number of modes. The results show that the error of the solution decreases when the number of modes increases.
Figure 7a to Figure 8a show the surface plots of the exact, UPFD and EUPFD solutions of the non-linear ADR equation. Figure 8b show the absolute errors between the exact solution and the EHUPFD solution for the non-linear ADR against time.

6.3. Results Obtained When Solving the Schnakenberg Model

The convergence rates of the UPFD, HUPFD and EHUPFD with respect to the time and spatial steps are shown in Table 7. The results show that the convergence rate of the EHUPFD methods is slightly higher than the UPFD when solving the Schnakenberg model. Table 8 and Table 9 compare the solutions of the UPFD, HUPFD and EHUPFD solutions at ( t , x ) = 1.33334 , 0.0002 , 2.6668 , 0.0004 and 4.0002 , 0.0006 , as well as the computational times. The results show that the EHUPFD takes less computational time to converge to the solution than the UPFD, HUPFD and EHUPFD methods. Figure 9a to Figure 10b show the surface plots of the UPFD, HUPFD and EHUPFD of the Schnakenberg model. The figures show that the developed methods achieve positivity of the solutions. Figure 11 compares the solutions obtained by the exact HUPFD, EHUPFD and UPFD. The results show that the proposed methods are comparable to the known methods of the UPFD.
The convergence rates of the UPFD, HUPFD and EHUPFD with respect to the time and spatial steps are shown in Table 7. The results show that the convergence rate of the EHUPFD methods is slightly higher than the UPFD when solving the Schnakenberg model. Table 8 and Table 9 compare the UPFD, HUPFD and EHUPFD solutions at ( x , t ) = 1.33334 , 0.0002 , 2.6668 , 0.0004 and 4.0002 , 0.0006 , as well as the computational times. The results show that the EHUPFD takes less computational time to converge to the solution than the UPFD, HUPFD and EHUPFD methods. Figure 9a to Figure 10b show the surface plots of the UPFD, HUPFD and EHUPFD of the Schnakenberg model. The figures show that the developed methods achieve positivity of the solutions. Figure 11 compares the solutions obtained by the exact HUPFD, EHUPFD and UPFD. The results show that the proposed methods are comparable to the know methods of the UPFD.

7. Conclusions

This paper uses the EHUPFD to find solutions to the linear, non-linear and system advection diffusion reaction equations. The goal was to investigate the performance of the EHUPFD in terms of computational time, absolute error and convergence rate. The results show that the EHUPFD gives solutions that are a good approximation to solve the linear, non-linear and system ADR equations, which validates the method’s efficiency because it confirms the theoretical results. When solving the ADR equations by EHUPFD, a few modes containing the predominant information are needed to satisfy the theoretical results. The theoretical results are achieved using only 4 POD bases from the HUPFD’s 40 classical solutions as a snapshot. This implies that the EHUPFD saves computational time, reduces the degree of freedom and alleviates truncation error. The EHUPFD method also preserves the positivity of the solutions.

Author Contributions

Formal analysis, N.N.; Investigation, N.N.; Supervision, P.D. and B.A.J.; Writing—original draft, N.N.; Writing—review and 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 Venda.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data was created in this paper.

Conflicts of Interest

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

Abbreviations

EHUPFDEnhanced higher-order unconditionally positive finite difference methods
ADRAdvection diffusion reaction
PODProper orthogonal decomposition
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. Liang, Y.; Lee, H.; Lim, S.; Lin, W.; Lee, K.; Wu, C. Proper orthogonal decomposition and its applications—Part I: Theory. J. Sound Vib. 2002, 252, 527–544. [Google Scholar] [CrossRef]
  2. Luo, Z.; Teng, F.; Xia, H. A reduced-order extrapolated Crank–Nicolson finite spectral element method based on POD for the 2D non-stationary Boussinesq equations. J. Math. Anal. Appl. 2019, 471, 564–583. [Google Scholar] [CrossRef]
  3. Luo, Z.; Li, H.; Sun, P. A reduced-order Crank–Nicolson finite volume element formulation based on POD method for parabolic equations. Appl. Math. Comput. 2013, 219, 5887–5900. [Google Scholar] [CrossRef]
  4. Zhou, Y.; Luo, Z. An optimized Crank–Nicolson finite difference extrapolating model for the fractional-order parabolic-type sine-Gordon equation. Adv. Differ. Equ. 2019, 2019, 1. [Google Scholar] [CrossRef]
  5. Britton, N. Others Reaction-Diffusion Equations and Their Applications to Biology; Academic Press: New York, NY, USA, 1986. [Google Scholar]
  6. 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]
  7. Ndou, N.; Dlamini, P.; Jacobs, B. Enhanced Unconditionally Positive Finite Difference Method for Advection–Diffusion–Reaction Equations. Mathematics 2022, 10, 2639. [Google Scholar] [CrossRef]
  8. Liu, B.; Allen, M.; Kojouharov, H.; Chen, B. Others Finite-element solution of reaction-diffusion equations with advection. Comput. Methods Water Resour. 1996, 1, 3–12. [Google Scholar]
  9. Khan, H.; Gómez-Aguilar, J.; Khan, A.; Khan, T. Stability analysis for fractional order advection–reaction diffusion system. Phys. Stat. Mech. Its Appl. 2019, 521, 737–751. [Google Scholar] [CrossRef]
  10. 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]
  11. 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]
  12. Appadu, A. Performance of UPFD scheme under some different regimes of advection, diffusion and reaction. Int. J. Numer. Methods Heat Fluid Flow 2017, 27, 1412–1429. [Google Scholar] [CrossRef]
  13. 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]
  14. Schwendt, M.; Pötz, W. Transparent boundary conditions for higher-order finite-difference schemes of the Schrödinger equation in (1+ 1) D. Comput. Phys. Commun. 2020, 250, 107048. [Google Scholar] [CrossRef]
  15. 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]
  16. Georgakopoulos, S.; Birtcher, C.; Balanis, C.; Renaut, R. Higher-order finite-difference schemes for electromagnetic radiation, scattering, and penetration. 1. Theory. IEEE Antennas Propag. Mag. 2002, 44, 134–142. [Google Scholar] [CrossRef]
  17. Bilbao, S.; Hamilton, B. Higher-order accurate two-step finite difference schemes for the many-dimensional wave equation. J. Comput. Phys. 2018, 367, 134–165. [Google Scholar] [CrossRef]
  18. Chung, Y.; Tucker, P. Accuracy of higher-order finite difference schemes on nonuniform grids. AIAA J. 2003, 41, 1609–1611. [Google Scholar] [CrossRef]
  19. Morinishi, Y.; Lund, T.; Vasilyev, O.; Moin, P. Fully conservative higher order finite difference schemes for incompressible flow. J. Comput. Phys. 1998, 143, 90–124. [Google Scholar] [CrossRef]
  20. 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]
  21. Duressa, G.; Bullo, T.; Kiltu, G. Fourth order compact finite difference method for solving one dimensional wave equation. Int. J. Eng. Appl. Sci. 2016, 8, 30–39. [Google Scholar] [CrossRef]
  22. Cui, M. Compact finite difference method for the fractional diffusion equation. J. Comput. Phys. 2009, 228, 7792–7804. [Google Scholar] [CrossRef]
  23. 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]
  24. 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]
  25. 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]
  26. Duressa, G.; Bullo, T. Higher-Order Numerical Method for Singularly Perturbed Delay Reaction-Diffusion Problems. Pure Appl. Math. 2021, 10, 68–76. [Google Scholar] [CrossRef]
  27. Mohanty, R.; Gopal, V. A fourth-order finite difference method based on spline in tension approximation for the solution of one-space dimensional second-order quasi-linear hyperbolic equations. Adv. Differ. Equ. 2013, 2013, 70. [Google Scholar] [CrossRef]
  28. 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]
  29. 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]
  30. 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]
  31. Kolev, M.; Koleva, M.; Vulkov, L. An Unconditional Positivity-Preserving Difference Scheme for Models of Cancer Migration and Invasion. Mathematics 2022, 10, 131. [Google Scholar] [CrossRef]
  32. Luo, Z.; Chen, G. Proper Orthogonal Decomposition Methods for Partial Differential Equations; Academic Press: New York, NY, USA, 2018. [Google Scholar]
Figure 1. (a) Plot for the number of the POD modes against the eigenvalues with significant POD modes. (b) Plot for the number of POD modes against computational time when solving the linear ADR equation.
Figure 1. (a) Plot for the number of the POD modes against the eigenvalues with significant POD modes. (b) Plot for the number of POD modes against computational time when solving the linear ADR equation.
Mathematics 12 01009 g001
Figure 2. (a) Plots for the UPFD, EHUPFD, Crank–Nicolson and NSFD solutions when solving the linear ADR equation. (b) shows the plot for the number of POD modes against the absolute error.
Figure 2. (a) Plots for the UPFD, EHUPFD, Crank–Nicolson and NSFD solutions when solving the linear ADR equation. (b) shows the plot for the number of POD modes against the absolute error.
Mathematics 12 01009 g002
Figure 3. (a) Surface plot of the UPFD solution when solving the linear ADR equation. (b) shows the surface plot for the EHUPFD solution when solving the linear ADR equation.
Figure 3. (a) Surface plot of the UPFD solution when solving the linear ADR equation. (b) shows the surface plot for the EHUPFD solution when solving the linear ADR equation.
Mathematics 12 01009 g003
Figure 4. (a) Surface plot of the exact solution when solving the linear ADR equation. (b) Plot for the time against the absolute difference.
Figure 4. (a) Surface plot of the exact solution when solving the linear ADR equation. (b) Plot for the time against the absolute difference.
Mathematics 12 01009 g004
Figure 5. (a) Plot of the number of POD modes against the eigenvalues with significant POD modes. (b) Number of the POD modes against the computational time.
Figure 5. (a) Plot of the number of POD modes against the eigenvalues with significant POD modes. (b) Number of the POD modes against the computational time.
Mathematics 12 01009 g005
Figure 6. (a) Plot for the solutions obtained by the exact, HUPFD, Crank–Nicolson and NSFD methods when solving the non-linear ADR equation. (b) Plot of the number of POD modes against the absolute error.
Figure 6. (a) Plot for the solutions obtained by the exact, HUPFD, Crank–Nicolson and NSFD methods when solving the non-linear ADR equation. (b) Plot of the number of POD modes against the absolute error.
Mathematics 12 01009 g006
Figure 7. (a) Surface plot of the exact solution when solving the the non-linear ADR equation. (b) Surface plot of the UPFD solutions when solving the non-linear ADR equation.
Figure 7. (a) Surface plot of the exact solution when solving the the non-linear ADR equation. (b) Surface plot of the UPFD solutions when solving the non-linear ADR equation.
Mathematics 12 01009 g007
Figure 8. (a) Surface plot of the EHUPFD solution when solving the non-linear ADR equation. (b) Plot of the time against the absolute error.
Figure 8. (a) Surface plot of the EHUPFD solution when solving the non-linear ADR equation. (b) Plot of the time against the absolute error.
Mathematics 12 01009 g008
Figure 9. (a) Exact solutions obtained when solving the Schnakenberg system equations. (b) UPFD solutions obtained when solving the Schnakenberg system equations.
Figure 9. (a) Exact solutions obtained when solving the Schnakenberg system equations. (b) UPFD solutions obtained when solving the Schnakenberg system equations.
Mathematics 12 01009 g009
Figure 10. (a) HUPFD solutions obtained when solving the Schnakenberg system equations. (b) EHUPFD solutions when solving the Schnakenberg system equations.
Figure 10. (a) HUPFD solutions obtained when solving the Schnakenberg system equations. (b) EHUPFD solutions when solving the Schnakenberg system equations.
Mathematics 12 01009 g010
Figure 11. Comparison of plots for the solutions obtained by the UPFD, HUPFD and EHUPFD methods when solving the Schnakenberg system equations.
Figure 11. Comparison of plots for the solutions obtained by the UPFD, HUPFD and EHUPFD methods when solving the Schnakenberg system equations.
Mathematics 12 01009 g011
Table 1. Convergence rate results obtained from the UPFD, NSFD, Crank–Nicolson and HUPFD methods when solving the ADR equation.
Table 1. Convergence rate results obtained from the UPFD, NSFD, Crank–Nicolson and HUPFD methods when solving the ADR equation.
qtk)
ΔtkΔxk UPFD HUPFD EHUPFD Crank–Nicolson NSFD
0.00005 0.6667 2.1675 2.5233 2.5233 4.5606 2.0579
0.0001 0.6667 0.9973 1.0054 1.0054 1.1996 1.9427
0.0002 0.6667 1.8201 1.9990 1.9990 2.7372 1.8871
0.0004 0.6667
Table 2. Numerical results obtained when solving the linear ADR equation by using the UPFD, HUPFD, EHUPFD, NSFD and Crank–Nicolson methods.
Table 2. Numerical results obtained when solving the linear ADR equation by using the UPFD, HUPFD, EHUPFD, NSFD and Crank–Nicolson methods.
(x, t) Exact SolutionUPFD SolutionAbsolute ErrorEHUPFDAbsolute ErrorCrank–NicolsonAbsolute ErrorNSFDAbsolute Error
( 1.3334 , 0.0002 ) 0.5136 0.5137 9.3425 × 10 5 0.5125 0.0011 0.5138 1.6145 × 10 4 0.5137 9.4142 × 10 5
( 2.6668 , 0.0004 ) 0.1355 0.1356 8.1664 × 10 5 0.1366 0.0011 0.1356 6.7682 × 10 5 0.1356 7.4518 × 10 5
( 4.0002 , 0.0006 ) 0.0357 0.0358 5.1090 × 10 4 0.0363 3.2521 × 10 5 0.0358 2.9601 × 10 5 0.0358 3.2769 × 10 5
Table 3. Infinity norm results obtained by using the exact, UPFD, HUPFD, EHUPFD, Crank–Nicolson and NSFD methods to solve the linear ADR equation.
Table 3. Infinity norm results obtained by using the exact, UPFD, HUPFD, EHUPFD, Crank–Nicolson and NSFD methods to solve the linear ADR equation.
(x, t) ExactUPFDErrorHUPFDErrorEHUPFDErrorCrank–NicolsonErrorNSFDError
u ( : , 0.0002 ) 1.1654 1.1654 0.0025 1.1654 0.0025 1.1655 0.0025 1.1654 0.0025 1.1654 0.0025
u ( : , 0.0004 ) 1.1657 1.1657 0.0025 1.1657 0.0025 1.1657 0.0025 1.1657 0.0025 1.1657 0.0025
u ( : , 0.0006 ) 1.1659 1.1660 0.0025 1.1660 0.0025 1.1660 0.0025 1.1660 0.0025 1.1657 0.0025
Time 0.001230 0.002748 0.002717 0.003089 0.002689
Table 4. Convergence rate results obtained when solving the non-linear ADR reaction equation by the UPFD, HUPFD, Crank–Nicolson, NSFD and EHUPFD methods.
Table 4. Convergence rate results obtained when solving the non-linear ADR reaction equation by the UPFD, HUPFD, Crank–Nicolson, NSFD and EHUPFD methods.
qtk)
Δ t k Δ x k UPFDHUPFDEHUPFDCrank–NicolsonNSFD
0.0005 0.6667 1.7166 1.7197 1.7197 1.7197 1.71652
0.001 0.6667 1.7160 1.7150 1.7151 1.7150 1.71594
0.002 0.6667 1.7168 1.7195 1.7195 1.7195 1.71671
0.004 0.6667
Table 5. Solutions obtained by the exact, UPFD, HUPFD and EHUPFD methods when solving the non-linear ADR equation.
Table 5. Solutions obtained by the exact, UPFD, HUPFD and EHUPFD methods when solving the non-linear ADR equation.
(x, t) Exact SolutionUPFD SolutionAbsolute ErrorHUPFD SolutionAbsolute ErrorEHUPFD SolutionAbsolute Error
( 0.03338 , 0.002 ) 6.5475 × 10 9 5.4779 × 10 9 1.0697 × 10 9 5.4763 × 10 9 1.0712 × 10 9 5.4908 × 10 9 1.0567 × 10 9
( 0.0676 , 0.004 ) 8.1536 × 10 8 6.8787 × 10 8 1.2749 × 10 8 6.8751 × 10 8 1.2784 × 10 8 6.9355 × 10 8 1.2181 × 10 8
( 0.1014 , 0.006 ) 8.4478 × 10 7 7.1808 × 10 7 1.2670 × 10 7 7.1759 × 10 7 1.2718 × 10 7 7.2659 × 10 7 1.1819 × 10 7
Table 6. Infinity norm results obtained when solving the non-linear ADR equation by using the exact, UPFD, NSFD and HUPFD methods.
Table 6. Infinity norm results obtained when solving the non-linear ADR equation by using the exact, UPFD, NSFD and HUPFD methods.
(x, t) ExactUPFDErrorHUPFDErrorEHUPFDErrorCrank–NicolsonErrorNSFDError
u ( : , 0.002 ) 2.0590 2.0127 0.0512 2.0127 0.0512 2.0126 0.0513 2.0127 0.0512 2.0127 0.0512
u ( : , 0.004 ) 2.0591 2.0127 0.0512 2.0126 0.0514 2.0123 0.0514 2.0126 0.0514 2.0127 0.0512
u ( : , 0.004 ) 2.0592 2.0128 0.0512 2.0126 0.0515 2.0121 0.0515 2.0126 0.0515 2.0128 0.0512
Time 0.001227 0.004136 0.000492 0.004205 0.003954
Table 7. Convergence rate when solving the Schnakenberg model by using the UPFD, HUPFD and EHUPFD methods.
Table 7. Convergence rate when solving the Schnakenberg model by using the UPFD, HUPFD and EHUPFD methods.
qtk)
ΔtkΔxk UPFD HUPFD EHUPFD
0.00005 0.6667 1.0002 1.0029 1.0027
0.0001 0.6667 1.0004 1.0058 1.0055
0.0002 0.6667 1.0007 1.0116 1.0109
0.0004 0.6667
Table 8. Solutions u obtained by the UPFD, HUPFD and EHUPFD methods when solving the Schnakenberg system equation.
Table 8. Solutions u obtained by the UPFD, HUPFD and EHUPFD methods when solving the Schnakenberg system equation.
(x, t) UPFDHUPFDEHUPFD
u ( : , 0.0002 ) 0.0048 0.0050 0.0050
u ( : , 0.0004 ) 0.0182 0.0198 0.0198
u ( : , 0.0006 ) 0.0690 0.0788 0.0788
Time 0.005959 0.024321 0.004894
Table 9. Solutions v obtained by the UPFD, HUPFD and EHUPFD methods when solving the Schnakenberg system equation.
Table 9. Solutions v obtained by the UPFD, HUPFD and EHUPFD methods when solving the Schnakenberg system equation.
(x, t) UPFDHUPFDEHUPFD
u ( : , 0.0002 ) 0.0048 0.5310 0.5310
u ( : , 0.0004 ) 0.0186 0.1429 0.1429
u ( : , 0.0006 ) 0.0711 0.0390 0.0390
Time 0.005312 0.021700 0.004220
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. Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method. Mathematics 2024, 12, 1009. https://doi.org/10.3390/math12071009

AMA Style

Ndou N, Dlamini P, Jacobs BA. Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method. Mathematics. 2024; 12(7):1009. https://doi.org/10.3390/math12071009

Chicago/Turabian Style

Ndou, Ndivhuwo, Phumlani Dlamini, and Byron Alexander Jacobs. 2024. "Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method" Mathematics 12, no. 7: 1009. https://doi.org/10.3390/math12071009

APA Style

Ndou, N., Dlamini, P., & Jacobs, B. A. (2024). Solving the Advection Diffusion Reaction Equations by Using the Enhanced Higher-Order Unconditionally Positive Finite Difference Method. Mathematics, 12(7), 1009. https://doi.org/10.3390/math12071009

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