Convergence Analysis and Approximate Optimal Temporal Step Sizes for Some Finite Difference Methods Discretising Fisher's Equation

1. Introduction

The most enthralling recent progress of nonlinear science in particular mathematical science of partial differential equations, theoretical physics, chemistry, and engineering sciences has been a growth of strategy or procedure to try to find exact solutions for nonlinear differential equations. This is substantial due to the fact that countless mathematical models are described by nonlinear differential equations. To mention few among others the inverse scattering transform [1], the singular manifold method [2], the transformation method [3], the tanh-function method [4], and the Weierstrass function method [5] are subservient in many applications and known as stunning techniques to look for solutions of exactly solvable nonlinear partial differential equations.

In Kudryashov [6] developed a new numerical method for the solution of nonlinear partial differential equations. Amajor novelty of that technique is the utilization of finite Fourier series for the numerical approximation of the spatial derivative terms of the equations. It was proved that the precision of this method is that the order of accuracy is greater than that found by approximating the spatial derivative terms by finite-difference methods [6]. The numerical performance of this method, which was called the Accurate Space Derivatives (ASD) approach, can be run accurately by the utilization of the Fast Fourier Transform (FFT) algorithm [7]. In Kudryashov [6], the ASD approach was used to obtain the solution of nonlinear hyperbolic equations depicting convective fluid flows. Furthermore, the ASD approach was used in Gazdag and Canosa [8] to solve Fisher's equation, a nonlinear diffusion equation portraying the rate of advance of a new advantageous gene within a population of constant density inhabiting a one dimensional habitat [9]. An outstanding and compendious debate of Fisher's equation in the framework of the genetic problem can be found in Moran [10] and Kendall [11].Kolmogorov [12] used the traveling waves with wave speed c to solve the problem. They showed that by assuming that when the initial condition belongs to the interval [0, 1], the speed of propagation c of the waves is superior to two and the solution is in the form u(x, t) = w(ξ) where ξ = x − ct. They further demonstrated that there are no solutions for c ∈ [0, 1).

Fisher's equation has a limitless number of traveling wave solutions and each wave propagates at a characteristic speed, c > 2. This result appears to point out that the velocity of gene advance is undefined. Gazdag and Canosa [8] studied the problem of the indeterminacy of the diffusion speed of Fisher's equation which has not been plainly investigated in Kolmogorov [12]. Furthermore, a modification was made by Fisher [9] to his original model, he demonstrated that the rate of gene advance became the minimum one when c = 2. Kendall [11] investigated a linear model portraying a population that undergoes a Brownian motion and spreads geometrically at the same time. Canosa [13] demonstrated that all waves are stable against local perturbations but are linearly unstable against general perturbations of limitless magnitude. It is worthy to emphasize that the traveling wave profiles of Fisher's equation are similar to some of the steady-state solutions of the Korteweg-de Vries-Burgers equation which is a third-order nonlinear partial differential equation integrating diffusive and dispersive effects which have been found useful to represent blood flow through an artery, shallow water waves and plasma shocks disseminating perpendicularly to a magnetic field [13, 14].

Kudryashov [6] showed that a simple stability analysis enables us to see the estimation which is unstable against the roundoff errors growing up at the right tail of the waves. This is due to the physical nature of the problem depicted by the equation, not to the numerical method utilized and moreover entailed an exponential growth of the solution when roundoff errors are exponentially small. This simple issue makes it hard to do a strict simulation of the solutions of Fisher's equation. Kudryashov [6] went on with the removal of the forward tail of the wave of advance. This removal is necessary for the numerical stability of the Accurate Space Derivatives (ASD) approach and is physically conclusive because it is approximately equivalent to assuming that the role of long-distance dissipation in the spread of the gene is insignificant and probably effective for some species but not for others. Other numerical computations present how fast the asymptotic minimum speed wave is reached from an initial step function and confirm the local stability analysis of Kudryashov [6] which unveils that local perturbations are flattened very rapidly, even from superspeed waves. Another amazing result of the estimation is obtained for an initial dispensation localized in space which further gives rise to two identical waves of minimum speed evolution cases, one disseminating to the right and the other to the left.

1.1. Some Generalized and Conserved Fisher's Equation

Fisher's equation can be represented as generalized or conserved forms. Fisher's equation is the elementary model of spatial dynamics, in which competitive interactions between individuals happen locally. In Kudryashov and Zakharchenko [15], the generalized form is written as

$\begin{array}{ll}{u}_t=u_x(u^l{u}_x)+u^a(1-u^b)& (1)\end{array}$

where t stands for time, x stands for spatial coordinate, u is a population density, and a, b, and l are all positive parameters.

Fisher's equation can also forecast circumstances where population regulation happens globally due to the existence of a secondary agent (the controller agent is itself dispersed over a scale significantly greater than the dispersal distance of the individuals themselves). It is, thus, of notice to envisage a simple model of spatial population dynamics in which the total population size is controlled via a nonlocal mechanism. In that case in Newman et al. [16], the conserved equation is written as

$\begin{array}{ll}{u}_t=D{\nabla }^{2}u+r(t)u(t)-K(t)u{(t)}^2& (2)\end{array}$

where D stands for the mobility of the individuals, r stands for the reproduction rate in the absence of competition. K is a parameter representing the carrying capacity of the system and regulating the population density through competition. The auxiliary equation related to Equation (2) is

$\begin{array}{ll}r(t)=K(t){\displaystyle \int d^{d}xu(x,t)}& (3)\end{array}$

It is worth mentioning that Fisher's equation belongs to the class of partial differential equations called Reaction-Diffusion equations. This class of equations has broad applications in science and engineering for instance transport of air, adsorption of pollutants in soil, food processing, and modeling of biological and ecological systems [17, 18]. Several reaction-diffusion equations involve traveling waves fronts yielding a fundamental role in the understanding of physical, chemical, and biological phenomena [19]. Reaction-diffusion systems clarify how the condensation of one or more substances diffused in space varies by the impact of two operations: first, it is local chemical reactions in which the substances are modified into each other, and second, it is the diffusion that sustains the substances to smear over a surface in space [20]. Reaction-diffusion systems are regularly used in chemistry. Nonetheless, the system can also portray the dynamical processes of non-chemical nature. Reaction-diffusion systems have mathematically the form of semi-linear parabolic partial differential equations. They are often written in the form of

$\begin{array}{ll}{u}_t=D{\nabla }^{2}u+R(u),& (4)\end{array}$

where each component of vector u(x, t) stands for the concentration of one substance, x is the space variable, and t is the time. D is the diffusion coefficient and R represents the reaction term. The solution of reaction-diffusion equations shows an ample scale of behaviors, enclosing the formation of traveling waves. These waves are like phenomena and self-organized patterns which are e.g., stripes, hexagons, or more complicated fabrics such as dissipative solutions [20]. Reaction-diffusion equations are also grouped as one component, two components, or more component diffusion equations counting upon the component involved in the reaction. The basic reaction-diffusion equation regarding the concentration u of a mere substance in one spatial dimension is

$\begin{array}{ll}{u}_t=D{u}_{xx}+R(u).& (5)\end{array}$

If the reaction R(u) term goes away, then the equation gives a pure diffusion process and if the thermal diffusivity term appears instead of diffusion term D then the equation will turn into a parabolic partial differential equation in one dimensional space [20]. The choice R(u) = u(1 − u) yields Fisher's equation. Those reaction-diffusion equations arise also in flame propagation, the branching Brownian motion process, and nuclear reactor theory [20]. Many methods such as Adomian Decomposition [21], Variational Iteration [22], Factorization [23], Nonstandard Finite Difference, and Exponential Finite Difference methods [24, 25] are used to solve Fisher's equation.

Anguelov et al. [26] investigated the same Fisher's equation by the means of a periodic initial data with θ-non standard approach and found that the Nonstandard Finite Difference approach is elementary stable in the limit case of space independent variable, stable in regard to the boundedness and positivity property. Finally also stable in regard to the conservation of energy in the stationary case.

Let us consider simple Fisher's equation given by

$\begin{array}{ll}{u}_t=u_{xx}+R(u),& (6)\end{array}$

where R(u) = u(1 − u) and x ∈ ℝ, t positive. The boundary and initial conditions are as follows

$\begin{array}{ll}{\displaystyle \underset{x\to +\left|ϵ\right|\infty }{lim}}u(x,t)=\{\begin{array}{l}1\text{ if }ϵ=1\\ 0\text{ if }ϵ=-1,\end{array}& (7)\end{array}$

$\begin{array}{ll}u(x,0)=u_0(x).& (8)\end{array}$

Hagstrom and Keller [27] revealed that when a positive function is taken as an initial condition satisfying

$\begin{array}{ll}{u}_0(x)~exp(-\alpha )\text{ when }x\to \infty ,& (9)\end{array}$

then the solution u develops a traveling wave speed in function of α which is

$\begin{array}{ll}c(\alpha )=\{\begin{array}{c}\alpha +\frac{1}{\alpha },\alpha \le 1,\\ 2,\alpha \ge 1.\end{array}& (10)\end{array}$

2. Organization of the Article

The organization of this article is as follows. In Section 3, we present the general form of the exact solution of Fisher's equation and in Section 4, we describe the numerical experiment [28]. In Section 5, we make use of Forward in Time Central Space (FTCS) in order to discretize Fisher's equation, study the stability and consistency and we also obtain error estimates. Sections 6, 7 discuss stability, consistency, and error estimates for NSFD and EEFD schemes. In Section 8, we conclude by presenting the important highlights of this article. The computations are performed by making use of MATLAB R2014a software on an intel core2 as CPU.

3. Exact Solution

In this section, we present the exact solution of generalized Fisher's equations as described in Kudryashov and Zakharchenko [15]. The nonlinear evolutional equation of that generalized Fisher's equation gives one dimensional diffusion models (for insect, animal dispersal, and invasion) as

$\begin{array}{ll}{u}_t=u_x(u^l{u}_x)+u^a(1-u^b),& (11)\end{array}$

where t stands for time, x stands for spatial coordinate, u is population density, and a, b, and l are positive parameters. The first term, $u_x(u^l{u}_x)$ on the right-hand side of Equation (11) stands for the growth of population. The term u^l represents the diffusion process depending on the population density.

Let us consider l ≠ 0 and u(x, t) = v(ξ) and ξ = sx − ct, s ≠ 0. Equation (11) gives the following nonlinear ordinary differential equation

$\begin{array}{ll}{s}^2\frac{d}{d\xi }\left(v^l\frac{dv}{d\xi }\right)+v^a-v^{a+b}+c\frac{dv}{d\xi }=0.& (12)\end{array}$

For l ≠ 0, v^l = w. Replacing v by $w^{\frac{1}{l}}$ in Equation (12) gives

$\begin{array}{ll}\frac{s^2}{l}{w}_{\xi }^2+s^{2}w{w}_{\xi \xi }-l{w}^{\frac{a+b+l-1}{l}}+w^{\frac{a+l-1}{l}}+\xi w_{\xi }=0.& (13)\end{array}$

Using the Q function method as it is in Kudryashov [29], one has

$\begin{array}{ll}w(\xi )={\displaystyle \sum _{j=0}^P}{P}_j{Q}^j(\xi ),Q(\xi )=\frac{1}{1+e^{\xi -{\xi }_0}}& (14)\end{array}$

where P stands for the pole order and ξ₀ stands as an arbitrary constant. Q(ξ) is the solution of

$\begin{array}{ll}{Q}_{\xi }=Q-Q^2.& (15)\end{array}$

Using Equation (15), we obtain w_ξ and w_ξξ by using polynomials of Q. Replacing w ≃ Q^P into Equation (15), we have for b ≥ 0, the following equality

$\begin{array}{ll}\frac{a+b+l-1}{l}=2+\frac{2}{P}.& (16)\end{array}$

To have an integer value $\frac{a+b+l-1}{l},$ P should be 1 or 2. In that case

$\begin{array}{ll}w(\xi )=\{\begin{array}{l}{P}_0+P_{1}Q(\xi )\\ P_0+P_{1}Q(\xi )+P_2{Q}^2(\xi )\end{array}& (17)\end{array}$

For P = 1, a = 3l + 1 − b, Equation (13) becomes

$\begin{array}{ll}\frac{s^2}{l}{w}_{\xi }^2+s^{2}w{w}_{\xi \xi }-l{w}^4+w^{4-\frac{b}{l}}+\xi w_{\xi }=0.& (18)\end{array}$

We have the following solutions

$\begin{array}{ll}w(x,t)=\{\begin{array}{l}\text{No exact solutions},b=l\text{ and }b=3l,\\ \pm 1\pm 2Q(\pm \frac{2l}{\sqrt{2l+1}}x\pm \frac{2l}{2l+1}t),b=2l,\\ \pm 1\pm 2Q(\pm \frac{2l}{\sqrt{2l+1}}x\pm \frac{4l(l+1)}{2l+1}t),b=4l,\\ \text{ or}\\ \pm i\pm 2iQ(\pm \frac{2l}{\sqrt{2l+1}}x\pm \frac{4l(l+1)}{2l+1}it),b=4l.\end{array}& (19)\end{array}$

We finally get as exact solution

$\begin{array}{ll}u\left(x,t\right)=\{\begin{array}{l}\sqrt[l]{\pm 1\pm 2Q\left(\pm \frac{2l}{\sqrt{2l+1}}x\pm \frac{2l}{2l+1}t\right)},V_1=\pm \frac{1}{\sqrt{2l+1}},\\ \sqrt[l]{\pm 1\pm 2Q\left(\pm \frac{2l}{\sqrt{2l+1}}x\pm \frac{4l\left(l+1\right)}{2l+1}t\right)},V_2=\pm \frac{2\left(l+1\right)}{\sqrt{2l+1}},\\ \text{or}\\ \sqrt[l]{\pm i\pm 2iQ\left(\pm \frac{2l}{\sqrt{2l+1}}x\pm \frac{4l\left(l+1\right)}{2l+1}it\right)},V_2=\pm \frac{2\left(l+1\right)}{\sqrt{2l+1}},\end{array}& (20)\end{array}$

where V₁ and V₂ stand for velocity.

For P = 2, a = 2l + 1 − b. With the same reasoning as above, we have

$\begin{array}{ll}u\left(x,t\right)=\{\begin{array}{l}\text{No exact solutions},b=2l\text{ and }b=3l,\\ \sqrt[l]{\frac{2\left(3l+2\right)}{l+1}\left(Q\left(\pm \frac{l}{\sqrt{l+1}}ix\right)-Q^2\left(\pm \frac{l}{\sqrt{l+1}}ix\right)\right)},b=l.\end{array}& (21)\end{array}$

The case of l = 0, a = 1 and b = 2, one obtains the Burgers-Huxley equation and the case of l = 0, a = 1, b = 1, we have Fisher's equation. The exact solution of Equation (11) is described in Li et al. [28] as a scaled Fisher's equation in the form

$\begin{array}{ll}{u}_t=u_{xx}+\rho u\left(1-u\right),& (22)\end{array}$

with x ∈ ℝ, t positive, and ρ is a positive constant. The Equations (7) and (8) stand for boundary and initial conditions, respectively. The traveling exact solution to this problem as presented in Polyanim and Zaitsev [30] is

$\begin{array}{ll}u\left(x,t\right)={\left[1+cexp\left(\sqrt{\frac{\rho }{6}}x-\frac{5\rho }{6}t\right)\right]}^{-2},& (23)\end{array}$

where $c=5\sqrt{\rho /6}$ stands for the wave speed with the minimum value, $2\sqrt{\rho }.$

4. Numerical Experiments

We consider the following problem from Qiu and Sloan [31].

Solve

$\begin{array}{l}{u}_t=u_{xx}+\rho u\left(1-u\right),\end{array}$

for x ∈ [−0.2, 0.8] and t ∈ [0, T_max] where $T_{max}=2.5\times 10^{-3}.$ The initial data is given by

$\begin{array}{ll}u\left(x,0\right)={\left[1+exp\left(\sqrt{\frac{\rho }{6}}x\right)\right]}^{-2}.& (24)\end{array}$

The exact solution is given by

$\begin{array}{ll}u\left(x,t\right)={\left[1+exp\left(\sqrt{\frac{\rho }{6}}x-\frac{5\rho }{6}t\right)\right]}^{-2}.& (25)\end{array}$

The boundary conditions are as follows

$\begin{array}{l}u\left(-0.2,t\right)={\left[1+exp\left(-0.2\sqrt{\frac{\rho }{6}}-\frac{5\rho }{6}t\right)\right]}^{-2}\text{and}\\ u\left(0.8,t\right)={\left[1+exp\left(0.8\sqrt{\frac{\rho }{6}}-\frac{5\rho }{6}t\right)\right]}^{-2}.\end{array}$

We consider three cases ρ namely; 1, 10², 10⁴, and obtain a solution at time, t = T_max.

DEFINITION 1. Miyata and Sakai [32]. For a vector $\overline{x}\in {ℝ}^N,\parallel \overline{x}{\parallel }_1=\sum _{i=1}^N\left|{\overline{x}}^i\right|$ and $\parallel \overline{x}{\parallel }_{\infty }=max\left\{\left|{\overline{x}}^i\right|,i=1,···,N\right\}.$

DEFINITION 2. Sutton [33]. Suppose ${\left\{t^n\right\}}_0^N$ forms a partition of [0, T], with t_n = nΔt for n = 0, ···, N, where Δt = T/N. Suppose a vector $\overline{x}\in {ℝ}^N,$ defined by

$\begin{array}{ll}{\parallel \overline{x}\parallel }_{L^P\left(0,t^n\right)}=\{\begin{array}{l}{\left({\parallel \overline{x}\parallel }_{L^P\left(0,t^{n-1}\right)}+\tau {\left({\overline{x}}^n\right)}^p\right)}^{\frac{1}{p}}forp\in \left[0,\infty \right),\\ max\left\{{\parallel \overline{x}\parallel }_{L^P\left(0,t^{n-1}\right)},{\overline{x}}^n\right\}forp=\infty .\end{array}& (26)\end{array}$

The rate of convergence with respect to time is defined by

$rat{e}_i\left(t\right)=\frac{log\left({\overline{x}}^i\left(t\right)\right)-log\left({\overline{x}}^{i-1}\left(t\right)\right)}{log\left(\text{Δ}{t}^i\right)-log\left(\text{Δ}{t}^{i-1}\right)}.$

5. Forward in Time Central Space

The discretization of Equation (22) using the FTCS method gives [34].

$\begin{array}{ll}\frac{u_m^{n+1}-u_m^n}{k}=\frac{u_{m+1}^n-2u_m^n+u_{m-1}^n}{h^2}+\rho u_m^n\left(1-u_m^n\right),& (27)\end{array}$

which leads to

$\begin{array}{ll}{u}_m^{n+1}=\left(1-2R\right)u_m^n+k\rho u_m^n\left(1-u_m^n\right)+R\left(u_{m+1}^n+u_{m-1}^n\right),& (28)\end{array}$

or

$\begin{array}{ll}{u}_m^{n+1}=\left(1-2R+k\rho \right)u_m^n-k\rho {\left(u_m^n\right)}^2+R{u}_{m+1}^n+R{u}_{m-1}^n,& (29)\end{array}$

where $R=\frac{k}{h^2}.$

5.1. Stability

The investigation regarding the stability of the scheme given by Equation (28) was done in Agbavon et al. [35]. Nevertheless, we highlight briefly some points. The stability of finite difference methods discretizing nonlinear partial differential equations is not straightforward. Subsequently, freezing the coefficients is needed before using Von Neumann stability analysis [36].

THEOREM 1. Agbavon et al. [35]. The FTCS scheme given by Equation (28) is conditionally stable, and the stability region is

$\begin{array}{ll}k\le \frac{h^2}{2}& (30)\end{array}$

for given spatial step size h > 0 and the time step size k > 0. FTCS is first order and second order accurate in time and space, respectively.

Proof. The stability region of the Zabusky and Kruskal scheme using the Korteweg de Vries (KdV) equation was found by Taha and Ablowitz [37] by using the freezing coefficients method and Von Neumann Stability Analysis. The derived scheme by Zabusky and Kruskal [38] for the KdV equation, u_t + 6uu_x + u_xxx = 0 is

$\begin{array}{l}\frac{u_m^{n+1}-u_m^{n-1}}{2k}+6\left(\frac{u_{m+1}^n+u_m^n+u_{m-1}^n}{3}\right)\left(\frac{u_{m+1}^n-u_{m-1}^n}{2h}\right)\\ +\frac{1}{2h^3}\left(u_{m+2}^n-2u_{m+1}^n+2u_{m-1}^n-u_{m-2}^n\right)=0.& (31)\end{array}$

Taha and Ablowitz [37] express u u_x as u_max u_x and utilize the ansatz

$u_m^n={\xi }^n{e}^{Imw}$ where w stands for the phase angle. They obtain

$\begin{array}{l}(\xi -{\xi }^{-1}){(2k)}^{-1}+{(h)}^{-1}(6u_{max})\text{I}sin(w)+{(2h^3)}^{-1}(e^{2Iw}\\ -2e^{Iw}+2e^{-Iw}-e^{-2Iw})=0,\end{array}$

which can be rewritten as

$\begin{array}{l}\xi ={\xi }^{-1}-{(h)}^{-1}(12k{u}_{max})\text{I}sin(w)-{(h^3)}^{-1}k(e^{2Iw}\\ -2e^{Iw}+2e^{-Iw}-e^{-2Iw})& (32)\end{array}$

where u_max = max|u(x, t)|. The requirement for the linear stability is

$\begin{array}{ll}{(h)}^{-1}k|2u_{max}-{(h^2)}^{-1}|\le 2{(3\sqrt{3})}^{-1}.& (33)\end{array}$

For obtaining the stability region of the FTCS scheme discretizing Equation (28), we rewrite Equation (28) using the same idea as

$\begin{array}{l}{u}_m^{n+1}=(1-{(h^2)}^{-1}2k)u_m^n+{(h^2)}^{-1}(k)(u_{m+1}^n+u_{m-1}^n)\\ +k\rho u_m^n-k\rho {(u_m^n)}^2.& (34)\end{array}$

Utilization of Fourier series analysis on Equation (34), gives the amplification factor

$\begin{array}{ll}\xi =1-{(h^2)}^{-1}(2k)(1-cos(w))+k\rho (1-u_{max}),& (35)\end{array}$

where u_max is frozen coefficient. For the numerical experiment considered, we have u_max = 1, and therefore,

$\begin{array}{ll}\xi =1-{\left(h^2\right)}^{-1}\left(4k\right){sin}^2\left(\frac{w}{2}\right).& (36)\end{array}$

The stability is obtained by solving |ξ| ≤ 1 for w ∈ [−π, π], and we obtain $k\le \frac{h^2}{2}.$

Using Taylor series expansion about the point (n, m) of Equation (28), we get

$\begin{array}{l}u+k{u}_t+\frac{k^2}{2}{u}_{tt}+\frac{k^3}{6}{u}_{ttt}+O\left(k^4\right)\\ =\left(1-{\left(h^2\right)}^{-1}\left(2k\right)+k\rho \right)u-k\rho u^2\\ +{\left(h^2\right)}^{-1}\left(k\right)\left(u+h{u}_x+\frac{h^2}{2}{u}_{xx}+\frac{h^3}{6}{u}_{xxx}+\frac{h^4}{24}{u}_{xxxx}+O\left(h^5\right)\right)\\ +{\left(h^2\right)}^{-1}\left(k\right)\left(u-h{u}_x+\frac{h^2}{2}{u}_{xx}-\frac{h^3}{6}{u}_{xxx}+\frac{h^4}{24}{u}_{xxxx}+O\left(h^5\right)\right),& (37)\end{array}$

which can be written as

$\begin{array}{l}{u}_t-u_{xx}-\rho u\left(1-u\right)=-\frac{k}{2}{u}_{tt}-\frac{k^2}{6}{u}_{ttt}\\ +\frac{h^2}{12}{u}_{xxxx}+O\left(k^4\right)+O\left(h^5\right).& (38)\end{array}$

Hence, FTCS is first order and second order accurate in time and space, respectively.

5.2. Error Estimates

THEOREM 2. Let u ∈ C^{4, 2}(Q), Q defined by Q = {(x, t)/a ≤ x ≤ b, 0 < t ≤ T, a, b ∈ ℝ}. If spatial step size, h and time step size, k are such that the stability condition (30) holds, then the error estimate, $E_m^n$ for Equation (28) is given by

$\begin{array}{ll}{E}_m^n\le {\left(1+3k\rho \right)}^n{E}_m^0+\frac{1}{9}\frac{M{h}^2}{\rho k}\left[{\left(1+3k\rho \right)}^n-1\right]& (39)\end{array}$

where $M=\underset{\left\{(x,t)\right\}\in Q}{max}\left\{\left|u_{xxxx}(x,t)\right|,\left|u_{tt}(x,t)\right|\right\}$ and Θ such that Θ(k, h)u_ttt = O(k, h) → 0, for k, h → 0.

Proof.

Forward in Time Central Space scheme is given by

$\begin{array}{ll}{u}_m^{n+1}=\left(1-2R+k\rho \right)u_m^n-k\rho {\left(u_m^n\right)}^2+R{u}_{m+1}^n+R{u}_{m-1}^n.& (40)\end{array}$

Taylor series expansion about (n, m) gives

$\begin{array}{l}7v+k{v}_t+\frac{k^2}{2}{v}_{tt}+\frac{k^3}{6}{v}_{ttt}+O\left(k^4\right)\\ =\left(1-{\left(h^2\right)}^{-1}\left(2k\right)+k\rho \right)v-k\rho v^2\\ +{\left(h^2\right)}^{-1}\left(k\right)\left(v+h{v}_x+\frac{h^2}{2}{v}_{xx}+\frac{h^3}{6}{v}_{xxx}+\frac{h^4}{24}{v}_{xxxx}+O\left(h^5\right)\right)\\ +{\left(h^2\right)}^{-1}\left(k\right)\left(v-h{v}_x+\frac{h^2}{2}{v}_{xx}-\frac{h^3}{6}{v}_{xxx}+\frac{h^4}{24}{v}_{xxxx}+O\left(h^5\right)\right),& (41)\end{array}$

which can be rewritten as

$\begin{array}{ll}{v}_t-v_{xx}-\rho v\left(1-v\right)=-\frac{k}{2}{v}_{tt}-\frac{k^2}{6}{v}_{ttt}+\frac{h^2}{12}{v}_{xxxx}+...,& (42)\end{array}$

and let $\text{Θ}(k,h)v_{ttt}=O(k,h)=-\frac{k^2}{6}{v}_{ttt}\to 0$ for k, h → 0. The exact discrete equation is

$\begin{array}{l}{u}_m^{n+1}=\left(1-2R\right)u_m^n+k\rho u_m^n\left(1-u_m^n\right)+R\left(u_{m+1}^n+u_{m-1}^n\right)\\ +\frac{k}{2}{u}_{tt}\left(x,{\tau }_n\right)-\frac{h^2}{12}{u}_{xxxx}\left(X_m,t\right)& (43)\end{array}$

where x_m < X_m < x_{m + 1} and t_n < τ_n < t_{n + 1}. We define

$e_m^n=u_m^n-v_m^n\Rightarrow e_m^{n+1}=u_m^{n+1}-v_m^{n+1}.$

It follows that

$\begin{array}{l}{e}_m^{n+1}=\left(1-2R\right)\left(u_m^n-v_m^n\right)+k\rho u_m^n\left(1-u_m^n\right)\\ -k\rho v_m^n\left(1-v_m^n\right)+R\left(u_{m+1}^n+u_{m-1}^n\right)\\ -R\left(v_{m+1}^n+v_{m-1}^n\right)+\frac{k}{2}{u}_{tt}\left(x,{\tau }_n\right)-\frac{h^2}{12}{u}_{xxxx}\left(X_m,t\right).& (44)\end{array}$

We have

$\begin{array}{l}{e}_m^{n+1}=\left(1-2R\right)\left(u_m^n-v_m^n\right)+k\rho u_m^n-k\rho {\left(u_m^n\right)}^2-k\rho v_m^n\\ +k\rho {\left(v_m^n\right)}^2+R\left(u_{m+1}^n-v_{m+1}^n\right)\\ -R\left(u_{m-1}^n-v_{m-1}^n\right)+\frac{k}{2}{u}_{tt}\left(x,{\tau }_n\right)-\frac{h^2}{12}{u}_{xxxx}\left(X_m,t\right),& (45)\end{array}$

which can be rewritten as

$\begin{array}{l}{e}_m^{n+1}=\left(1-2R\right)\left(u_m^n-v_m^n\right)+k\rho \left(u_m^n-v_m^n\right)\\ -k\rho \left(u_m^n-v_m^n\right)\left(u_m^n+v_m^n\right)+R\left(u_{m+1}^n-v_{m+1}^n\right)\\ -R\left(u_{m-1}^n-v_{m-1}^n\right)+\frac{k}{2}{u}_{tt}\left(x,{\tau }_n\right)-\frac{h^2}{12}{u}_{xxxx}\left(X_m,t\right).& (46)\end{array}$

Using the properties of absolute values |a + b| ≤ |a| + |b| for a, b ∈ ℝ, we have

$\begin{array}{l}\left|e_m^{n+1}\right|\le |1-2R|\left|e_m^n\right|+\left|k\rho \right|\left|e_m^n\right|+\left|k\rho \right|\left|e_m^n\right||u_m^n+v_m^n|\\ +R\left|e_{m+1}^n\right|\\ +R\left|e_{m-1}^n\right|+M\left(\frac{k}{2}+\frac{h^2}{12}\right),& (47)\end{array}$

where $M=\underset{\left\{(x,t)\in Q\right\}}{max}\left\{\left|u_{xxxx}(x,t)\right|,\left|u_{tt}(x,t)\right|\right\}.$ Since $0\le u_m^n\le 1$ and $0\le v_m^n\le 1,$ based on numerical experiment chosen, we have

$\begin{array}{l}\left|e_m^{n+1}\right|\le |1-2R|\left|e_m^n\right|\\ +k\rho \left|e_m^n\right|+2k\rho \left|e_m^n\right|+R\left|e_{m+1}^n\right|\\ +R\left|e_{m-1}^n\right|+M\left(\frac{k}{2}+\frac{h^2}{12}\right).& (48)\end{array}$

Let $E_m^n=\underset{0<m<N}{max}\{\left|e_m^n\right|)\left|\right\}.$ We have

$\begin{array}{l}\left|e_m^{n+1}\right|\le \left(|1-2R|+k\rho +2k\rho +2R\right)\left|E_m^n\right|\\ +M\left(\frac{k}{2}+\frac{h^2}{12}\right),& (49)\end{array}$

and for stability R ≤ 1/2, therefore, |1 − 2R| = 1 − 2R ≥ 0. We finally obtain

$\begin{array}{l}\left|e_m^{n+1}\right|\le \left(1+3k\rho \right)E_m^n\\ +M\left(\frac{k}{2}+\frac{h^2}{12}\right).& (50)\end{array}$

Let $E_m^{n+1}=\left(1+3k\rho \right)E_m^n+\left(\frac{k}{2}+\frac{h^2}{12}\right)M.$ We have

For n = 0, $E_m^1=\left(1+3k\rho \right)E_m^0+\left(\frac{k}{2}+\frac{h^2}{12}\right)M.$

For n = 1, we have

$\begin{array}{ll}{E}_m^2=\left(1+3k\rho \right)E_m^1+\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (51)\end{array}$

$\begin{array}{l}={\left(1+3k\rho \right)}^2{E}_m^0+{\left(1+3k\rho \right)}^1(\frac{k}{2}\\ +\frac{h^2}{12})M+{\left(1+3k\rho \right)}^0\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (52)\end{array}$

$\begin{array}{l}={\left(1+3k\rho \right)}^2{E}_m^0+[{\left(1+3k\rho \right)}^1\\ +{\left(1+3k\rho \right)}^0]\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (53)\end{array}$

For n = 2, we have

$\begin{array}{ll}{E}_m^3={\left(1+3k\rho \right)}^3{E}_m^2+\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (54)\end{array}$

$\begin{array}{l}={\left(1+3k\rho \right)}^3{E}_m^0+{\left(1+3k\rho \right)}^2(\frac{k}{2}\\ +\frac{h^2}{12})M+{\left(1+3k\rho \right)}^1\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (55)\end{array}$

$\begin{array}{l}+{\left(1+3k\rho \right)}^0\left(\frac{k}{2}+\frac{h^2}{12}\right)M\\ ={\left(1+3k\rho \right)}^3{E}_m^0+[{\left(1+3k\rho \right)}^2\\ +{\left(1+3k\rho \right)}^1+{\left(1+3k\rho \right)}^0]\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (56)\end{array}$

For n, we have

$\begin{array}{l}{E}_m^n=\left(1+3k\rho \right)E_m^{n-1}+\left(\frac{k}{2}+\frac{h^2}{12}\right)M\\ ={\left(1+3k\rho \right)}^n{E}_m^0+[{\left(1+3k\rho \right)}^{n-1}\\ +{\left(1+3k\rho \right)}^{n-2}+··+{\left(1+3k\rho \right)}^1\\ +{\left(1+3k\rho \right)}^0]\left(\frac{k}{2}+\frac{h^2}{12}\right)M\\ ={\left(1+3k\rho \right)}^n{E}_m^0+\left[{\displaystyle \sum _{i=0}^{n-1}}{\left(1+3k\rho \right)}^i\right]\left(\frac{k}{2}+\frac{h^2}{12}\right)M\\ ={\left(1+3k\rho \right)}^n{E}_m^0+\left[\frac{1-{\left(1+3k\rho \right)}^n}{1-\left(1+3k\rho \right)}\right]\left(\frac{k}{2}+\frac{h^2}{12}\right)M\\ ={\left(1+3k\rho \right)}^n{E}_m^0-\frac{1}{3k\rho }\left[1-{\left(1+3k\rho \right)}^n\right]\left(\frac{k}{2}+\frac{h^2}{12}\right)M& (57)\end{array}$

Hence, for $k\le \frac{h^2}{2},$ we can also write

$\begin{array}{ll}{E}_m^n\le {\left(1+3k\rho \right)}^n{E}_m^0+\frac{1}{9}\frac{M{h}^2}{\rho k}\left[{\left(1+3k\rho \right)}^n-1\right]& (58)\end{array}$

6. Nonstandard Finite Difference Scheme (NSFD)

Over the past decade, the NSFD has been used extensively and often abbreviated as NSFD. The method was introduced by Mickens for the approximation of solutions of partial differential equations and is largely based on the concept of dynamical consistency [39] which plays a significant role in the construction of discrete models whose numerical solution can be complicated to compute. The dynamical consistency is bound to a precise property of a physical system (and varies according to the systems). To mention few among others these properties include the preservation of positivity, boundedness, monotonicity of the solutions, and stability of fixed-points [39]. The main advantage of this method was the dismissal of the primary numerical instabilities [40] caused by the use of standard methods. In order to reduce numerical sensitivities appearing using the classical finite difference methods, these NSFD were developed.For practical use, the construction of NSFD methods is based on the following basic rules [39]:

(1) The order of discrete derivatives should be equal to the order of corresponding derivatives appearing in the differential equation.

(2) Discrete representation for derivatives, in general, have non trivial denominator functions, e.g.,

$\begin{array}{ll}{u}_t\approx \frac{u_m^{n+1}-u_m^n}{\varphi \left(\text{Δ}t,\text{λ}\right)}& (59)\end{array}$

where

$\begin{array}{ll}\varphi \left(\text{Δ}t,\text{λ}\right)=\text{Δ}t+O\left(\text{Δ}{t}^2\right).& (60)\end{array}$

6.1. Example of the Definition of the Function ϕ

Consider the following decay equation and logistic growth equation, respectively as in Anguelov et al. [26]

$\begin{array}{ll}\{\begin{array}{ll}{u}^{\prime }=\text{λ}u,& u\left(0\right)=u_0,\text{λ}\ne 0,\\ u^{\prime }=\text{λ}u\left(1-u\right),& u\left(0\right)=u_0,\text{λ}>0,\end{array}& (61)\end{array}$

and the respective solutions at the time t = t_n+1 are

$\begin{array}{ll}\{\begin{array}{l}u\left(t_{n+1}\right)=u_0{e}^{\text{λ}{t}_{n+1}},\\ u\left(t_{n+1}\right)=\frac{u_0}{e^{-\text{λ}{t}_{n+1}}+\left(1-e^{-\text{λ}{t}_{n+1}}\right)u_0}.\end{array}& (62)\end{array}$

Let $u(t_n)=u^n.$ We have

$\begin{array}{ll}\{\begin{array}{l}\frac{u^{n+1}-u^n}{\left(e^{\text{λ Δ}t}-1\right){\text{λ}}^{-1}}=\text{λ }{u}^n,\\ \frac{u^{n+1}-u^n}{\left(e^{\text{λ Δ}t}-1\right){\text{λ}}^{-1}}=\text{λ }{u}^n\left(1-u^n\right).\end{array}& (63)\end{array}$

The Equation (63) is called the exact scheme. The function ϕ can be then defined as

$\varphi \left(\text{Δ}t,\text{λ}\right)=\frac{e^{\text{λ Δ}t}-1}{\text{λ}}\text{or }\varphi \left(\text{Δ}t,\text{λ}\right)=\frac{1-e^{-\text{λ Δ}t}}{\text{λ}}$

(3) Nonlocal discrete representations of nonlinear terms. For instance

$\begin{array}{ll}{u}_m^2\approx u_m{u}_{m+1},,u_m^2\approx \left(\frac{u_{m-1}+u_m+u_{m-1}}{3}\right)u_m,& (64)\end{array}$

and

$\begin{array}{ll}{u}^3\approx 2u_m^3-u_m^2{u}_{m+1},u_m^3\approx u_{m-1}{u}_m{u}_{m+1}.& (65)\end{array}$

In Agbavon et al. [35], followed by the rule of Mickens [39] a NSFD

for Equation (22) is

$\begin{array}{l}\frac{u_m^{n+1}-u_m^n}{\varphi \left(\text{Δ}t\right)}=\frac{u_{m+1}^n-2u_m^n+u_{m-1}^n}{{\left(\text{Δ}x\right)}^2}\\ +\rho u_m^n-\rho \left(\frac{u_{m+1}^n+u_m^n+u_{m-1}^n}{3}\right)u_m^{n+1},& (66)\end{array}$

where

$\begin{array}{ll}\varphi \left(\text{Δ}t\right)=\varphi \left(k\right)=\frac{1-e^{-\text{λ}k}}{\text{λ}};{\left(\text{Δ}x\right)}^2=h^2.& (67)\end{array}$

The Equation (66) gives the following single expression

$\begin{array}{ll}{u}_m^{n+1}=\frac{\left(1-\frac{2\varphi \left(k\right)}{h^2}+\rho \varphi \left(k\right)\right)u_m^n+\frac{\varphi \left(k\right)}{h^2}\left(u_{m+1}^n+u_{m-1}^n\right)}{1+\frac{\rho \varphi \left(k\right)}{3}\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}.& (68)\end{array}$

6.2. Positivity and Boundedness

In this section, the dynamical consistency and some useful relationship between time and space step-sizes of NSFD are presented.

THEOREM 3. The dynamical consistency (positivity and boundedness) of NSFD constructed in Equation (68) holds for Equation (22) and for relevant time step k, spatial step h if the following conditions hold

(a) $\varphi (k)\le \frac{h^2}{2-\rho h^2}\left[1-\text{Γ}\right]$ with Γ = 1 − 2^(h2)−1ϕ(k) + ρϕ(k),

(b) For $u_m^i\in \left[0,1\right]$, ∀ i. $\text{Γ}={(h^2)}^{-1}(\varphi (k))=\frac{1}{2}\left[\frac{1}{1-\frac{\rho h^2}{2}}\right]$ and Γ′ = σ Γ.

Proof.

We assume u(x, 0) = h(x) ∈ [0, 1]. We have, therefore, u(x, t) ∈ [0, 1] [24]. We assume also $u_m^n\ge 0.$ The quantity $u_m^{n+1}$ from Equation (68) is positive ($u_m^{n+1}\ge 0$) if only

$\begin{array}{ll}\text{Γ}=1-2{(h^2)}^{-1}\varphi (k)+\rho \varphi (k)\ge 0,& (69)\end{array}$

It follows that

$\begin{array}{ll}0\le 1-\text{Γ}=(2{(h^2)}^{-1}-\rho )\varphi (k)\le 1.& (70)\end{array}$

Hence, in Mickens [24], the condition required for positivity is

$\begin{array}{ll}\varphi \left(k\right)\le \frac{h^2}{2-\rho h^2}\left[1-\text{Γ}\right]\text{ and }0\le \text{Γ}<1,\rho h^2\ne 2.& (71)\end{array}$

We investigate next the boundedness by assuming $u_m^n\in \left[0,1\right].$ Equation (68) is rewritten as follows

$\begin{array}{ll}{u}_m^{n+1}=\frac{\text{Γ}{u}_m^n+R\left(u_{m+1}^n+u_{m-1}^n\right)}{1+\left(\frac{\rho \varphi \left(k\right)}{3}\right)\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}.& (72)\end{array}$

where Γ = 1 − 2^(h2)−1ϕ(k) + ρϕ(k), $R=\frac{\varphi (k)}{h^2}.$

Following the idea of Mickens [24], Equation (72) takes the symmetric form if Γ = R. Therefore, it follows that

$\begin{array}{ll}\text{Γ}=\frac{\varphi \left(k\right)}{h^2}=\frac{1}{3}+\frac{\rho \varphi \left(k\right)}{3}.& (73)\end{array}$

We also have from Equation (71)

$\begin{array}{ll}\varphi \left(k\right)\le \frac{h^2}{2-\rho h^2}\left[1-\text{Γ}\right]\Rightarrow \frac{\varphi \left(k\right)}{h^2}\le \frac{1}{2}\left[\frac{1}{1-\frac{\rho h^2}{2}}\right]& (74)\end{array}$

Based on the symmetric condition, we can take

$\begin{array}{ll}\text{Γ}=\frac{\varphi \left(k\right)}{h^2}=\frac{1}{2}\left[\frac{1}{1-\frac{\rho h^2}{2}}\right]& (75)\end{array}$

With regard to the symmetric condition (Equations 73), Equation (72) can be rewritten as

$\begin{array}{ll}{u}_m^{n+1}=\frac{\text{Γ}\left(u_m^n+u_{m+1}^n+u_{m-1}^n\right)}{1+\left(\frac{\rho \varphi \left(k\right)}{3}\right)\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}.& (76)\end{array}$

We know by the assumption that $u_j^n\in \left[0,1\right],$ ∀j. We have

$\begin{array}{ll}0\le \frac{u_m^n+u_{m+1}^n+u_{m-1}^n}{3}\le 1.& (77)\end{array}$

By multiplying Equation (77) by $1-\frac{\rho h^2}{2}$ and dividing by $1-\frac{\rho h^2}{2}$ and expanding, we have

$\begin{array}{ll}\frac{u_m^n+u_{m+1}^n+u_{m-1}^n}{3\left[1-\frac{\rho h^2}{2}\right]}-\left[\frac{\rho h^2}{2}\right]\frac{u_m^n+u_{m+1}^n+u_{m-1}^n}{3\left[1-\frac{\rho h^2}{2}\right]}\le 1.& (78)\end{array}$

From Equation (78), we have

$\begin{array}{l}\frac{u_m^n+u_{m+1}^n+u_{m-1}^n}{3\left[1-\frac{\rho h^2}{2}\right]}\le 1+\rho \frac{h^2}{2}\frac{1}{3}\left[\frac{1}{1-\frac{\rho h^2}{2}}\right](u_m^n+u_{m+1}^n\\ +u_{m-1}^n).& (79)\end{array}$

Let Γ′ = σΓ, σ ≠ 0. Then

$\begin{array}{ll}{\text{Γ}}^{\prime }=\frac{\sigma }{2}\left[\frac{1}{1-\frac{\rho h^2}{2}}\right]=\sigma \frac{\varphi \left(k\right)}{h^2}& (80)\end{array}$

If $\frac{\sigma }{2}=\frac{1}{3}$, then $\sigma =\frac{2}{3}$ and using Equation (79), we have

$\begin{array}{ll}{\text{Γ}}^{\prime }\left(u_m^n+u_{m+1}^n+u_{m-1}^n\right)\le 1+\frac{\rho \varphi \left(k\right)}{3}\left(u_m^n+u_{m+1}^n+u_{m-1}^n\right).& (81)\end{array}$

Hence,

$\begin{array}{ll}0\le u_m^{n+1}=\frac{{\text{Γ}}^{\prime }\left(u_m^n+u_{m+1}^n+u_{m-1}^n\right)}{1+\left(\frac{\rho \varphi \left(k\right)}{3}\right)\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}\le 1.& (82)\end{array}$

Thus, the boundedness of $u_m^{n+1}.$

6.3. Error Estimate

THEOREM 4. Let u ∈ C^{4, 2}(Q) where Q is defined by

$Q=\left\{\left(x,t\right)/a\le x\le b,0<t\le T,T>0,a,b\in ℝ\right\}.$

Assume h and k are such that the Theorem 3 is satisfied and $e_m^n=u_m^n-v_m^n,$ is the defined error of the scheme constructed. NSFD is consistent, and the estimate error $e_m^n$ is defined by

$\begin{array}{l}\left|e_m^n\right|\le E_m^n={\left(3\text{Γ}\right)}^n{E}_m^0+\left(\frac{1-{\left(3\text{Γ}\right)}^n}{1-3\text{Γ}}\right)\\ \left[\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\right]& (83)\end{array}$

where Γ defined in 3 M is defined by $M=\underset{\left\{(x,t)\right\}\in Q}{max}\left\{\left|u_{xxxx}(x,t)\right|,\left|u_{tt}(x,t)\right|\right\},$ and Θ_i, i = 1, 2, 3 such that

$\begin{array}{l}{\text{Θ}}_1\left(\varphi \left(k\right),h\right)=\rho \frac{h^2}{3}u-\rho h^2\frac{\varphi {\left(k\right)}^2}{6}{u}_{tt}\\ {\text{Θ}}_2\left(\varphi \left(k\right),h\right)=-\frac{\varphi {\left(k\right)}^2}{6}+\rho \frac{\varphi \left(k\right)}{3}\\ \left[\frac{{\varphi }^2\left(k\right)}{2}u+h^2\frac{{\varphi }^2}{6}{u}_{xx}+h^4\frac{\varphi {\left(k\right)}^2}{72}{u}_{xxxx}\right]\\ {\text{Θ}}_3\left(\varphi \left(k\right),h\right)=-\varphi \left(k\right)v-\rho h^2\frac{\varphi \left(k\right)}{3}-\rho h^4\frac{\varphi \left(k\right)}{36}{u}_{xxxx}& (84)\end{array}$

and

Θ₁(ϕ(k), h)u_xx + Θ₂(ϕ(k), h)u_ttt + Θ₃(ϕ(k), h)u_t = O(ϕ(k), h) → 0 when ϕ(k) → 0 and h → 0.

Proof.

$\begin{array}{ll}{v}_m^{n+1}=\frac{\text{Γ}{v}_m^n+R\left(v_{m+1}^n+v_{m-1}^n\right)}{1+\left(\frac{\rho \varphi \left(k\right)}{3}\right)\left(v_{m+1}^n+v_m^n+v_{m-1}^n\right)},& (85)\end{array}$

where $\text{Γ}=1-2\frac{\varphi (k)}{h^2}+\rho \varphi (k),R=\frac{\varphi (k)}{h^2}.$ Taylor series expansion of Equation (85) about (t_n, x_m) gives

$\begin{array}{l}\left(v+\varphi \left(k\right)v_t+\frac{{\left(\varphi \left(k\right)\right)}^2}{2}{v}_{tt}+\frac{{\left(\varphi \left(k\right)\right)}^3}{6}{v}_{ttt}+O\left({\left(\varphi \left(k\right)\right)}^4\right)\right)\\ \left(1+\rho \frac{\varphi \left(k\right)}{3}\left\{\begin{array}{l}v+v+h{v}_x+\frac{h^2}{2}{v}_{xx}+\frac{h^3}{6}{v}_{xxx}+\frac{h^4}{24}{v}_{xxxx}\\ +v-h{v}_x+\frac{h^2}{2}{v}_{xx}-\frac{h^3}{6}{v}_{xxx}+\frac{h^4}{24}{v}_{xxxx}\end{array}\right\}\right)\\ =\left(1-\frac{2\varphi \left(k\right)}{h^2}+\rho \varphi \left(k\right)\right)v\\ +\frac{\varphi \left(k\right)}{h^2}\left\{\begin{array}{l}v+h{v}_x+\frac{h^2}{2}{v}_{xx}+\frac{h^3}{6}{v}_{xxx}+\frac{h^4}{24}{v}_{xxxx}\\ +v-h{v}_x+\frac{h^2}{2}{v}_{xx}-\frac{h^3}{6}{v}_{xxx}+\frac{h^4}{24}{v}_{xxxx}\end{array}\right\}\end{array}$

This gives

$\begin{array}{l}v+\varphi \left(k\right)v_t+\frac{{\left(\varphi \left(k\right)\right)}^2}{2}{v}_{tt}+\frac{{\left(\varphi \left(k\right)\right)}^3}{6}{v}_{ttt}\\ +\rho \frac{\varphi \left(k\right)}{3}\left(v+\varphi \left(k\right)v_t+\frac{{\left(\varphi \left(k\right)\right)}^2}{2}{v}_{tt}+\frac{{\left(\varphi \left(k\right)\right)}^3}{6}{v}_{ttt}\right)\\ \left(3v+h^2{v}_{xx}+\frac{h^4}{12}{v}_{xxxx}\right)\\ =v+\left(-\frac{2\varphi \left(k\right)}{h^2}+\rho \varphi \left(k\right)\right)v+\frac{\varphi \left(k\right)}{h^2}\left(2v+h^2{v}_{xx}++\frac{h^4}{12}{v}_{xxxx}\right).& (86)\end{array}$

It follows after division by ϕ(k) and simplification, we have

$\begin{array}{l}{v}_t+\frac{\left(\varphi \left(k\right)\right)}{2}{v}_{tt}+\frac{{\left(\varphi \left(k\right)\right)}^2}{6}{v}_{ttt}\\ +\frac{\rho }{3}\left(\varphi \left(k\right)v_t+\frac{{\left(\varphi \left(k\right)\right)}^2}{2}{v}_{tt}+\frac{{\left(\varphi \left(k\right)\right)}^3}{6}{v}_{ttt}\right)\\ \left(\begin{array}{l}3v+h^2{v}_{xx}+\frac{h^4}{12}{v}_{xxxx}\end{array}\right)\\ +\frac{\rho }{3}v\left(3v\right)+\frac{\rho }{3}v\left(\begin{array}{l}{h}^2{v}_{xx}+\frac{h^4}{12}{v}_{xxxx}\end{array}\right)\\ =v_{xx}+\rho v+\frac{h^2}{12}{v}_{xxxx}.& (87)\end{array}$

It follows that

$\begin{array}{l}{v}_t-v_{xx}+\rho v^2-\rho v\\ =\left\{\begin{array}{l}\left(\frac{h^2}{12}-\rho \frac{h^4}{36}v\right)v_{xxxx}-\left(\frac{\varphi \left(k\right)}{2}+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{v}_{xxxx}+\rho \frac{{\varphi }^2\left(k\right)}{2}v\right)v_{tt}\\ +\left(\rho \frac{h^2}{3}v-\rho h^2\frac{\varphi {\left(k\right)}^2}{6}{v}_{tt}\right)v_{xx}\\ +\left(-\frac{\varphi {\left(k\right)}^2}{6}+\rho \frac{\varphi \left(k\right)}{3}\left[\frac{\varphi {\left(k\right)}^2}{2}v+h^2\frac{{\varphi }^2\left(k\right)}{6}{v}_{xx}+h^4\frac{\varphi {\left(k\right)}^2}{72}{v}_{xxxx}\right]\right)v_{ttt}\\ +\left(-\varphi \left(k\right)v-\rho h^2\frac{\varphi \left(k\right)}{3}-\rho h^4\frac{\varphi \left(k\right)}{36}{v}_{xxxx}\right)v_t\end{array}\right\}& (88)\end{array}$

If ϕ(k) → 0 and h → 0, Equation (88) reduces to $v_t-v_{xx}+\rho v^2-\rho v\to 0.$ Hence, the consistency.

For the simplicity of the proof, we consider the function Θ_i, i = 1, 2, 3 such that

$\begin{array}{l}{\text{Θ}}_1\left(\varphi \left(k\right),h\right)=\rho \frac{h^2}{3}v-\rho h^2\frac{\varphi {\left(k\right)}^2}{6}{v}_{tt}\\ {\text{Θ}}_2\left(\varphi \left(k\right),h\right)=-\frac{\varphi {\left(k\right)}^2}{6}+\rho \frac{\varphi \left(k\right)}{3}\\ \left[\frac{{\varphi }^2\left(k\right)}{2}v+h^2\frac{{\varphi }^2}{6}{v}_{xx}+h^4\frac{\varphi {\left(k\right)}^2}{72}{v}_{xxxx}\right]\\ {\text{Θ}}_3\left(\varphi \left(k\right),h\right)=-\varphi \left(k\right)v-\rho h^2\frac{\varphi \left(k\right)}{3}-\rho h^4\frac{\varphi \left(k\right)}{36}{v}_{xxxx}& (89)\end{array}$

and Θ₁(ϕ(k), h)v_xx + Θ₂(ϕ(k), h)v_ttt + Θ₃(ϕ(k), h)v_t = O(ϕ(k), h) → 0 when ϕ(k) → 0 and h → 0.

The exact discrete equation is

$\begin{array}{l}{u}_m^{n+1}=\frac{\text{Γ}{u}_m^n+R\left(u_{m+1}^n+u_{m-1}^n\right)}{1+\frac{\rho \varphi \left(k\right)}{3}\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}\\ +\left(\frac{\varphi \left(k\right)}{2}+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{u}_{xxxx}\left({\epsilon }_m,t\right)+\rho \frac{{\varphi }^2\left(k\right)}{2}u\right)u_{tt}\left(x,{\tau }_n\right)\\ -\left(\frac{h^2}{12}-\rho \frac{h^4}{36}u\right)u_{xxxx}\left({\epsilon }_m,t\right)& (90)\end{array}$

where $\text{Γ}=1-2\frac{\varphi (k)}{h^2}+\rho \varphi (k),R=\frac{\varphi (k)}{h^2},$ and x_m < ε_m < x_{m + 1} and t_n < τ_n < t_{n + 1}.

We define $e_m^n=u_m^n-v_m^n\equiv e_m^{n+1}=u_m^{n+1}-v_m^{n+1}.$ It follows by considering symmetry condition Γ = R

$\begin{array}{l}{u}_m^{n+1}-v_m^{n+1}\\ =\left\{\begin{array}{l}\frac{\text{Γ}\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}{1+\frac{\rho \varphi \left(k\right)}{3}\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)}-\frac{\text{Γ}\left(v_{m+1}^n+v_m^n+v_{m-1}^n\right)}{1+\frac{\rho \varphi \left(k\right)}{3}\left(v_{m+1}^n+v_m^n+v_{m-1}^n\right)}\\ +\left(\frac{\varphi \left(k\right)}{2}+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{u}_{xxxx}\left({\epsilon }_m,t\right)+\rho \frac{{\varphi }^2\left(k\right)}{2}u\right)u_{tt}\left(x,{\tau }_n\right)\\ -\left(\frac{h^2}{12}-\rho \frac{h^4}{36}u\right)u_{xxxx}\left({\epsilon }_m,t\right)\end{array}\right\}.& (91)\end{array}$

It follows

$\begin{array}{l}{e}_m^{n+1}\\ =\frac{\text{Γ }\left(e_{m+1}^n+e_m^n+e_{m-1}^n\right)}{\left(1+\frac{\rho \varphi \left(k\right)}{3}\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)\right)\left(1+\frac{\rho \varphi \left(k\right)}{3}\left(v_{m+1}^n+v_m^n+v_{m-1}^n\right)\right)}\\ +\left(\frac{\varphi \left(k\right)}{2}+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{u}_{xxxx}\left({\epsilon }_m,t\right)+\rho \frac{\varphi {\left(k\right)}^2}{2}u\right)u_{tt}\left(x,{\tau }_n\right)\\ -\left(\frac{h^2}{12}-\rho \frac{h^4}{36}u\right)u_{xxxx}\left({\epsilon }_m,t\right)& (92)\end{array}$

Let $M=\underset{\left\{(x,t)\right\}\in Q}{max}\left\{\left|u_{xxxx}(x,t)\right|,\left|u_{tt}(x,t)\right|\right\}$ and $E_m^n=\underset{0<m<N}{max}\left\{\left|e_m^n\right|\right\}.$ We have

$\begin{array}{l}\left(1+\frac{\rho \varphi \left(k\right)}{3}\left(u_{m+1}^n+u_m^n+u_{m-1}^n\right)\right)\\ \left(1+\frac{\rho \varphi \left(k\right)}{3}\left(v_{m+1}^n+v_m^n+v_{m-1}^n\right)\right)>1,\end{array}$

$\forall u_i^n,v_i^n\in \left[0,1\right],i=m-1,m,m+1$

and by using Theorem 3, we have

$\begin{array}{l}\left|e_m^{n+1}\right|\le \left|\text{Γ}\right|\left(\left|e_{m+1}^n\right|+\left|e_m^n\right|+\left|e_{m-1}^n\right|\right)\\ +\left|\left(\frac{\varphi \left(k\right)}{2}+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{u}_{xxxx}\left({\epsilon }_m,t\right)+\rho \frac{\varphi {\left(k\right)}^2}{2}u\right)\right|\\ \left|u_{tt}\left(x,{\tau }_n\right)\right|+|-\left(\frac{h^2}{12}-\rho \frac{h^4}{36}u\right)|\left|u_{xxxx}\left({\epsilon }_m,t\right)\right|\\ \le 3\text{ Γ }{E}_m^n+\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M\\ +\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2& (93)\end{array}$

Let

$\begin{array}{l}{E}_m^{n+1}=3\text{ Γ }{E}_m^n+\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M\\ +\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2.\end{array}$

We have

$\begin{array}{l}\left|e_m^{n+1}\right|\le E_m^{n+1}=3\text{ Γ }{E}_m^n+(\frac{\varphi \left(k\right)}{2}\\ +\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36})M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2& (94)\end{array}$

For n = 0, $E_m^1=3\text{Γ}{E}_m^0+\left(\frac{\varphi (k)}{2}+\rho \frac{\varphi {(k)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {(k)}^2}{72}{M}^2$

For n = 1, we have

$\begin{array}{l}{E}_m^2=3\text{ Γ }{E}_m^1+\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M\\ +\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\\ =3\text{Γ}(3\text{ Γ }{E}_m^0+\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M\\ +\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2)\\ +\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\\ ={\left(3\text{ Γ}\right)}^2{E}_m^0+\left(1+3\text{ Γ}\right)\\ \left[\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\right]& (95)\end{array}$

For n = 2, we have

$\begin{array}{l}{E}_m^3=3\text{ Γ }{E}_m^2+\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\\ =3\text{ Γ}(3^2{\text{Γ}}^2+\left(1+3\text{ Γ}\right)[\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M\\ +\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\left]\right)\\ +\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\\ ={\left(3\text{ Γ}\right)}^3{E}_m^0+\left(1+3\text{ Γ}+{\left(3\text{ Γ}\right)}^2\right)[\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M\\ +\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2]& (96)\end{array}$

By recurrence for n, we have

$\begin{array}{l}{E}_m^{n+1}={\left(3\text{ Γ}\right)}^n{E}_m^0+\left({\displaystyle \sum _{i=0}^{i-1}}{\left(3\text{ Γ}\right)}^i\right)\\ \left[\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\right]\\ ={\left(3\text{ Γ}\right)}^n{E}_m^0+\left(\frac{1-{\left(3\text{ Γ}\right)}^n}{1-3\text{ Γ}}\right)\\ \left[\left(\frac{\varphi \left(k\right)}{2}+\rho \frac{\varphi {\left(k\right)}^2}{2}+\frac{h^2}{12}+\rho \frac{h^4}{36}\right)M+\rho h^4\frac{\varphi {\left(k\right)}^2}{72}{M}^2\right]& (97)\end{array}$

7. Explicit Exponential Finite Difference Scheme

The EFFD method was developed by Bhattacharya [41] (primarily called the Bhattacharya method) for the numerical solution of the heat equation. The Exponential Finite Difference method was utilized to solve Burgers' equation and generalized Huxley and Burgers-Huxley equations [42–44]. Later, Macías-Díaz and Ĩnan [45] used modified exponential methods to obtain the solution of the Burgers' equation. Furthermore, Inan et al. [46] utilized the EEFD method for numerical solutions for the Newell-Whitehead-Segel type equations which are very useful in biomathematics. They showed convergence, consistency, and stability of the method.

In this section, we obtain numerical solutions of the equation by EEFD method. The solution domain are discretized into cells as (x_m, t_n) in which x_m = mh, (m = 0, 1, 2, ..., N) and t_n = nk, (n = 0, 1, 2, ...), $h=\text{Δ}x=\frac{b-a}{N}$ is the spatial mesh size and k = Δt is the time step, $u_m^n$ denotes the EEFD approximation and u(x, t) denotes the exact solution.

Dividing Equation (6) by u gives

$\begin{array}{ll}\frac{\partial lnu}{\partial t}=\frac{1}{u}\left(u\left(1-u\right)+\frac{{\partial }^{2}u}{\partial x^2}\right).& (98)\end{array}$

Using the finite difference approximations for derivatives, Equation (98) gives

$\begin{array}{ll}{u}_m^{n+1}=u_m^{n}exp\left\{k\left(1-u_m^n\right)+R\left(\frac{u_{m+1}^n-2u_m^n+u_{m-1}^n}{u_m^n}\right)\right\}& (99)\end{array}$

where $R=\frac{k}{h^2}.$ Equation (99) gives the expression for the EEFD method for Fisher's equation.

7.1. Convergence and Estimate Error

For stability analysis, we require non-iterative exponential time-linearization and iterative exponential quasilinearization techniques for Equation (6) which are found in the discretization of the time derivative, the freezing of the coefficients of the resulting linear ordinary differential equations, and the piecewise analytical solution of these ordinary differential equations. These techniques give three-point finite difference expressions that depend in an exponential manner on either the diffusion, reaction, and transient terms or the diffusion and reaction terms. Following the idea of Ramos [25], we transform (Equations 6, 99) into a linear ordinary differential equation by discretizing the time derivative by means of θ-method [25] and linearizing the nonlinear source term, u(1 − u), with respect to either the previous time level or the previous iteration with Jacobian, $J=\frac{d(u(1-u))}{du}=1-2u$:

a) If the linearization is performed with respect to the previous time level, one obtains

$\begin{array}{l}\frac{u_m^{n+1}-u_m^n}{k}=\theta \frac{u_{m-1}^{n+1}-2u_m^{n+1}+u_{m+1}^{n+1}}{h^2}\\ +\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}\\ +u_m^n\left(1-u_m^n\right)\\ +\theta J_m^n\left(u_m^{n+1}-u_m^n\right)& (100)\end{array}$

which yields a non-iterative time linearization method.

b) If the linearization is performed with respect to the previous iteration, one obtains

$\begin{array}{l}\frac{u_m^{i+1}-u_m^n}{k}=\theta \frac{u_{m-1}^{i+1}-2u_m^{i+1}+u_{m+1}^{i+1}}{h^2}\\ +\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}\\ +\left(1-\theta \right)u_m^n\left(1-u_m^n\right)+\theta u_m^{n+1}\left(1-u_m^{i+1}\right)\\ +\theta J_m^{i+1}\left(u_m^{i+1}-u_m^n\right)& (101)\end{array}$

which corresponds to an iterative quasilinear technique and i = 1, 2, ···, n.

Equations (100) and (101) can be solved in closed form in (x_m−1, x_m+1) subject to the following conditions:

$\begin{array}{ll}u\left(x_{m-1}\right)=u_{m-1},u\left(x_m\right)=u_m,u\left(x_{m+1}\right)=u_{m+1}& (102)\end{array}$

and yield exponential solutions in (x_m−1, x_m+1) which are analytical in that interval and continuous everywhere. Since Equations (100) and (101) are very similar, we will only present in detail exponential methods for Equation (100) in the following subsections.

7.1.1. Time-Linearized Full Exponential Techniques

The piecewise analytical solution of Equation (100) can be rewritten as

$\begin{array}{l}\left(\theta J_m^n-\frac{1}{k}\right)u_m^{n+1}+\theta \frac{u_{m-1}^{n+1}-2u_m^{n+1}+u_{m+1}^{n+1}}{h^2}\\ =\left(\theta J_m^n-\frac{1}{k}\right)u_m^n-\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}\\ -u_m^n\left(1-u_m^n\right).& (103)\end{array}$

Let $\left(\theta J_m^n-\frac{1}{k}\right)u_m^{n+1}+\theta \frac{u_{m-1}^{n+1}-2u_m^{n+1}+u_{m+1}^{n+1}}{h^2}=z_m^{n+1}.$ Equation (103) becomes

$\begin{array}{l}{z}_m^{n+1}=\left(\theta J_m^n-\frac{1}{k}\right)u_m^n-\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}\\ -u_m^n\left(1-u_m^n\right)& (104)\end{array}$

which accounts for diffusion, reaction, and transient effects in the differential operator [25]. The solution of Equation (104) depends on the sign of $D_m^n\equiv \theta J_m^n-\frac{1}{k}.$ We have the following three cases:

a) If $D_m^n=0,$ the solution of Equation (104) subject to the condition (102) gives the finite difference expression

$\begin{array}{ll}{z}_m^{n+1}=\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}.& (105)\end{array}$

b) If $D_m^n=-\theta {({\lambda }_m)}^2<0$, the solution of Equation (104) subject to the condition (102) gives the three-point finite difference expression

$\begin{array}{ll}{z}_m^{n+1}=\frac{D_m^n}{2}\left(\frac{u_{m-1}^n-2cosh\left({\text{λ}}_{m}h\right)u_m^n+u_{m+1}^n}{1-cosh\left({\text{λ}}_{m}h\right)}\right).& (106)\end{array}$

c) If $D_m^n=\theta {({\lambda }_m)}^2>0$, the solution of Equation (104) subject to the condition (102) gives the three-point finite difference expression

$\begin{array}{ll}{z}_m^{n+1}=\frac{D_m^n}{2}\left(\frac{u_{m-1}^n-2cos\left({\text{λ}}_{m}h\right)u_m^n+u_{m+1}^n}{1-cos\left({\text{λ}}_{m}h\right)}\right).& (107)\end{array}$

REMARK 1. We can make the following remarks

a) The values $\theta =\frac{1}{2},1,$ corresponding to the time linearization methods.

b) The quasilinear full exponential corresponds to the iterative solution of

$\begin{array}{l}{z}_m^{i+1}=-\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}\\ -\left(1-\theta \right)u_m^n\left(1-u_m^n\right)-\theta u_m^i\left(1-u_m^i\right)+\theta J_m^i{u}_m^i-\frac{u_m^n}{k}& (108)\end{array}$

for i = 1, 2, ··, n − 1, have analogous solutions to those reported in time-linearized full exponential techniques section and the cases $\theta =\frac{1}{2}$ and 1, corresponding to the quasilinear methods.

7.1.2. Time-Linearized Exponential Techniques

The time-linearized exponential techniques presented in this section consider the following differential operator.

$\begin{array}{l}\theta J_m^n{u}_m^{n+1}+\theta \frac{u_{m-1}^{n+1}-2u_m^{n+1}+u_{m+1}^{n+1}}{h^2}=\theta J_m^n{u}_m^n\\ -\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}+\frac{u_m-u_m^n}{k}\\ -u_m^n\left(1-u_m^n\right)& (109)\end{array}$

Let $\theta J_m^n{u}_m^{n+1}+\theta \frac{u_{m-1}^{n+1}-2u_m^{n+1}+u_{m+1}^{n+1}}{h^2}=z_m^{n+1}.$ The Equation (109) becomes

$\begin{array}{l}{z}_m^{n+1}=\theta J_m^n{u}_m^n-\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}\\ +\frac{u_m-u_m^n}{k}-u_m^n\left(1-u_m^n\right)& (110)\end{array}$

which only accounts for reaction and diffusion processes and whose solutions depend on the sign of $J_m^n$. We have therefore the solution of Equation (110) subject to the condition (Equation 102) gives the solution in form of:

a) If $J_m^n=0$

$\begin{array}{l}{u}_{m-1}-\left(2+\frac{h^2}{k}\right)u_m+u_{m+1}\\ =-\frac{h^2}{k}\left(u_m^n+k{u}_m^n\left(1-u_m^n\right)\right),t^n<t<t^{n+1}.& (111)\end{array}$

b) If $J_m^n=-{({\lambda }_m)}^2<0$

$\begin{array}{l}{u}_{m-1}-\frac{2+\left(k{J}_m^n-1\right)\left(e^{-{\text{λ}}_{m}h}+e^{{\text{λ}}_{m}h}\right)}{k{J}_m^n}{u}_m\\ +u_{m+1}\\ =\frac{2-\left(e^{-{\text{λ}}_{m}h}+e^{{\text{λ}}_{m}h}\right)}{k{J}_m^n}\\ \left[-u_m+k\left(J_m^n{u}_m-u_m^n\left(1-u_m^n\right)\right)\right],\\ t^n<t<t^{n+1}.& (112)\end{array}$

c) If $J_m^n={({\lambda }_m)}^2>0$

$\begin{array}{l}{u}_{m-1}-2\frac{1+\left(k{J}_m^n-1\right)cos\left({\text{λ}}_{m}h\right)}{k{J}_m^n}{u}_m+u_{m+1}\\ =2\frac{1-(cos({\text{λ}}_{m}h)}{k{J}_m^n}\\ \left[-u_m+k\left(J_m^n{u}_m-u_m^n\left(1-u_m^n\right)\right)\right],\\ t^n<t<t^{n+1}.& (113)\end{array}$

REMARK 2. The quasilinear full exponential corresponds to the iterative solution of

$\begin{array}{l}\theta J_m^i{u}_m^{n+1}+\theta \frac{u_{m-1}^{n+1}-2u_m^{n+1}+u_{m+1}^{n+1}}{h^2}=\theta J_m^n{u}_m^n\\ -\left(1-\theta \right)\frac{u_{m-1}^n-2u_m^n+u_{m+1}^n}{h^2}+\frac{u_m-u_m^n}{k}\\ -\left(1-\theta \right)u_m^n\left(1-u_m^n\right)-\theta u_m^i\left(1-u_m^i\right)& (114)\end{array}$

for i = 1, 2···n − 1, have analogous solutions to those reported in time-linearized exponential techniques Section 7.1.2.

REMARK 3. The principal primacy of the full exponential techniques is that they reckon for the reaction, diffusion, and transient terms in finding the homogeneous solution of Equation (103). Notwithstanding, this technique also has the disadvantage that the time step cannot be chosen carelessly because, if $k\ll \frac{1}{\theta J_i^n}$ then the reaction terms do not influence the values of λ_m despite the fact that they do influence the particular solution of Equation (103) through $z_m^{n+1}.$ On top of that, the portion of the transient term influences the ordinary differential operator and, therefore, the homogeneous solution, though the other part influences the particular solution. These obstacles are fully suppressed with the exponential techniques displayed in the time-linearized exponential techniques and the quasilinear full exponential in Section 7.1.2.

THEOREM 5. Ramos [25]. The schemes displayed in Equations (100) and (101) are convergent and convergence is reached when

$\begin{array}{ll}{\left(e_m^n\right)}^2=\frac{1}{N}{\displaystyle \sum _{j=1}^N}{\left(u_j^{i+1}-u_j^i\right)}^2\le 10^{-\left|\alpha \right|}& (115)\end{array}$

where i = 1, 2···n and |α| is an integer obtained from numerical computation, N denotes the number of grid points, and $e_m^n$ is the error which is defined by $e_m^n=u_m^n-v_m^n$.

Proof. The full proof of this theorem is detailed in Ramos [25].

8. Numerical Results

The stability region of FTCS is $k\le \frac{h^2}{2}.$ For h = 0.01, we obtain k ≤ 5 × 10⁻⁵. In the case of NSFD, the condition for positivity gives $\varphi (k)\le \frac{h^2}{2}$ where $\varphi (k)=\frac{1-e^{\lambda k}}{\lambda }.$ We tabulate L₁ and L_∞ errors at certain values of k using ρ = 1, h = 0.01, at time, $T_{max}=2.5\times 10^{-3}$ using FTCS, NSFD, and EEFD schemes at certain various values of time-step size as chosen from T_max/52, T_max/100, T_max/200, ⋯ , T_max/1,800, T_max/1,900 and T_max/2,000. The errors are shown in Table 1.

TABLE 1

Table 1. L₁ and L_∞ errors at some different values of time-step size, k with ρ = 1 at time, $T_{max}=2.5\times 10^{-3}$ with spatial mesh size, h = 0.01 using three methods.

In the case of FTCS, NSFD, minimum L₁, and L_∞ errors occur at k = T_max/2,000 while in the case of EEFD, the errors are least k = T_max/1,600. The least error is of order 10⁻⁹ and 10⁻¹⁰ in the case of NSFD and FTCS respectively while the least error is of order 10⁻¹¹ in the case of EEFD. EEFD is a better scheme than FTCS at all values of k used.

We obtain values for L₁ and L_∞ errors at certain values of k using ρ = 10², h = 0.01, at time, $T_{max}=2.5\times 10^{-3}$ using the three methods in Table 2. The schemes behave differently. In all the three methods FTCS, NSFD, and EEFD, the L₁ and L_∞ errors keep on decreasing as the values of k are decreased gradually from k = T_max/52 to k = T_max/2,000.

TABLE 2

Table 2. L₁ and L_∞ errors at some different values of time-step size, k with ρ = 10² at time, $T_{max}=2.5\times 10^{-3}$ with spatial mesh size, h = 0.01 using three methods.

L₁ and L_∞ errors for the third case using ρ = 10⁴, h = 0.01 are displayed in Table 3. Again, the schemes behave differently from each other. Optimal k using FTCS occurs when k ≊ T_max/1,300. The optimal k using NSFD and EEFD are k ≊ T_max/1,500 and k ≊ T_max/2,000, respectively. Once that optimal is reached the error starts increasing again. The least L₁ and L_∞ errors using NSFD are 4.5235 × 10⁻⁴ and 4.5623 × 10⁻³. The corresponding errors are 1.9347 × 10⁻⁴ and 2.3239 × 10⁻³ when FTCS is used. In the case of EEFD, least L₁ is 6.9852 × 10⁻³ and least L_∞ is 8.3399 × 10⁻².

TABLE 3

Table 3. L₁ and L_∞ errors at some different values of time-step size, k with ρ = 10⁴ at time, $T_{max}=2.5\times 10^{-3}$ with spatial mesh size, h = 0.01 using three methods.

We obtain plots of numerical solution vs. x at time, $T_{max}=2.5\times 10^{-3}$ using three methods FTCS, NSFD, and EEFD in Figure 1.

FIGURE 1

Figure 1. Plot of initial solution, numerical solution of (FTCS), NSFD, and EEFD respectively in (A–C) and exact solution against x at time, $T_{max}=2.5\times 10^{-3}$ using h = 0.01 and optimal value of time step size of k for three different values of ρ namely 1, 10², and 10⁴.

9. Conclusion

We have investigated in this paper the spectral analysis and optimal step sizes for some finite difference methods discretising Fisher's equation. We used three methods namely; FTCS, NSFD, and EEFD in order to solve Fisher's equation with a coefficient of reaction being 10, 10², and 10⁴. We studied the properties of the methods such stability, positivity, and boundedness. This is the one of rare article which includes the estimate errors for the methods studied. Numerical results are displayed at optimal time step size with h = 0.01 for the three cases for the three methods used. We also obtained numerically the rate of convergence as shown in Tables 4–6. We have shown from Tables 1–3 that all the three methods (FTCS, NSFD, and EEFD) perform well for the small coefficient of reaction. It is worthy mentioning that freezing coefficient technique with Von Neumann Stability Analysis only present an approximate stability region for standard methods discretising non linear partial differential equations which might give reason to the standard method in Table 1 to perform better than the NSFD. Furthermore the NSFD in regard to the discrete representation derivative in Mickens [39] rule has nontrivial denominator function and make use of positivity and boundedness conditions. Finally the results are dependent on initial conditions used. For ρ = 1, the difference in L₁ and L_∞ errors from the three methods is very small which lead to conclude that the best methods are FTCS and NSFD. Also for ρ = 10² and ρ = 10⁴, the best method are EEFD and NSFD, respectively. Our results matched with the one found in Lubuma and Roux [47] for numerical experiment for small reaction term. Moreover, Lubuma and Roux [47] proved that NSFD is elementary stable. As NSFD methods, the EEFD displayed in this article do not require any knowledge of the exact solution of the differential equation. Contrast to that, the best finite difference scheme is stable for large grid sizes but costly in inaccuracies at the propagation front.

TABLE 4

Table 4. Rate of convergence in time for FTCS using different values of ρ at time, $T_{max}=2.5\times 10^{-3}$ with h = 0.01.

TABLE 5

Table 5. Rate of convergence in time for NSFD using different values of ρ at time, $T_{max}=2.5\times 10^{-3}$ with h = 0.01.

TABLE 6

Table 6. Rate of convergence in time for EEFD using different values of ρ at time, $T_{max}=2.5\times 10^{-3}$ with h = 0.01.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

The plan of the paper was prepared by AA. The coding was done by KA and BI. All authors helped in writing the manuscript, contributed to the article, and approved the submitted version.

Funding

KA is supported by Postdoctoral Research Fellowship (PDRF) from North-West University (Source of funding: NWU PDRF Fund NW.1G01487). AA is supported by Nelson Mandela University for this research.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

KA is grateful to the North-West University for the Postdoctoral Research Fellowship (PDRF). AA thank Nelson Mandela University for the financial support. Finally, all the authors are grateful to the two reviewers who provided feedback which enabled to improve the paper considerably.

References

1. Gardner CS, Greene JM, Kruskal MD, Miura RM. Method for solving the Korteweg-de Vries equation. Front Neurosci Phys Rev Lett. (1967) 19:1095. doi: 10.1103/PhysRevLett.19.1095

CrossRef Full Text | Google Scholar

2. Weiss J, Tabor M, Carnevale G. The Painlevé property for partial differential equations. J Math Phys. (1983) 155:522–6. doi: 10.1063/1.525721

CrossRef Full Text | Google Scholar

3. Kudryashov NA. On types of nonlinear nonintegrable equations with exact solutions. Phys Lett A. (1991) 24:269–75. doi: 10.1016/0375-9601(91)90481-M

CrossRef Full Text | Google Scholar

4. Lou SY, Huang GX, Ruan HY. Exact solitary waves in a convecting fluid. J Phys A Math Gen. (1991) 24:L587. doi: 10.1088/0305-4470/24/11/003

CrossRef Full Text | Google Scholar

5. Kudryashov NA. Exact solutions of a family of Fisher equations. Theor Math Phys. (1993) 94:211–218. doi: 10.1007/BF01019332

CrossRef Full Text | Google Scholar

6. Kudryashov NA. Numerical convective schemes based on accurate computation of space derivatives. J Comput Phys. (1973) 13:100–13. doi: 10.1016/0021-9991(73)90128-9

CrossRef Full Text | Google Scholar

7. Cooley JW, Lewis P, Welch P. The finite Fourier transform. IEEE Trans Audio Electroacoust. (1969) 17:77–85. doi: 10.1109/TAU.1969.1162036

CrossRef Full Text | Google Scholar

8. Gazdag J, Canosa J. Numerical solution of Fisher's equation. J Appl Probab. (1974) 11:445–57. doi: 10.2307/3212689

CrossRef Full Text | Google Scholar

9. Fisher RA. The wave of advance of advantageous genes. Ann Hum Genet. (1937) 17:355–69. doi: 10.1111/j.1469-1809.1937.tb02153.x

CrossRef Full Text | Google Scholar

10. Moran PAP. The Statistical Processes of Evolutionary Theory. Oxford: Clarendon Press; Oxford University Press (1962).

11. Kendall DG. A Form of Wave Propagation Associated With the Equation of Heat Conduction. Cambridge: Cambridge University Press (1948).

Google Scholar

12. Kolmogorov AN. Étude de l'équation de la diffusion avec croissance de laquantité de matière et son application à un problème biologique. Bull Univ Moskow Ser Internat A. (1937) 1:1–25.

Google Scholar

13. Canosa J. On a nonlinear diffusion equation describing population growth. IBM J Res Dev. (1937) 17:307–13. doi: 10.1147/rd.174.0307

CrossRef Full Text | Google Scholar

14. Jeffrey A, Kakutani T. Weak nonlinear dispersive waves: a discussion centered around the Korteweg-de Vries equation. Siam Rev SIAM. (1972) 14:582–643. doi: 10.1137/1014101

CrossRef Full Text | Google Scholar

15. Kudryashov NA, Zakharchenko AS. A note on solutions of the generalized Fisher equation. Appl Math Lett. (2014) 32:53–6. doi: 10.1016/j.aml.2014.02.009

CrossRef Full Text | Google Scholar

16. Newman TJ, Kolomeisky EB, Antonovics J. Population dynamics with global regulation: the conserved Fisher equation. Phys Rev Lett. (2004) 92:228103. doi: 10.1103/PhysRevLett.92.228103

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Chen BM, Kojouharov HV. Non-standard numerical methods applied to subsurface biobarrier formation models in porous media. Bull Math Biol. (1999) 61:779–98. doi: 10.1006/bulm.1999.0113

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yatat V, Couteron P, Dumont Y. Spatially explicit modelling of tree-grass interactions in fire-prone savannas: a partial differential equations framework. Ecol Complex. (2018) 36:290–313. doi: 10.1016/j.ecocom.2017.06.004

CrossRef Full Text | Google Scholar

19. Wang XY. Exact and explicit solitary wave solutions for the generalised Fisher equation. Phys Lett A. (1988) 131:277–9. doi: 10.1016/0375-9601(88)90027-8

CrossRef Full Text | Google Scholar

20. Chandraker V, Awasthi A, Jayaraj S. A numerical treatment of Fisher equation. Procedia Eng. (2015) 131:1256–62. doi: 10.1016/j.proeng.2015.11.481

CrossRef Full Text | Google Scholar

21. Ismail HN, Raslan K, Rabboh AAA. Adomian decomposition method for Burger's-Huxley and Burger's-Fisher equations. Appl Math Comput. (2004) 159:291–301. doi: 10.1016/j.amc.2003.10.050

CrossRef Full Text | Google Scholar

22. Moghimi M, Hejazi FSA. Variational iteration method for solving generalized Burger-Fisher and Burger equations. Chaos Solitons Fractals. (2007) 33:1756–61. doi: 10.1016/j.chaos.2006.03.031

CrossRef Full Text | Google Scholar

23. Fahmy ES. Travelling wave solutions for some time-delayed equations through factorizations. Chaos Solitons Fractals. (2008) 38:1209–16. doi: 10.1016/j.chaos.2007.02.007

CrossRef Full Text | Google Scholar

24. Mickens RE. Relation between the time and space step-sizes in nonstandard finite-difference schemes for the Fisher equation. Numer Methods Partial Differ Equ. (1997) 13:51–5. doi: 10.1002/(SICI)1098-2426(199701)13:1<51::AID-NUM4>3.0.CO;2-L

CrossRef Full Text | Google Scholar

25. Ramos JI. Exponential methods for one-dimensional reaction-diffusion equations. Math Methods Appl Sci. (2005) 170:380–98. doi: 10.1016/j.amc.2004.12.003

CrossRef Full Text | Google Scholar

26. Anguelov R, Kama P, Lubuma JMS. On non-standard finite difference models of reaction-diffusion equations. J Comput Appl Math. (2005) 175:11–29. doi: 10.1016/j.cam.2004.06.002

CrossRef Full Text | Google Scholar

27. Hagstrom T, Keller HB. The numerical calculation of traveling wave solutions of nonlinear parabolic equations. SIAM J Scientific Stat Comput. (1986) 7:978–88. doi: 10.1137/0907065

CrossRef Full Text | Google Scholar

28. Li S, Petzold L, Ren Y. Stability of moving mesh systems of partial differential equations. SIAM J Scientific Comput. (1998) 20:719–38. doi: 10.1137/S1064827596302011

CrossRef Full Text | Google Scholar

29. Kudryashov NA. One method for finding exact solutions of nonlinear differential equations. Commun Nonlinear Sci Numer Simulat. (2012) 17:2248–53. doi: 10.1016/j.cnsns.2011.10.016

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Polyanim AD, Zaitsev ZF. Handbook of Nonlinear Partial Differential Equations. 2nd ed. Chapman and Hall/CRC Press. (2011). p. 1912.

Google Scholar

31. Qiu Y, Sloan DM. Numerical solution of Fisher's equation using a moving mesh method. J Comput Phys. (1998) 146:726–746. doi: 10.1006/jcph.1998.6081

CrossRef Full Text | Google Scholar

32. Miyata T, Sakai Y. Vectorized total variation defined by weighted L infinity norm for utilizing inter channel dependency. In: 2012 19th IEEE International Conference on Image Processing. Orlando, FL: IEEE (2012).

Google Scholar

33. Sutton OJ. Long time L^∞(L²) a posteriori error estimates for fully discrete parabolic problems. arXiv[Preprint].arXiv:180303207. (2018). doi: 10.1093/imanum/dry078

CrossRef Full Text | Google Scholar

34. Chen-Charpentier BM, Kojouharov HV. An unconditionally positivity preserving scheme for advection-diffusion reaction equations. Math Comput Model. (2013) 57:2177–85. doi: 10.1016/j.mcm.2011.05.005

CrossRef Full Text | Google Scholar

35. Agbavon K, Appadu A, Khumalo M. On the numerical solution of Fisher's equation with coefficient of diffusion term much smaller than coefficient of reaction term. Adv Diff Equat. (2019) 2019:1–33. doi: 10.1186/s13662-019-2080-x

CrossRef Full Text | Google Scholar

36. Durran DR. Numerical Methods for Fluid Dynamics (With Applications to Geophysics). Springer Science & Business Media (2010).

Google Scholar

37. Taha TR, Ablowitz MI. Analytical and numerical aspects of certain nonlinear evolution equations, III. Numerical, Korteweg-de Vries equation. J Comput Phys. (1984) 55:231–53. doi: 10.1016/0021-9991(84)90004-4

CrossRef Full Text | Google Scholar

38. Zabusky NJ, Kruskal MD. Interaction of solitons in a collisionless plasma and the recurrence of initial states. Phys Rev Lett. (1965) 15:240. doi: 10.1103/PhysRevLett.15.240

CrossRef Full Text | Google Scholar

39. Mickens RE. Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations. J Diff Equat Appl. (2005) 11:645–53. doi: 10.1080/10236190412331334527

CrossRef Full Text | Google Scholar

40. Kumar Verma A, Kayenat S. On the stability of Micken's type NSFD schemes for generalized Burgers Fisher equation. Chaos Solitons Fractals. (2019) 25:1706–37. doi: 10.1080/10236198.2019.1689236

CrossRef Full Text | Google Scholar

41. Bhattacharya MC. An explicit conditionally stable finite difference equation for heat conduction problems. Int J Numer Methods Eng. (1985) 21:239–65. doi: 10.1002/nme.1620210205

CrossRef Full Text | Google Scholar

42. Handschuh RF, Keith TG. Applications of an exponential finite-difference technique. Numer Heat Transfer A Appl. (1992) 22:363–78. doi: 10.1080/10407789208944773

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Ĩnan B, Bahadır AR. An explicit exponential finite difference method for the Burger's equation. Eur Int J Sci Technol. (2013) 2:61–72.

Google Scholar

44. Ĩnan B, Bahadır AR. Finite difference methods for the generalized Huxley and Burgers-Huxley equations. Kuwait J Sci. (2017) 44:20–7.

Google Scholar

45. Macías-Díaz JE, Ĩnan B. Numerical efficiency of some exponential methods for an advection-diffusion equation. Int J Comput Math. (2019) 96:1005–29. doi: 10.1080/00207160.2018.1478416

CrossRef Full Text | Google Scholar

46. Inan B, Osman M, Ak T, Baleanu D. Analytical and numerical solutions of mathematical biology models: the newell-whitehead-segel and allen-cahn equations. Math Methods Appl Sci. (2020) 43:2588–600. doi: 10.1002/mma.6067

CrossRef Full Text | Google Scholar

47. Lubuma J, Roux A. An improved theta-method for systems of ordinary differential equations. J Diff Equat Appl. (2003) 9:1023–35. doi: 10.1080/1023619031000146904

CrossRef Full Text | Google Scholar