Numerical study and stability of the Lengyel–Epstein chemical model with diffusion

*Correspondence: zahir.sha@kmutt.ac.th; poom.kum@kmutt.ac.th 2Center of Excellence in Theoretical and Computational Science (TaCS-CoE), SCL 802 Fixed Point Laboratory, Science Laboratory Building, King Mongkut’s University of Technology Thonburi (KMUTT), 126 Pracha-Uthit Road, Bang Mod, Thrung Khru, Bangkok 10140, Thailand 4KMUTT Fixed Point Research Laboratory, Room SCL 802 Fixed Point Laboratory, Science Laboratory Building, Department of Mathematics, Faculty of Science, King Mongkut’s University of Technology Thonburi (KMUTT), 126 Pracha-Uthit Road, Bang Mod, Thrung Khru, Bangkok 10140, Thailand Full list of author information is available at the end of the article Abstract


Introduction
The chemical reactions as the Belousov-Zhabotinsky reaction and the Briggs-Rauscher reaction are exceptional. The mathematical specimens of these chemical reactions are scrutinize mathematically; then again; these replicas are cluttered. The Lengyel-Epstein reaction associating iodine (I 2 ), malonic acid (MA), and chlorine dioxide (ClO 2 ) is a straightforward reaction as follows [1]: The iodide ion oxidation by free chlorine dioxide radical is given in the second equation of (1). (c) A reaction of iodide and chlorite ions created in the first and second equation of (1) generates iodine as given in the third equation of (1). The rate equations for the so-called ClO 2 -I 2 -MA (chlorine dioxide-iodine-malonic acid) are given by - and where k 1 , k 3 , k 4 , and k 5 are reaction rate constants, whereas k 2 and α denote saturation levels. Furthermore, the last term in (iii) characterizes the autocatalytic consequence of I 2 and the self-inhibitory outcome of Ion the chlorite-iodide reaction [2]. This framed term disappears when [I -] → 0, where no reaction can happen since no iodide is accessible, and in the limit and [MA]. Nonetheless, the iodization of Malonic acid helps largely as a cradle of iodine ions, and MA can be swapped by ethyl acetoacetate [2]. Furthermore, it is experimentally seen that the concentrations of chlorite and iodide ions fluctuate over numerous orders of level through an oscillation, whereas the concentrations of chlorine dioxide and malonic acid vary slowly. These concentrations may consequently be considered as constants, and actions of the system may be estimated by a two-variable specimen including only the concentrations of iodide and chlorine ions. For a flow reactor with appropriate serving, it is conceivable to retain the concentrations of malonic acid, chlorine dioxide and iodine approximately constant, and oscillations can still be perceived in the suitable ranges of temperature and concentrations. Thus, we conclude that MA, ClO 2 and I 2 vary much more sluggishly than the intermediate Arguing as in [3], an improved application in scientific modeling, it is vital that models are to be written dimensionless. This state is achieved if we make the transforma- Under these transformations, we have the following Lengyel-Epstein model [5]: Kinetic applications, in elucidation, are a controlling gizmo in the study of the reaction structure, allowing to gather vital facts of the processes that transpire before the influential phase of the speed. Through kinetic investigation one can decide the speed law of a reaction just as its steady rate. One of the methodologies for this is the utilization of integrated equations. In this methodology, one checks whether the adjustment in the concentration of one of the reactants or items follows first or second order kinetics or, all the more once in a while, kinetics with higher orders or even zero order. In a reaction, reactants go through a progress state zone along the reaction organize among products and reactants, where chemical bonds are broken and changed. This change state was first projected by Eyring and Polanyi in the mid-1930s. Numerous scientists controlled pretty much every part of a chemical reaction, and have assumed a focal visionary aspect in the advancement of chemistry as a part of science. Along these lines, watching and understanding the progress state have been viewed as the "Holy Grail" of chemistry [6]. Model (6) is a nonfractional order system, i.e., it contains the first order derivative with respect to time variable τ . The first-order derivative with respect to the variable τ infers the transient change pace of these chemical reactions. Though, as a result of the intricacy of reactions which were biochemical in nature, chemical reaction practices are regularly influenced by or rely on the historical background of chemical reactions. Many researchers have solved the Lengyel-Epstein chemical model [7][8][9][10]. A number of researchers were using fractional order techniques and they claimed these to be more suitable than classical ones [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26].
In this paper, we are keen on an old style version of the Lengyel-Epstein reaction diffusion structure as a specimen of the chlorite-iodide-malonic-acid (CIMA) reaction. The deliberated model has been pulled into the light of legitimate concern for some analysts since its origin in 1991. The purpose behind this intrigue is the way that CIMA reaction is perhaps the most punctual trial that restricted the theoretical suggestions of Alan Turing in 1952, concerning the chemical basis for morphogenesis and, all the more by and large, pattern formation [27].
This paper is systematized into five sections. The introduction is the first section in which we intricate some history of the Lengyel-Epstein model in kinetic studies. In Sect. 2, we will study the Lengyel-Epstein model with diffusion. Besides this, we will study the steady states with and without diffusion of the Lengyel-Epstein model. In Sect. 3, we ponder the global steady states for the system. In Sects. 4 and 5, numerical finding are offered for two systems using three well known techniques to validate the main outcomes. Results and Discussion are presented in Sect. 6, and conclusions are drawn in Sect. 7.

Model with diffusion
Here, we consider one-dimensional coupled Lengyel-Epstein model: where Υ is a bounded domain in R N with sufficiently smooth boundary ∂Υ . Here, u = u(x, τ ) and v = v(x, τ ) denote the concentration of the inhibitor chlorite (ClO 2 ) and the activator iodide (I -), respectively, at time τ > 0 and point x ∈ Υ . The constants l and m are restrictions depending on the concentration of the starch, broadening the diffusion ratio to be effective by κ 1 and κ 2 . The constants κ 1 , κ 2 , l and m are nonnegative. The function ψ is supposed to be positive and continuously differentiable on R + such that ψ(0) = 0, and for u ∈ (0, l),

Invariant regions
Now, we study the invariant zones for the system (7). , a = 0.

Steady states of model (6) without diffusion
We shall consider system (6) without diffusion.
Proof An equilibrium node (u * , v * ) of (6) solves the system Simple calculations from (12) provide the solution given below: Next, we ponder the stability of the equilibrium. The Jacobian matrix for the structure (6) at the equilibrium point (u * , v * ) is given by We have and Thus, the Jacobian of (14) has eigenvalues with negative real parts. Hence, the equilibrium node of (14) is asymptotically stable, which is bound by condition (13). Hence, the proof of Proposition 2 is completed.
is satisfied then we call u(x, t) an activator, v(x, t) an inhibitor, and the structure (6) is an activator-inhibitor structure. Since ψ(μ) is always positive, multiplying both sides of the inequity by ψ(μ), we get Comment 2 Merging the activator-inhibitor condition (18) with the steady state condition given in Proposition 2, we discover that the condition makes the system (6) a diffusion free stable activator-inhibitor structure.

Steady states of model (6) with diffusion
In this subsection, we shall debate the basic properties of the nonhomogeneous steady states of the Lengyel-Epstein structure. The steady state points satisfy the following system: with the homogeneous Neumann boundary conditions

Lemma 1 Conditions
Hence, ψ(u) u is a decreasing function. Now, for some s ∈ (0, u), we have This yields .
Proof If at some node inῩ the function u reaches its extreme overῩ , then by (19) and Proposition 3, at this point we have It implies that u < l. Likewise, if v achieves an extreme overῩ at some node, then by (19) and Proposition 3, . Since (11) warranties that u(ψ(u)) -1 is growing for u between 0 and ε 1 . If u achieves its minimum overῩ at some point, then This, along with condition (20), yields Therefore, ε 2 < u. Also, if v has a minimum overῩ at some point, then vψ(u) ≥ u. Therefore, we get u ψ(u) ≤ v which leads to ε 2 ψ(ε 2 ) ≤ v.

Global asymptotic stability
In this section, we study the global asymptotic stability of the system (7). The reason is to find some adequate conditions for the global stability of the steady state equations. First, assume that h a (u) = (ψ(u)) -1 (lu), which results in h a (u * ) = μ(ψ(u)) -1 . Now, system (7) can be rewritten as , then there exist a real θ between μ and u and a nonnegative function Proof Letting s = u μ and w(s) = sμ ψ(μs) , we have Using the mean value theorem, for suitable δ 1 we have Also where μs = u. Now, suppose that , then (30) reads w (δ 1 ) = μχ(θ ), which guarantees that χ(θ ) ≥ 0.

Numerical schemes
In this section, we look at two specific examples and use numerical analysis based on forward Euler, Crank-Nicolson, and nonstandard methods to investigate the solution of the system and its stability.

The CIMA model
In this first example, we suppose ψ(u) = u(1 + u 2 ) -1 , then the system (7) takes the form To discretize system (20) using the finite difference method, firstly [0, ] 2 × [0, T] is partitioned into M 2 × N parts with spatial and temporal step sizes dx = M and dτ = T N . Then, the grid points are x i = i dx and τ j = j dτ , where i = 0, 1, 2, . . . , M and j = 0, 1, 2, . . . , N . We represent u j i and v j i as the finite difference approximations of u(i dx, j dτ ) and v(i dx, j dτ ), respectively. First order temporal derivative and second order spatial derivative finite difference formulas are: (33)

Forward Euler method
Substituting the values of u τ | After simple calculation, we have where

Crank-Nicolson method
Substituting the values of u τ | After simple calculation, we have the following scheme for system (32): This scheme is unconditionally stable for the system (32).

Nonstandard finite difference method
Here, we will design a nonstandard finite difference method for the system (20). Substi- After simple calculation, we obtain the following explicit nonstandard finite difference scheme: where R 1 = κ 1 dτ (dx) 2 and R 2 = κ 2 dτ (dx) 2 . This technique is completely stable for the structure (32).

Stability analysis of nonstandard finite difference method
To find the stability bounds, the von Neumann stability technique is used. After linearizing equation (33) and then using the von Neumann stability bounds, we get dτ e iαx T(t + t), Similarly, following the same lines, (33) gives From (40) and (41), it is clear that the proposed nonstandard finite difference technique for system (32) is completely stable.

Consistency of nonstandard finite difference method
To check the uniformity of nonstandard finite difference technique, Taylor series is used. The series of u| j+1 i , u| j i+1 and u| j i-1 are given below: Consider the unconditional proposed nonstandard finite difference scheme (43) Substituting the values of u| j+1 i , u| j i+1 and u| j i-1 into the above equation and after simplification, we have Putting dτ = (dx) 3 and letting dx → 0, the above equation gives Similarly, Taylor series expansion of v| j+1 i , v| j i+1 and v| j i-1 are given below: Substituting the values of v| Putting dτ = (dx) 3 and letting dx → 0, the above equation gives (48)

Second model
In this second example, we assume ψ(u) = u(1 + e u ) -1 , then the system (7) takes the form The equilibrium point of this system is (μ, 1 + e μ ). To discretize the system (45)

Proposed nonstandard finite difference method
Here, we will present a scheme the nonstandard finite difference method for the system (45). Substituting the values of After simple calculation, we have where R 1 = κ 1 dτ (dx) 2 and R 2 = κ 2 dτ (dx) 2 . This technique is completely stable for the structure (53).

Stability analysis of nonstandard finite difference method
To find the stability bounds, the von Neumann stability technique is used. After linearizing equation (31) and then using the von Neumann stability bounds, we get

Consistency of nonstandard finite difference method
To check the uniformity of the nonstandard finite difference technique, Taylor series is used. The series of u| j+1 i , u| j i+1 and u| j i-1 are given below: Consider the unconditional proposed nonstandard finite difference scheme Putting dτ = (dx) 3 and letting dx → 0, the above equation gives Similarly, Taylor series expansion of v| j+1 i , v| j i+1 and v| j i-1 are given below: Substituting the values of v| j+1 i , v| j i+1 and v| j i-1 into equation (32) and after simplification, we get Putting dτ = (dx) 3 and letting dx → 0, the above equation gives v τ = κ 2 v + m(u -uv 1+e u ).
6 Results and discussion 6.1 Test Problem 1 Here, the proposed (i) forward Euler scheme, (ii) Crank-Nicolson scheme, and (iii) nonstandard finite difference scheme are tested on the model considered for the onedimensional CIMA problem Numerical simulations are carried out to confirm the efficiency and effectiveness of the nonstandard finite difference method. Figures 1(a)

Test Problem 2
Here, the proposed (i) forward Euler scheme, (ii) Crank-Nicolson scheme, and (iii) nonstandard finite difference scheme are tested on a model described by the one-dimensional problem as with u 0 = 1 + cos(x) and v 0 = 2 + sin(x).

Conclusion
In this article, we consider a nonlinear model with diffusion to review the dynamics of Lengyel-Epstein reaction model describing oscillating chemical reactions. The model demonstrates the connection between malonic acid, iodine, and chlorine dioxide. When simulating the model with the given methodology, we have recognized that the techniques are capturing to the true equilibrium nodes. The stability and uniformity of suggested techniques are also confirmed with the aid of von Neumann stability bounds and Taylor series, respectively. Besides, suggested nonstandard finite difference technique is unconditionally consistent regarding positivity property. The diagrams show that the suggested finite difference method is dynamically consistent with the performance of continuous systems. Conversely, two other well-known methods failed to preserve the positivity property and shown divergence on diverse step size dτ .
In the future, we will try to solve this system by using stochastic theory or by adding a delaying factor and then solving it using numerical techniques.