Novel numerical analysis for nonlinear advection–reaction–diffusion systems

Naveed Shahid; Nauman Ahmed; Dumitru Baleanu; Ali Saleh Alshomrani; Muhammad Sajid Iqbal; Muhammad Aziz-ur Rehman; Tahira Sumbal Shaikh; Muhammad Rafiq

doi:10.1515/phys-2020-0011

Article Open Access

Novel numerical analysis for nonlinear advection–reaction–diffusion systems

Naveed Shahid , Nauman Ahmed , Dumitru Baleanu , Ali Saleh Alshomrani , Muhammad Sajid Iqbal , Muhammad Aziz-ur Rehman , Tahira Sumbal Shaikh and Muhammad Rafiq

Published/Copyright: May 20, 2020

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

From the journal Open Physics Volume 18 Issue 1

Abstract

In this article, a numerical model for a Brusselator advection–reaction–diffusion (BARD) system by using an elegant numerical scheme is developed. The consistency and stability of the proposed scheme is demonstrated. Positivity preserving property of the proposed scheme is also verified. The designed scheme is compared with the two well-known existing classical schemes to validate the certain physical properties of the continuous system. A test problem is also furnished for simulations to support our claim. Prior to computations, the existence and uniqueness of solutions for more generic problems is investigated. In the underlying system, the nonlinearities depend not only on the desired solution but also on the advection term that reflects the pivotal importance of the study.

Keywords: advection; diffusion; reaction; consistency; stability; structure preserving; equi-continuity

1 Introduction

Nonlinear partial differential equations describe the various physical phenomena in applied sciences. So more researchers are directed to find the solutions of these equations in recent years. Nowadays, reaction diffusion systems of partial differential equations have launched multitudinous applications in chemical and biological phenomena to finance, physics, medicine, genetics, weather predictions and so on [1,2,3,4,5,6]. In this paper, we are engrossed to analyze the mathematical models, which comprise the mixture of advection, diffusion and reaction in the framework of partial differential equations. The advection–reaction–diffusion equation is one of the most pertinent areas in the applied mathematics, like chemical reactions in chemistry, biology, meteorology, epidemiology, fluid dynamics, and other fields of applied sciences. The inclusion of the combination of these three components in a mathematical model of partial differential equations puts a strong impact on the theory of partial differential equations and gives a rise to rethink the corresponding models in the light of the physical properties of a governing phenomenon [7,8,9,10,11]. Most of the mathematical models associated with the mixture of the reactive and advective processes can be observed in various directions in real life, for instance, in meteorological pollution control models, dynamics of age-structured population and problems consisting of the enhancement of oil recovery, etc. [12,13,14,15,16].

If u(x,t) and v(x,t) are the concentrations of two chemical species, then the terms ∂ u ∂ x and ∂ v ∂ x represent the advection of the concentration of quantities u and v, respectively, by the velocity field [17]. Clearly, the terms ∂ 2 u ∂ x 2 and ∂ 2 v ∂ x 2 represent the diffusion (random movement of species) of the concentration of u and v, respectively, in one dimension, which describe the random movement of Ref. [18]. Often, the exact solutions of an advection–reaction–diffusion problem are not available when such systems are nonlinear in nature, so one encounters a number of complications to gauge their approximations. We are then bound to find the approximate solutions of the corresponding problem. In spite of all these, we face great complications when the system of advection–reaction–diffusion equations is nonlinear.

This article is organized as follows: Section 2 presents the general problem of the advection–reaction–diffusion system and results on existence of a solution for the proposed model. Section 3 presents a Brusselator advection–reaction–diffusion (BARD) model. Here, we analyze the problem by known numerical schemes: upwind implicit and Crank Nicolson numerical scheme. In Section 4, we apply a newly constructed scheme on BADR model and simulate the results and observe the unconditional natural properties such as consistency and stability of our proposed scheme. Also, the results about positivity preserving algorithms are formulated. Section 5 presents an example confirming the stated results. Also, we observe the behavior of our constructed scheme with the comparison of the other two schemes discussed earlier in this section. Section 6 concludes this study with an overview of the obtained results.

2 Existence of solutions

In this contemprary section, the existence of the solution for the advection–reaction–diffusion system is discussed. A couple of results for the existence of the solution to the said system are established.

Consider a system of differential equations:

(1) ∂ u ∂ t + ∂ u ∂ x = d 1 ∂ 2 u ∂ x 2 + f 1 ( u , v )

(2) ∂ v ∂ t + ∂ v ∂ x = d 2 ∂ 2 v ∂ x 2 + f 2 ( u , v ) .

To understand the mathematical analysis of the advection–reaction–diffusion system system, which is a well-known first-order nonlinear dynamical system of advection–diffusion partial differential equations, we use the Schauder fixed point theorem for the existence of the solution. An important question in the existence theory is to guarantee if the advection–reaction–diffusion system system possesses some solutions. The junction of the fixed-point, operator theory helps us to reduce the initial value problem for th advection–reaction–diffusion system system in the corresponding fixed point operator. Since both the partial differential equations are of first order with respect to the time variable t, it is quite simple to invert the differential operator ∂ / ∂ t . But before we reduce the given system to the fixed-point problem, we rewrite equations (1) and (2) in a more compact form, that is,

(3) u t = F 1 ( u . v , ∂ u ∂ x , Δ x )

(4) v t = F 2 ( u . v , ∂ v ∂ x , Δ x ) .

Here, it is important to note that both u and v being the concentrations of the quantities should be nonnegative. Importantly, the functions F ₁ and F ₂ on the right-hand side of equations (3) and (4) can be nonlinear in not only the solution pair (u,v) but also the the first- and second-order space derivatives u _x, v _x, u _xx and v _xx of the desired solution pair.

In view of the operator theory, equations (1) and (2) can be written in the following form:

(5) u ( x , t ) = u 0 + ∫ 0 t F 1 ( u , v , u x , u x x ) ( τ , x ) d t

(6) v ( x , t ) = v 0 + ∫ 0 t F 2 ( u , v , v x , v x x ) ( τ , x ) d t

with setting F ₂(u, v, v _x, v _xx) ≡ F ₂.

Such inversion of the first-order partial derivatives is obvious, but for the inversion of more general operators, one may need special kernel functions. It is clear from equations (5) and (6) that the solutions u and v depend on the same functions appearing in the right-hand sides of the same equations, respectively. So the operator theory allows us to write equations (5) and (6) in an operator form [19,20,21]. But in the case of a system of equations under the contemporary study, we can rewrite them into the following single equation.

(7) u = u 0 + ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ .

The integral equation (7) can be written as a fixed point problem:

(8) S u ( t ) = u ( t ) ,

where

( u ( t ) , v ( t ) ) ∈ C 2 [ a , b ] × C 1 [ 0 , T ] = X

and let S:E → E be a self-mapping defined by

(9) S u ( t ) = u 0 + ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ

where E is supposed to be a Banach space.

We now establish a lemma, which shows the compactness of S.

Lemma 2.1

Let S:E → E be any map defined by

(10) S u ( t ) = u 0 + ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ

with u(0) = u ₀ and suppose that

k = max { F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) }

with 0 ≤ u ≤ N ₁, 0 ≤ v ≤ N ₂ , where N ₁ and N ₂ are finite positive numbers. Then, S is relatively compact.

Proof

Since both F ₁ and F ₂ are assumed to be nonnegative and continuous functions, so, the operator S is continuous. Let B be a bounded subset of E and there exists a number N ₀ such that for any u ∈ B, we have

∥ u ∥ ≤ N 0 .

Now, take

∥ S u ∥ = | u 0 + ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ | ≤ | u 0 | + | ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ | ≤ | u 0 | + | ∫ 0 t k d τ | = | u 0 | + t k

(11) ⇒ ∥ S u ∥ ≤ | u 0 | + t k

hence, S(B) is uniformly bounded.

Next, we have to show the equi-continuity of S. For each u ∈ B and for ε > 0 and t ₁,t ₂ ∈ [0,T] such that t ₁ < t ₂ then |t ₂ − t ₁| < δ and let δ = ε k .

For this, let

| S i u ( x , t 1 ) − S i u ( x , t 2 ) | = | ∫ 0 t 1 F 1 ( u i , v i 0 + ∫ 0 τ F 2 ( s ) d s , u i x , u i x x ) d τ − ∫ 0 t 2 F 1 ( u i , v i 0 + ∫ 0 τ F 2 ( s ) d s , u i x , u i x x ) d τ | = | ∫ 0 t 1 F 1 ( u i , v i 0 + ∫ 0 τ F 2 ( s ) d s , u i x , u i x x ) d τ − ∫ 0 t 1 F 1 ( u i , v i 0 + ∫ 0 τ F 2 ( s ) d s , u i x , u i x x ) d τ − ∫ t 1 t 2 F 1 ( u i , v i 0 + ∫ 0 τ F 2 ( s ) d s , u i x , u i x x ) d τ | = | ∫ t 1 t 2 F 1 ( u i , v i 0 + ∫ 0 τ F 2 ( s ) d s , u i x , u i x x ) d τ | = k | t 2 − t 1 | .

Clearly, |S _i u(x,t ₁) − S _i u(x,t ₂)| approaches zero as |t ₁ − t ₂| → 0. This implies that S _i is equicontinuous. By the Arzela Ascoli theorem, there exists a subsequence S i j ⊆ S i , so S i j converges uniformly to a point x* ∈ C[0,T]. This shows that S(B) is equicontinuous. Therefore, S is relatively compact.

Also, as S is in space of continuous functions, S:C[0,T] → C[0,T] is self-map. F ₁ and F ₂ are nonnegative, and u ₀ ≥ 0 because u ₀ represents the initial concentration that cannot be negative. That is,

S u ( x , t ) = u 0 + ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ ≥ 0 .

Hence, according to the statement of the Schauder fixed-point theorem, the operator S has at least one fixed point u(x,t) ∈ E, where t ∈ [0,T]. Hence, equation (1) has at least one solution u(x,t) in C ²[0,T].

Now we can establish a theorem for the existence of the solution of equation (7).

Theorem 2.2

If u, v are twice continuously differentiable with respect to the space variable x and continuously differentiable with respect to time t, then the initial boundary value problem for equations (1) and (2) are solvable by the Schauder fixed point theorem for any continuous initial values.

2.1 Uniqueness of the solution

Theorem 2.2 guarantees the existence of the solution, which itself gives importance to the considered problem for physical purpose. For the uniqueness of the solution of equations (1) and (2),

S ≡ S u , v ( x , t ) = u 0 + ∫ 0 t F 1 ( u , v 0 + ∫ 0 τ F 2 ( s ) d s , u x , u x x ) d τ .

Consider

S 1 − S 2 = ∫ 0 t { F 1 1 ( u 1 , v 0 + ∫ 0 τ F 2 1 ( s ) d s , u x 1 , u x x 1 ) − F 1 2 ( u 2 , v 0 + ∫ 0 τ F 2 2 ( s ) d s , u x 2 , u x x 2 ) } d τ .

Now suppose F ₁ is Lipschitz continuous in u, u _x, u _x x and ∫ 0 t F 2 ( s ) d s . The operator S turns out to be contractive under certain conditions on the bounds of the functions F ₁ and F ₂. Suppose that the absolute value of the integrand is bounded by

L ρ ∥ u 1 − u 2 ∥ C 2 ( Ω ¯ ) × c 2 [ 0 , T ]

where ρ is the length of interval of continuity. Then for the contractivity of S guarantees if

L ρ < 1 L < 1 ρ , ρ = T .

Finally, we establish the following theorem with unique existence of solution.

Theorem 2.3

Let u , v ∈ C 2 ( Ω ¯ ) × c 1 [ 0 , T ] and then the initial value problem (1) and (2) is uniquely solvable by contraction mapping principle, provided the functions F ₁ and F ₂ satisfy the Lipschitz condition in functions u, v, u _x, v _x, u _xx, v _xx.

3 Numerical methods

The primary aim of the study of finite difference schemes to find the solution of linear and nonlinear partial differential equations is to discretize the given continuous system by approximating the partial derivatives occurring in the continuous formulations by the finite number of function values at some selected finite number of points in the domain. In this respect, Taylor’s series is the finest way to obtain these approximations.

For the rest of the paper, let K and L be positive integers and τ be any positive real number. To find the approximate solution of system of equations (12) and (13) in the spatial interval [a,b] over the time period [0,τ], we make a partitions a = x ₀ < x ₁ < x ₂ < … < x _L = b and 0 = t ₀ < t ₁ < t ₂ < … < t _K = τ of [a,b] and [0,τ] respectively, with the norm Δ x = b − a L and Δ t = τ K .

Let u i n be the approximate value of the function u at a grid point (x _i,t _n).

The generalized advection–reaction–diffusion model is as follows:

(12) ∂ u ∂ t + ∂ u ∂ x = d 1 ∂ 2 u ∂ x 2 + B − ( A + 1 ) u + u 2 v

(13) ∂ v ∂ t + ∂ v ∂ x = d 2 ∂ 2 v ∂ x 2 + A u − u 2 v

with initial and boundary conditions, i.e.,

u ( x , 0 ) = f ( x ) v ( x ,0) = g ( x ) ∂ u (0, t ) ∂ x = ∂ u (1, t ) ∂ t = 0 ∂ v (0, t ) ∂ x = ∂ v (1, t ) ∂ t = 0

where u and v are the concentration of two different reactants, A and B are constant terms representing constant concentration during reaction operation, d ₁ and d ₂ are diffusion coefficients and all A, B, d ₁, and d ₂ are positive constants. The equilibrium point of the system (12) and (13) is ( u ⁎ , v ⁎ ) = ( B , A B ) . Like most of the other mathematical models, it is quit hard to find the exact solution of the system (12) and (13), so we use some numerical methods to solve it. It should be kept in mind that these numerical schemes are effective, which preserve the basic physical properties of the model. For instance, positivity preserving and structural preserving are the most important physical properties of the numerical schemes. Mickens in Ref. [22] described the criteria of developing the structural preserving finite difference schemes named as nonstandard finite difference scheme. Twizell et al. in Ref. [23] examined the attractive fixed points for the Brusselator model and established a condition for the convergence of the concentration profile (u,v) to (u*,v*), that is, (u,v) → (u*,v*), which is 1 − A + B ² > 0. In this paper, we use a finite difference scheme developed by Chen-Charpentier and Kojouharov [24], which is unconditionally positivity preserving the scheme for the solution of advection reaction–diffusion equation. To show the efficiency of our proposed scheme, we apply some other reliable methods discussed in the literature and compare the results of these methods with the results developed by the proposed scheme.

3.1 The upwind implicit scheme

According to this scheme, both time and space derivatives are established as follows:

∂ u ∂ t = u i n + 1 − u i n k ∂ u ∂ x = u i n + 1 − u i − 1 n + 1 h ∂ 2 u ∂ x 2 = u i + 1 n + 1 − 2 u i n + 1 + u i − 1 n + 1 h 2 .

Substituting theses values in equations (12) and (13), we get the upwind implicit scheme as follows:

− ( ξ 1 + τ ) u i − 1 n + 1 + ( 1 + τ + 2 ξ 1 ) u i n + 1 − ξ 1 u i + 1 n + 1 = ( 1 − k ( A + 1 ) ) u i n + k ( u i n ) 2 v i n + B k − ( ξ 2 + τ ) v i − 1 n + 1 + ( 1 + τ + 2 ξ 2 ) v i n + 1 − ξ 2 v i + 1 n + 1 = ( 1 − k A − k ( u i n ) 2 ) v i n

where τ = k h , ξ 1 = k d 1 h 2 and ξ 2 = k d 2 h 2 .

The aforementioned upwind implicit schemes indicate that the system of algebraic equations developed by this scheme is diagonally dominant and so its solution exists. The behavior of this scheme is found to be Von Neumann stable.

3.2 The Crank Nicolson scheme

In this scheme, we replace both the time and the spatial partial derivatives by their central difference approximations at a grid point ( i h , 2 n + 1 2 k ) , which is the center point of the consecutive grid points (i,n) and (i,n + 1). That is,

∂ u ∂ t | i n = u i n + 1 − u i n k ∂ u ∂ x | i n = u i + 1 n − u i − 1 n 4 h + u i + 1 n + 1 − u i − 1 n + 1 4 h ∂ 2 u ∂ x 2 | i n = 1 2 ( u i + 1 n − 2 u i n + u i − 1 n h 2 + u i + 1 n + 1 − 2 u i n + 1 + u i − 1 n + 1 h 2 ) .

The difference equations associated by equations (12) and (13), after using the above approximations, we get,

( τ 4 − ξ 1 2 ) u i + 1 n + 1 + ( 1 + ξ 1 ) u i n + 1 − ( τ 4 + ξ 1 2 ) u i − 1 n + 1 = ( ξ 1 2 − τ 4 ) u n + 1 n + ( 1 + ξ 1 ) u i n + ( τ 4 + ξ 1 2 ) u i − 1 n + k B − k ( A + 1) u i n + k ( u i n ) 2 v i n ( τ 4 − ξ 2 2 ) v i + 1 n + 1 + ( 1 + ξ 2 ) v i n + 1 − ( τ 4 + ξ 2 2 ) v i − 1 n + 1 = ( ξ 2 2 − τ 4 ) v i + 1 n + ( 1 + ξ 2 ) v i n + ( τ 4 + ξ 2 2 ) v i − 1 n + A k v i n − ( u i n ) 2 v i n

where τ = k h and ξ 1 = k d 1 h 2 , ξ 2 = k d 2 h 2 .

Now, we are adopting an efficient finite difference scheme, which has positivity preserving and structural preserving properties. We apply this scheme on the model (12) and (13).

4 Proposed structure preserving numerical scheme

Numerous numerical schemes are developed and used to solve the mathematical models so far, and some of them are implicit and some are explicit schemes. In explicit schemes, dependent variables are expressed as a function of some known quantities at the previous time unit (say n time step), whereas in case of implicit schemes, dependent variables are determined by the coupled system of multiple simultaneous algebraic equations and used to obtain the solution either in the form of the matrix or some iterative processes. In this method, all unknown quantities are evaluated at the future time step, (say (n + 1) time step), as well as we are allowed to take a large time step size in each iteration of the implicit scheme. Although computations through implicit schemes are intensively expensive due to their formations, there is only less error in the simulation process to explicit methods. As stability represents the behavior of the solution with the increase of the time step size, if the solution shows the controlled behavior for a large time step size, the numerical scheme is unconditionally stable. But such a situation does not appear in case of explicit schemes, which are usually always conditionally stable. In this paper, we use an explicit scheme to solve the system (12) and (13). But we observe that this proposed scheme, despite of being explicit, is very effective and behaved well, consistent and unconditionally stable. However, other implicit schemes in this paper do not behave well. We demonstrate such a situation in this paper through simulation.

The proposed scheme, which we are adopting for the model (12) and (13) developed in Ref. [24], is explicit in nature and unconditionally positivity preserving.

The formulas for temporal and spatial derivatives according to our proposed scheme are given by

∂ u ∂ t = u i n + 1 − u i n k ∂ u ∂ x = u i n + 1 − u i − 1 n h ∂ 2 u ∂ x 2 = u i − 1 n − 2 u i n + 1 + u i + 1 n h 2 .

The numerical modulation for equations (12) and (13) are as follows:

(14) u i n + 1 = ξ 1 u n + 1 n + ( ξ 1 + τ ) u i − 1 n + u i n + k B + k ( u i n ) v i n u i n 1 + τ + 2 ξ 1 + k ( A + 1 )

(15) v i n + 1 = v i n + τ v i − 1 n + ξ 2 { v i − 1 n + v i + 1 n } + A k u i n 1 + τ + 2 ξ 2 + k ( u i n ) 2

where τ = k h , ξ 1 = k d 1 h , ξ 2 = k d 2 h 2 .

We will check the consistency and the stability of our proposed scheme by equations (12) and (13).

Theorem 4.1

Suppose that u i n and v i n are the approximate values of the functions u and v, respectively, at a grid point (x _i,t _n). Then, the proposed scheme (14) and (15) assures the positivity, in whole domain, of positive solutions under the hypothesis of positive initial functions, That is,

u i n ≥ 0 , v i n ≥ 0 ⇒ u i n + 1 ≥ 0 , v i n + 1 ≥ 0 .

4.1 Remark

Equations (14) and (15) give the positivity of the solution because the known functions occurred in the initial conditions are nonnegative, so the right-hand sides of equations (14) and (15) bear no negative terms for all discretization parameters i, n = 0, 1, 2, 3,….

4.2 Consistency of proposed scheme

For equation (12)

∂ u ∂ t + ∂ u ∂ x = d 1 ∂ 2 u ∂ x 2 + B − ( A + 1 ) u + u 2 v .

By using the approximations of the proposed scheme, we obtain

(16) u i n + 1 − u i n + k h ( u i n + 1 − u i − 1 n − 1 ) = d 1 k h 2 ( u i − 1 n − 2 u i n + 1 + u i + 1 n ) + k B − k ( A + 1) u i n + 1 + k ( u i n ) 2 v i n ( 1 + k h + 2 2 d 1 k h 2 + k ( A + 1) ) u i n + 1 = d 1 k h 2 u i + 1 n + ( k h + d k h 2 ) u i − 1 n + u i n + k B + k ( u i n ) 2 v i n

By using Taylor’s series and simplifications, we obtain

( 1 + k h + 2 d 1 k h 2 + k ( A + 1 ) ) ( ∂ u ∂ t + k 2 ! ∂ 2 u ∂ t 2 + ⋯ ) = d 1 ( 2 2 ! ∂ 2 u ∂ t 2 + 2 h 2 4 ! ∂ 4 u ∂ t 4 + ⋅ ⋅ ⋅ ) − ( ∂ u ∂ x − h 2 ! ∂ 2 u ∂ x 2 + h 2 3 ! ∂ 3 u ∂ x 3 + ⋯ ) − ( A + 1) u i n + B + ( u i n ) 2 v i n

By replacing k by h ³, we obtain

( 1 + h 3 h + 2 d 1 h 3 h 2 + h 3 ( A + 1 ) ) ( ∂ u ∂ t + h 3 2 ! ∂ 2 u ∂ t 2 + ⋯ ) = d 1 ( 2 2 ! ∂ 2 u ∂ t 2 + 2 h 2 4 ! ∂ 4 u ∂ t 4 + ⋯ ) − ( ∂ u ∂ x − h 2 ! ∂ 2 u ∂ x 2 + h 2 3 ! ∂ 3 u ∂ x 3 ) − ( A + 1 ) u i n + B + ( u i n ) 2 v i n

By taking the limit as h → 0

∂ u ∂ t + ∂ u ∂ x = d 1 ∂ 2 u ∂ x 2 + B − ( A + 1 ) u + u 2 v

which proves the consistency.

Now from equation (13)

(17) v i n + 1 − v i n k + v i n + 1 − v i − 1 n h = d 2 v i − 1 n − 2 v i n + 1 + v i + 1 n h 2 + A u i n − ( u i n + 1 ) 2 ( v i n + 1 )

v i n + 1 − v i n + k h ( v i n + 1 − v i − 1 n ) = k d 2 h 2 ( v i − 1 n − 2 v i n + 1 + v i + 1 n ) + A k u i n − k ( u i n ) 2 ( v i n + 1 )

v i n + 1 + k h v i n + 1 + 2 k d 2 h 2 v i n + 1 + k ( u i n ) 2 ( v i n + 1 ) = v i n + k h v i − 1 n + k d 2 h 2 ( v i − 1 n + v i + 1 n ) + A k u i n

( 1 + k h + 2 k d 2 h 2 + k ( u i n ) 2 ) v i n + 1 = v i n + k h v i − 1 n + 2 k d 2 h 2 ( v i − 1 n + v i + 1 n ) + A k u i n

Taylor’s formulae are as follows:

v i n + 1 = v i n + k ∂ u ∂ t + k 2 2 ! ∂ 2 u ∂ t 2 + k 3 3 ! ∂ 3 u ∂ t 3 + ⋯ v i − 1 n = v i n − h ∂ v ∂ x + h 2 2 ! ∂ 2 v ∂ x 2 − h 3 3 ! ∂ 3 v ∂ x 3 + ⋯ v i + 1 n = v i n + h ∂ v ∂ x + h 2 2 ! ∂ 2 v ∂ x 2 + h 3 3 ! ∂ 3 v ∂ x 3 + ⋯

After substituting these formulas in the aforementioned equation, we obtain

( 1 + k h + 2 k d 2 h 2 + k ( u i n ) 2 ) ( ∂ u ∂ t + k 2 ! ∂ 2 u ∂ t 2 + k 2 3 ! ∂ 3 u ∂ t 3 + ⋯ ) = − ∂ v ∂ x + h 2 ! ∂ 2 v ∂ x 2 − h 2 3 ! ∂ 3 v ∂ x 3 + h 3 4 ! ∂ 4 v ∂ x 4 + ⋯ + d 2 ( 2 2 ! ∂ 2 v ∂ x 2 + 2 h 2 4 ! ∂ 4 v ∂ x 4 + ⋯ ) + A u i n − ( u i n ) 2 v i n

By replacing k by h ³ and h → 0, we obtain

∂ u ∂ t = − ∂ v ∂ x + d 2 ∂ 2 v ∂ x 2 + A u − u 2 v .

Which is (13). It verifies that the proposed scheme is consistent.

Now we investigate the stability of the proposed scheme referred to given system of equations.

4.3 Stability of proposed scheme

To verify the stability of our proposed scheme, we consider both the equations of given mathematical model of advection–reaction–diffusion. First, we select equation (12) and incorporate τ = k h , ξ 1 = d 1 k h 2 in the associated equation and after linearizing it by equation (16), we obtain

u i n + 1 − u i n + τ ( u i n + 1 − u i − 1 n ) = ξ 1 ( u i − 1 n − 2 u i n + 1 + u i + 1 n ) + k B − k ( A + 1) u i n + 1 + k ( u i n ) 2 v i n

( 1 + τ + 2 ξ 1 + k ( A + 1 ) ) u i n + 1 = ξ 1 u n + 1 n + ( ξ 1 + τ ) u i − 1 n + u i n + k B + k ( u i n ) v i n u i n

Let

u i n = ψ ( t ) e i α x u i n + 1 = ψ ( t + Δ t ) e i α x u i + 1 n = ψ ( t ) e i α ( x + Δ x ) u i − 1 n = ψ ( t ) e i α ( x − Δ x ) .

For the sake of linearization, let v i n = 1 , B = 0 , we have

( 1 + τ + 2 ξ 1 + k ( A + 1 ) ) ψ ( t + Δ t ) e i α x = ξ 1 ψ ( t ) e i α ( x + Δ x ) + ψ ( t ) e i α ( x + Δ x ) + ψ ( t ) e i α x + k ψ ( t ) e i α x .

This implies that

| ψ ( t + Δ t ) ψ ( t ) | ≤ 2 ξ 1 + τ + 1 + k 1 + τ + 2 ξ 1 + k ( A + 1 )

which is always less than 1.

Now set τ = k h , ξ 2 = k d 2 h 2 in equation (17), we obtain

v i n + 1 + τ v i n + 1 + 2 ξ 2 v i n + 1 + k ( u i n ) 2 ( v i n + 1 ) = v i n + τ v i − 1 n + ξ 2 ( v i − 1 n + v i + 1 n ) + A k u i n .

Let

v i n = ϕ ( t ) e i α x v i n + 1 = ϕ ( t + Δ t ) e i α x v i + 1 n = ϕ ( t ) e i α ( x + Δ x ) v i − 1 n = ϕ ( t ) e i α ( x − Δ x )

( 1 + τ + 2 ξ 2 + k ) ϕ ( t + Δ t ) e i α x = ϕ ( t ) e i α x + τ ϕ ( t ) e i α ( x − Δ x ) + ξ 2 ( ϕ ( t ) e i α ( x − Δ x ) + ϕ ( t ) e i α ( x + Δ x ) )

( 1 + τ + 2 ξ 2 + k ) ϕ ( t + Δ t ) = ϕ ( t ) ( 1 + τ e − i α x + ξ 2 ( e − i α x + e i α x ) )

( 1 + τ + 2 ξ 2 + k ) ϕ ( t + Δ t ) = ϕ ( t ) ( 1 + τ cos α Δ x − i τ sin α Δ x + 2 ξ 2 cos α Δ x )

( 1 + τ + 2 ξ 2 + k ) | ϕ ( t + Δ t ) ϕ ( t ) | = | ( 1 + ( τ + 2 ξ 2 ) cos α Δ x − i τ sin α Δ x ) |

| ϕ ( t + Δ t ) ϕ ( t ) | = ( 1 + ( τ + 2 ξ 2 ) cos α Δ x ) 2 + ( τ sin α ) 2 1 + τ + 2 ξ 2 + k ≤ 1

According to Von Neumann stability criteria, our proposed scheme is unconditionally stable.

5 Numerical example and simulations

(18) ∂ u ∂ t + ∂ u ∂ x = d 1 ∂ 2 u ∂ x 2 + B − ( A + 1 ) u + u 2 v

(19) ∂ v ∂ t + ∂ v ∂ x = d 2 ∂ 2 v ∂ x 2 + A u − u 2 v

with boundary value conditions

u ( x , 0 ) = 1 / 2 v ( x , 0 ) = 1 + 5 x ∂ u ( 0 , t ) ∂ x = ∂ u ( 1 , t ) ∂ t = 0 .

Clearly, in Figure 1(b) and (c) drawn from upwind implicit scheme, by choosing A and B such that 1 − A + B ² > 0, we observe that the graph is lying in the negative side, which is not possible for any type of concentration.

Figure 1

(a) Mesh graph of u (concentration profile) using upwind implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.006, τ = 3 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using the upwind implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.006, τ = 3 and d ₁ = d ₂ = 10⁻⁴. (c) Combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

On the other hand, in Figure 2, when we observe the graph of concentrations u and v from the proposed scheme at the same values of the parameters, it remains in the positive side, which shows that the proposed scheme is well behaved or positivity preserving with the aforementioned parameter values.

Figure 2

(a) Mesh graph of u (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.006, τ = 3 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.006, τ = 3 and d ₁ = d ₂ = 10⁻⁴. (c) Combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

Now we look at the comparison of proposed scheme with Crank Nicolson implicit scheme.

In Figure 3, after choosing the values of A and B such that 1 − A + B ² > 0, the graphs of u and v fall in the negative side, which is naturally not possible.

Figure 3

(a) Mesh graph of u (concentration profile) using Crank Nicolson implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.032, τ = 1.6 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using Crank Nicolson implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.008, τ = 1.6 and d ₁ = d ₂ = 10⁻⁴. (c) Combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

In Figure 4, the graph of the proposed scheme with same values of parameters used in Crank Nicolson implicit scheme shows the positivity of the concentration variables, that is, the scheme is positivity preserving.

Figure 4

(a) Mesh graph of u (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.032, τ = 1.6 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.008, τ = 1.6 and d ₁ = d ₂ = 10⁻⁴. (c) The combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

In Figure 5, when we incorporate the same values of ξ ₁ and ξ ₂ in both concentration profiles with comparatively large time scale, we obtain the graphs of u and v that are tilted in the negative sides. This shows the unstability of the scheme.

Figure 5

(a) Mesh graph of u (concentration profile) using upwind implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.012, τ = 5.90 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using upwind implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.012, τ = 5.90 and d ₁ = d ₂ = 10⁻⁴. (c) Combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

But in the graph of proposed scheme, with the same values of the parameters, we observe the stability of our proposed scheme (Figure 6).

Figure 6

(a) Mesh graph of u (concentration profile) using proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.012, τ = 5.90 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.012, τ = 5.90 and d ₁ = d ₂ = 10⁻⁴. (c) Combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

In Figure 7, when we incorporate the same values to ξ ₁ = ξ ₂ = 0.004, with same τ = 2, the well-known Crank Nicolson scheme gives the unstable behavior.

Figure 7

(a) Mesh graph of u (concentration profile) using Crank Nicolson implicit FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.004, τ = 2 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using Crank Nicolson FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.004, τ = 2 and d ₁ = d ₂ = 10⁻⁴. (c) Combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

On the other hand, our proposed scheme shows the stability at the same values of parameters used in the graph of Figure 8.

In Figure 9, when the condition for convergence of concentration profile is violated, that is, even our scheme shows the positivity, which is one of the basic requirements of any numerical scheme. Finally, we observe that our proposed scheme for the proposed model is unconditionally stable.

Figure 8

(a) Mesh graph of u (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.004, τ = 2 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.004, τ = 2 and d ₁ = d ₂ = 10⁻⁴. (c) The combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

Figure 9

(a) Mesh graph of u (concentration profile) using the proposed FD method at h = 0.05, B = 3.4, A = 1, ξ ₁ = ξ ₂ = 0.004, τ = 2 and d ₁ = d ₂ = 10⁻⁴. (b) Mesh graph of v (concentration profile) using the proposed FD method at h = 0.05, B = 1, A = 3.4, ξ ₁ = ξ ₂ = 0.004, τ = 2 and d ₁ = d ₂ = 10⁻⁴. (c) The combined plot of u, v at x = 1 vs time with same parameters as discussed in (a) and (b).

6 Conclusion

In this paper, we demonstrated the existence of solution of a nonlinear advection–reaction–diffusion model by using the Schauder fixed point theorem, which is a strong tool to prove the existence of fixed points of nonlinear operators in the fixed point theory. First, we converted the system of equations to a single nonlinear differential equation and then constructed a fixed-point differential operator. We proved here that the fixed points of this operator exist under some suitable conditions. After obtaining the existence of solution of our proposed model, we designed an unconditionally positivity preserving finite difference numerical scheme, which is explicit in nature, for a nonlinear advection–reaction–diffusion model. Despite of explicitly, our proposed scheme is more efficient than the existing techniques due to positivity preserving property. The results of the proposed scheme are compared with the well-known upwind implicit scheme and Crank Nicolson scheme and can be observed that these schemes are unable to preserve the positivity at certain time interval. Since the solutions functions of the differential equations are not globally bounded, it is better to consider the closed optimal balls in function spaces, which give the explicit estimates for the solutions of the system. So for the future work, we can optimize the convergence of the solution. In case of the numerical analysis, we can use our proposed scheme in two- or three-dimensional advection–reaction–diffusion system to obtain the numerical solution for the future perspective. Also, the work can be extended to the system of stochastic differential equations.

dumitru.baleanu@gmail.com

Conflict of interest: The authors declare that they have no conflict of interest.
Funding: This study has no funding.

References

[1] Arson DG. A comparison method for stability analysis of nonlinear parabolic problems. SIAM Rev. 1978;20:245–64.10.1137/1020038Search in Google Scholar

[2] Bebemes JW, Schmitt K. On the existence of maximal and minimal solutions for parabolic partial differential equations. Pro Amer Math Soc. 1979;73(21):1–218.Search in Google Scholar

[3] Escobedo M, Herrero MA. Boundedness and blow up for a semilinear reaction-diffusion system. J Differential Equ. 1991;89:176–202.10.1016/0022-0396(91)90118-SSearch in Google Scholar

[4] Liu BC, Dong MZ, Li FJ. Asymptotic properties of blow-up solutions in reaction-diffusion equations with nonlocal boundary flux. Z Angew Math Phys. 2018;69:22.10.1007/s00033-018-0920-2Search in Google Scholar

[5] Ahmed N, Rafiq M, Rehman MA, Iqbal MS, Ali M. Numerical modeling of three dimensional Brusselator system. AIP Adv. 2019;9:015205.10.1063/1.5070093Search in Google Scholar

[6] He X, Ni WM. The effects of diffusion and spatial variation in Lotka-Volterra competition-diffusion system II: the general case. J Differential Equ. 2013;254:4088–108.10.1016/j.jde.2013.02.009Search in Google Scholar

[7] Zhou P, Xiao D. Global dynamics of a classical Lotka-Volterra competition-diffusion-advection system. J Funct Anal. 2018;275(2):356–80.10.1016/j.jfa.2018.03.006Search in Google Scholar

[8] Zhao XQ, Zhou P. On a Lotka-Volterra competition model: the effects of advection and spatial variation. Calc Var. 2016;55:137.10.1007/s00526-016-1021-8Search in Google Scholar

[9] Rui P, Zhao XQ. A nonlocal and periodic reaction-diffusion-advection model of a single phytoplankton species. J Math Biol. 2016;72:755–91.10.1007/s00285-015-0904-1Search in Google Scholar PubMed

[10] Yu X, Zhao XQ. Propagation phenomena for a reaction-advection-diffusion competition model in a periodic habitat. J Dyn Diff Equat. 2017;29:41, 10.1007/s10884-015-9426-1 Search in Google Scholar

[11] Zhang L, Wanga ZC, Zhao XQ. Threshold dynamics of a time periodic reaction-diffusion epidemic model with latent period. J Differential Equ. 2015;258(9):3011–36.10.1016/j.jde.2014.12.032Search in Google Scholar

[12] Lin H, Wang F. On a reaction-diffusion system modeling the Dengue transmission with nonlocal infections and crowding effects. Appl Math Comput. 2014;248:184–94.10.1016/j.amc.2014.09.101Search in Google Scholar

[13] Macufa MM, Mayer JF, Krindges A. Difussion of Malaria in Mozambique. Model Comput Simul, Biomatematica. 2015;25:161–84.Search in Google Scholar

[14] Lotfi EM, Maziane M, Hattaf K, Yousti N. Partial differential equations of an epidemic model with spatial difussion. Int J Partial Differential Equ. 2014;186437.10.1155/2014/186437Search in Google Scholar

[15] Ren X, Tian Y, Liu L, Liu X. A reaction difussion withen-host HIV model with cell-to-cell transmision. J Math Biol. 2018;76:1831–72.10.1007/s00285-017-1202-xSearch in Google Scholar PubMed

[16] Cosner C. Reaction-diffusion-advection models for the effects and evolution of dispersal. Discret Contin Dyn Syst. 2014;34:1701–45.10.3934/dcds.2014.34.1701Search in Google Scholar

[17] Cantrell RS, Cosner C. Spatial ecology via reaction-diffusion equations. Wiley Ser Math Comput Biol. Chichester, UK: John Wiley and Sons; 2003.10.1002/0470871296Search in Google Scholar

[18] Dockery J, Hutson V, Mischaikow K, Pernarowski M. The evolution of slow dispersal rates: areaction-diffusion model. J Math Biol. 1968;37:61–83.10.1007/s002850050120Search in Google Scholar

[19] Hsieh PF, Sibuya Y. Theory of ordinary differential equations. New York: Springer; 1999.10.1007/978-1-4612-1506-6Search in Google Scholar

[20] Agarwal RP, Lakshmikantham V. Uniqueness and nonuniqueness criteria for ordinary differential equations. Singapore: World Scientific; 1993.10.1142/1988Search in Google Scholar

[21] Djebali S, Górniewicz L, Ouahab A. Solution sets for differential equations and inclusions. Berlin/Boston: Walter de Gruyter GmbH; 2013.10.1515/9783110293562Search in Google Scholar

[22] Mickens RE. Nonstandard finite difference models of differential equations. Singapore: World Scientific; 1994. 10.1142/2081 Search in Google Scholar

[23] Twizell EH, Gumel AB, Cao Q. A second order scheme for the Brusselator reaction-diffusion system. J Math Chem. 1999;26:297–316.10.1023/A:1019158500612Search in Google Scholar

[24] Chen-Charpentier BM, Kojouharov HV. Unconditionally positive preserving scheme for advection-diffusion-reaction equations. Math Comput Modeling. 2013;57:2177–185.10.1016/j.mcm.2011.05.005Search in Google Scholar

Received: 2019-08-11

Accepted: 2020-02-11

Published Online: 2020-05-20

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

Articles in the same Issue

https://doi.org/10.1515/phys-2020-0011

Keywords for this article

advection; diffusion; reaction; consistency; stability; structure preserving; equi-continuity

Creative Commons

BY 4.0