Front. Appl. Math. Stat. Frontiers in Applied Mathematics and Statistics Front. Appl. Math. Stat. 2297-4687 Frontiers Media S.A. 10.3389/fams.2022.921170 Applied Mathematics and Statistics Original Research Convergence Analysis and Approximate Optimal Temporal Step Sizes for Some Finite Difference Methods Discretising Fisher's Equation Agbavon Koffi Messan 1 Appadu Appanah Rao 2 * Inan Bilge 3 Tenkam Herve Michel 1 1Department of Mathematics and Applied Mathematics, North-West University Potchefstroom Campus, Potchefstroom, South Africa 2Department of Mathematics and Applied Mathematics, Nelson Mandela University, Gqeberha, South Africa 3Department of Basic Sciences of Engineering, Faculty of Engineering and Natural Sciences, İskenderun Technical University, Hatay, Turkey

Edited by: Raluca Eftimie, University of Franche-Comté, France

Reviewed by: Soheil Salahshour, Bahcesehir University, Turkey; Oluwaseun F. Egbelowo, University of Texas at Austin, United States

 *Correspondence: Appanah Rao Appadu rao.appadu@mandela.ac.za

This article was submitted to Mathematical Biology, a section of the journal Frontiers in Applied Mathematics and Statistics

 13 07 2022 2022 8 921170 15 04 2022 09 06 2022 Copyright © 2022 Agbavon, Appadu, Inan and Tenkam. 2022 Agbavon, Appadu, Inan and Tenkam

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

In this study, we obtain a numerical solution for Fisher's equation using a numerical experiment with three different cases. The three cases correspond to different coefficients for the reaction term. We use three numerical methods namely; Forward-Time Central Space (FTCS) scheme, a Nonstandard Finite Difference (NSFD) scheme, and the Explicit Exponential Finite Difference (EEFD) scheme. We first study the properties of the schemes such as positivity, boundedness, and stability and obtain convergence estimates. We then obtain values of L1 and L∞ errors in order to obtain an estimate of the optimal time step size at a given value of spatial step size. We determine if the optimal time step size is influenced by the choice of the numerical methods or the coefficient of reaction term used. Finally, we compute the rate of convergence in time using L1 and L∞ errors for all three methods for the three cases.

 Fisher's equation FTCS NSFD EEFD optimal convergence estimate rate of convergence coefficient of reaction 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

ut=ux(ulux)+ua(1-ub)

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

ut=D∇2u+r(t)u(t)-K(t)u(t)2

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

r(t)=K(t)∫ddx u(x,t)

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

ut=D∇2u+R(u),

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

ut=Duxx+R(u).

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

ut=uxx+R(u),

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

limx→+|ϵ|∞u(x,t)={1    if    ϵ=10    if    ϵ=-1,     u(x,0)=u0(x).

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

u0(x)~exp(-α)    when    x→∞,

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

c(α)={α+1α,    α≤1,2,            α≥1. 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

ut=ux(ul ux)+ua(1-ub),

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, ux(ul ux) on the right-hand side of Equation (11) stands for the growth of population. The term ul 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

s2ddξ(vldvdξ)+va-va+b+cdvdξ=0.

For l ≠ 0, vl = w. Replacing v by w1l in Equation (12) gives

s2lwξ2+s2w wξξ-lwa+b+l-1l+wa+l-1l+ξ wξ=0.

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

w(ξ)=∑j=0PPjQj(ξ),      Q(ξ)=11+eξ-ξ0

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

Qξ=Q-Q2.

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

a+b+l-1l=2+2P.

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

w(ξ)={P0+P1Q(ξ)P0+P1Q(ξ)+P2Q2(ξ)

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

s2lwξ2+s2w wξξ-lw4+w4-bl+ξ wξ=0.

We have the following solutions

w(x,t)={No exact solutions,     b=l and b=3l,±1±2Q(±2l2l+1x±2l2l+1t),  b=2l,±1±2Q(±2l2l+1x±4l(l+1)2l+1t),  b=4l,  or±i±2iQ(±2l2l+1x±4l(l+1)2l+1i t),  b=4l.

We finally get as exact solution

u(x,t)={±1±2Q(±2l2l+1x±2l2l+1t)l,  V1=±12l+1,±1±2Q(±2l2l+1x±4l(l+1)2l+1t)l,  V2=±2(l+1)2l+1,or±i±2iQ(±2l2l+1x±4l(l+1)2l+1i t)l,  V2=±2(l+1)2l+1,

where V1 and V2 stand for velocity.

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

u(x,t)={No exact solutions,     b=2l and b=3l,2(3l+2)l+1(Q(±ll+1i x)-Q2(±ll+1i x))l,b=l.

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

ut=uxx+ρu(1-u),

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

u(x,t)=[1+cexp(ρ6x-5ρ6t)]-2,

where c=5ρ/6 stands for the wave speed with the minimum value, 2ρ.

 4. Numerical Experiments

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

Solve ut=uxx+ρ u(1-u),

for x ∈ [−0.2, 0.8] and t ∈ [0, Tmax] where Tmax=2.5×10-3. The initial data is given by u(x,0)=[1+exp(ρ6x)]-2.

The exact solution is given by u(x,t)=[1+exp(ρ6x-5ρ6t)]-2.

The boundary conditions are as follows u(-0.2,t)=[1+exp(-0.2ρ6-5ρ6t)]-2andu(0.8,t)=[1+exp(0.8ρ6-5ρ6t)]-2.

We consider three cases ρ namely; 1, 102, 104, and obtain a solution at time, t = Tmax.

Definition 1. Miyata and Sakai [32]. For a vector x¯∈ℝN,∥x¯∥1=∑i=1N|x¯i| and ∥x¯∥∞=max{|x¯i|,i=1,···,N}.

Definition 2. Sutton [33]. Suppose {tn}0N forms a partition of [0, T], with tn = nΔt for n = 0, ···, N, where Δt = T/N. Suppose a vector x¯∈ℝN, defined by

∥x¯∥LP(0,tn)={(∥x¯∥LP(0,tn-1)+τ(x¯n)p)1p  for p∈[0,∞),max{∥x¯∥LP(0,tn-1),x¯n}  for p=∞.

The rate of convergence with respect to time is defined by

ratei(t)=log(x¯i(t))-log(x¯i-1(t))log(Δti)-log(Δti-1). 5. Forward in Time Central Space

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

umn+1-umnk=um+1n-2umn+um-1nh2+ρ umn(1-umn),

which leads to

umn+1=(1-2R)umn+k ρ umn(1-umn)+R(um+1n+um-1n),

or

umn+1=(1-2R+k ρ)umn-k ρ(umn)2+Rum+1n+Rum-1n,

where R=kh2.

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

k≤h22

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, ut + 6uux + uxxx = 0 is

umn+1-umn-12k+6(um+1n+umn+um-1n3)(um+1n-um-1n2h)                      +12h3(um+2n-2um+1n+2um-1n-um-2n)=0.

Taha and Ablowitz [37] express u ux as umax ux and utilize the ansatz

umn=ξneImw where w stands for the phase angle. They obtain

(ξ-ξ-1)(2k)-1+(h)-1(6umax)Isin(w)+(2h3)-1(e2Iw-2eIw+2e-Iw-e-2Iw)=0,

which can be rewritten as

ξ = ξ-1-(h)-1(12kumax)Isin(w)-(h3)-1k(e2Iw-2eIw+2e-Iw-e-2Iw)

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

(h)-1k|2umax-(h2)-1|≤2(33)-1.

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

umn+1=(1-(h2)-12k)umn+(h2)-1(k)(um+1n+um-1n)+k ρ umn-k ρ(umn)2.

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

ξ=1-(h2)-1(2k)(1-cos(w))+k ρ(1-umax),

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

ξ=1-(h2)-1(4k)sin2(w2).

The stability is obtained by solving |ξ| ≤ 1 for w ∈ [−π, π], and we obtain k≤h22.

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

u+kut+k22utt+k36uttt+O(k4)=(1-(h2)-1(2k)+k ρ)u-k ρu2+(h2)-1(k)(u+hux+h22uxx+h36uxxx+h424uxxxx+O(h5))+(h2)-1(k)(u-hux+h22uxx-h36uxxx+h424uxxxx+O(h5)),

which can be written as

ut-uxx-ρu(1-u)=-k2utt-k26uttt+h212uxxxx+O(k4)+O(h5).

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

 5.2. Error Estimates

Theorem 2. Let u ∈ C4, 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, Emn for Equation (28) is given by

Emn≤(1+3k ρ)nEm0+19Mh2ρ k[(1+3k ρ)n-1]

where M=max{(x,t)}∈Q{|uxxxx(x,t)|,|utt(x,t)|} and Θ such that Θ(k, h)uttt = O(k, h) → 0, for k, h → 0.

Proof.

Forward in Time Central Space scheme is given by

umn+1=(1-2R+k ρ)umn-k ρ(umn)2+R um+1n+R um-1n.

Taylor series expansion about (n, m) gives

7v+kvt+k22vtt+k36vttt+O(k4)=(1-(h2)-1(2k)+kρ)v-kρv2+(h2)-1(k)(v+hvx+h22vxx+h36vxxx+h424vxxxx+O(h5))+(h2)-1(k)(v-hvx+h22vxx-h36vxxx+h424vxxxx+O(h5)),

which can be rewritten as

vt-vxx-ρv(1-v)=-k2vtt-k26vttt+h212vxxxx+...,

and let Θ(k,h)vttt=O(k,h)=-k26vttt→0 for k, h → 0. The exact discrete equation is

umn+1=(1-2R)umn+kρ umn(1-umn)+R(um+1n+um-1n)+k2utt(x,τn)-h212uxxxx(Xm,t)

where xm < Xm < xm + 1 and tn < τn < tn + 1. We define

emn=umn-vmn⇒emn+1=umn+1-vmn+1.

It follows that

emn+1=(1-2R)(umn-vmn)+k ρ umn(1-umn)-k ρ vmn(1-vmn)+R(um+1n+um-1n)-R(vm+1n+vm-1n)+k2utt(x,τn)-h212uxxxx(Xm,t).

We have

emn+1=(1-2R)(umn-vmn)+k ρ umn-k ρ(umn)2-k ρ vmn+k ρ(vmn)2+R(um+1n-vm+1n)-R(um-1n-vm-1n)+k2utt(x,τn)-h212uxxxx(Xm,t),

which can be rewritten as

emn+1=(1-2R)(umn-vmn)+k ρ(umn-vmn)-k ρ(umn-vmn)(umn+vmn)+R(um+1n-vm+1n)-R(um-1n-vm-1n)+k2utt(x,τn)-h212uxxxx(Xm,t).

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

|emn+1|≤|1-2R||emn|+|k ρ||emn|+|k ρ||emn||umn+vmn|+R|em+1n|+R|em-1n|+M(k2+h212),

where M=max{(x,t)∈Q}{|uxxxx(x,t)|,|utt(x,t)|}. Since 0≤umn≤1 and 0≤vmn≤1, based on numerical experiment chosen, we have

|emn+1|≤|1-2R||emn|+k ρ|emn|+2k ρ|emn|+R|em+1n|+R|em-1n|+M(k2+h212).

Let Emn=max0<m<N{|emn|)|}. We have

|emn+1|≤(|1-2R|+k ρ+2k ρ+2R)|Emn|+M(k2+h212),

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

|emn+1|≤(1+3k ρ)Emn              +M(k2+h212).

Let Emn+1=(1+3kρ)Emn+(k2+h212) M. We have

For n = 0, Em1=(1+3kρ)Em0+(k2+h212) M.

For n = 1, we have Em2=(1+3k ρ)Em1+(k2+h212) M =(1+3k ρ)2Em0+(1+3k ρ)1(k2+h212)M+(1+3k ρ)0(k2+h212) M =(1+3k ρ)2Em0+[(1+3k ρ)1    +(1+3k ρ)0](k2+h212) M

For n = 2, we have Em3=(1+3k ρ)3Em2+(k2+h212) M =(1+3k ρ)3Em0+(1+3k ρ)2(k2+h212)M+(1+3k ρ)1(k2+h212) M +(1+3k ρ)0(k2+h212) M=(1+3k ρ)3Em0+[(1+3k ρ)2+(1+3k ρ)1+(1+3k ρ)0](k2+h212) M

For n, we have       Emn=(1+3k ρ)Emn-1+(k2+h212) M=(1+3k ρ)nEm0+[(1+3k ρ)n-1+(1+3k ρ)n-2+··+(1+3k ρ)1+(1+3k ρ)0](k2+h212) M=(1+3k ρ)nEm0+[∑i=0n-1(1+3k ρ)i](k2+h212) M=(1+3k ρ)nEm0+[1-(1+3k ρ)n1-(1+3k ρ)](k2+h212) M=(1+3k ρ)nEm0-13k ρ[1-(1+3k ρ)n](k2+h212) M

Hence, for k≤h22, we can also write Emn≤(1+3k ρ)nEm0+19M h2ρ k[(1+3k ρ)n-1]

 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.,

 ut≈umn+1-umnϕ(Δt,λ)

where

ϕ(Δt,λ)=Δt+O(Δt2). 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]

{u′=λu,u(0)=u0,λ≠0,u′=λu(1-u),u(0)=u0,λ>0,

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

{u(tn+1)=u0eλtn+1,u(tn+1)=u0e-λtn+1+(1-e-λtn+1)u0.

Let u(tn)=un. We have

{un+1-un(eλ Δt-1)λ-1=λ un,un+1-un(eλ Δt-1)λ-1=λ un(1-un).

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

ϕ(Δt,λ)=eλ Δt-1λor ϕ(Δt,λ)=1-e-λ Δtλ

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

 um2≈umum+1, ,um2≈(um-1+um+um-13)um,

and

u3≈2um3-um2um+1,um3≈um-1umum+1.

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

for Equation (22) is

   umn+1-umnϕ(Δt)=um+1n-2umn+um-1n(Δx)2+ρumn-ρ(um+1n+umn+um-1n3)umn+1,

where

ϕ(Δt)=ϕ(k)=1-e-λkλ;(Δx)2=h2.

The Equation (66) gives the following single expression

umn+1=(1-2ϕ(k)h2+ρϕ(k))umn+ϕ(k)h2(um+1n+um-1n)1+ρϕ(k)3(um+1n+umn+um-1n). 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) ϕ(k)≤h22-ρh2[1-Γ] with Γ = 1 − 2(h2)−1ϕ(k) + ρϕ(k),

(b) For umi∈[0,1], ∀ i. Γ=(h2)-1(ϕ(k))=12[11-ρh22] and Γ′ = σ Γ.

Proof.

We assume u(x, 0) = h(x) ∈ [0, 1]. We have, therefore, u(x, t) ∈ [0, 1] [24]. We assume also umn≥0. The quantity umn+1 from Equation (68) is positive (umn+1≥0) if only

Γ=1-2(h2)-1ϕ(k)+ρϕ(k)≥0,

It follows that

0≤1-Γ=(2(h2)-1-ρ)ϕ(k)≤1.

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

ϕ(k)≤h22-ρh2[1-Γ] and 0≤Γ<1,ρ h2≠2.

We investigate next the boundedness by assuming umn∈[0,1]. Equation (68) is rewritten as follows

umn+1=Γumn+R(um+1n+um-1n)1+(ρϕ(k)3)(um+1n+umn+um-1n).

where Γ = 1 − 2(h2)−1ϕ(k) + ρϕ(k), R=ϕ(k)h2.

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

Γ=ϕ(k)h2=13+ρϕ(k)3.

We also have from Equation (71)

ϕ(k)≤h22-ρh2[1-Γ]⇒ϕ(k)h2≤12[11-ρh22]

Based on the symmetric condition, we can take

Γ=ϕ(k)h2=12[11-ρh22]

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

umn+1=Γ(umn+um+1n+um-1n)1+(ρϕ(k)3)(um+1n+umn+um-1n).

We know by the assumption that ujn∈[0,1], ∀j. We have

0≤umn+um+1n+um-1n3≤1.

By multiplying Equation (77) by 1-ρh22 and dividing by 1-ρh22 and expanding, we have

umn+um+1n+um-1n3[1-ρh22]-[ρh22]umn+um+1n+um-1n3[1-ρh22]≤1.

From Equation (78), we have

umn+um+1n+um-1n3[1-ρh22]≤1+ρh2213[11-ρh22](umn+um+1n+um-1n).

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

Γ′=σ2[11-ρh22]=σϕ(k)h2

If σ2=13, then σ=23 and using Equation (79), we have

Γ′(umn+um+1n+um-1n)≤1+ρϕ(k)3(umn+um+1n+um-1n).

Hence,

0≤umn+1=Γ′(umn+um+1n+um-1n)1+(ρϕ(k)3)(um+1n+umn+um-1n)≤1.

Thus, the boundedness of umn+1.

 6.3. Error Estimate

Theorem 4. Let u ∈ C4, 2(Q) where Q is defined by

Q={(x,t)/a≤x≤b,0<t≤T,T>0,a,b∈ℝ}.

Assume h and k are such that the Theorem 3 is satisfied and emn=umn-vmn, is the defined error of the scheme constructed. NSFD is consistent, and the estimate error emn is defined by

   |emn|≤Emn=(3Γ)nEm0+(1-(3Γ)n1-3Γ)[(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2]

where Γ defined in 3 M is defined by M=max{(x,t)}∈Q{|uxxxx(x,t)|,|utt(x,t)|}, and Θi, i = 1, 2, 3 such that

Θ1(ϕ(k),h)=ρh23u-ρ h2ϕ(k)26uttΘ2(ϕ(k),h)=-ϕ(k)26+ρϕ(k)3[ϕ2(k)2u+h2ϕ26uxx+h4ϕ(k)272uxxxx]Θ3(ϕ(k),h)=-ϕ(k)v-ρ h2ϕ(k)3-ρ h4ϕ(k)36uxxxx

and

Θ1(ϕ(k), h)uxx + Θ2(ϕ(k), h)uttt + Θ3(ϕ(k), h)ut = O(ϕ(k), h) → 0 when ϕ(k) → 0 and h → 0.

Proof.

vmn+1=Γvmn+R(vm+1n+vm-1n)1+(ρϕ(k)3)(vm+1n+vmn+vm-1n),

where Γ=1-2ϕ(k)h2+ρϕ(k),R=ϕ(k)h2. Taylor series expansion of Equation (85) about (tn, xm) gives

(v+ϕ(k)vt+(ϕ(k))22vtt+(ϕ(k))36vttt+O((ϕ(k))4))(1+ρϕ(k)3{v+v+h vx+h22vxx+h36vxxx+h424vxxxx+v-h vx+h22vxx-h36vxxx+h424vxxxx})=(1-2ϕ(k)h2+ρϕ(k))v+ϕ(k)h2{v+h vx+h22vxx+h36vxxx+h424vxxxx+v-h vx+h22vxx-h36vxxx+h424vxxxx}

This gives

v+ϕ(k)vt+(ϕ(k))22vtt+(ϕ(k))36vttt+ρϕ(k)3(v+ϕ(k)vt+(ϕ(k))22vtt+(ϕ(k))36vttt)(3v+h2vxx+h412vxxxx)=v+(-2ϕ(k)h2+ρϕ(k))v+ϕ(k)h2(2v+h2vxx++h412vxxxx).

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

vt+(ϕ(k))2vtt+(ϕ(k))26vttt+ρ3(ϕ(k)vt+(ϕ(k))22vtt+(ϕ(k))36vttt)(3v+h2vxx+h412vxxxx)+ρ3v(3v)+ρ3v(h2vxx+h412vxxxx)=vxx+ρv+h212vxxxx.

It follows that

vt-vxx+ρ v2-ρ v={(h212-ρh436v)vxxxx-(ϕ(k)2+ρ h4ϕ(k)272vxxxx+ρϕ2(k)2v)vtt+(ρh23v-ρ h2ϕ(k)26vtt)vxx+(-ϕ(k)26+ρϕ(k)3[ϕ(k)22v+h2ϕ2(k)6vxx+h4ϕ(k)272vxxxx])vttt+(-ϕ(k)v-ρ h2ϕ(k)3-ρ h4ϕ(k)36vxxxx)vt}

If ϕ(k) → 0 and h → 0, Equation (88) reduces to vt-vxx+ρv2-ρv→0. Hence, the consistency.

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

Θ1(ϕ(k),h)=ρh23v-ρ h2ϕ(k)26vttΘ2(ϕ(k),h)=-ϕ(k)26+ρϕ(k)3[ϕ2(k)2v+h2ϕ26vxx+h4ϕ(k)272vxxxx]Θ3(ϕ(k),h)=-ϕ(k)v-ρ h2ϕ(k)3-ρ h4ϕ(k)36vxxxx

and Θ1(ϕ(k), h)vxx + Θ2(ϕ(k), h)vttt + Θ3(ϕ(k), h)vt = O(ϕ(k), h) → 0 when ϕ(k) → 0 and h → 0.

The exact discrete equation is

umn+1=Γumn+R(um+1n+um-1n)1+ρϕ(k)3(um+1n+umn+um-1n)+(ϕ(k)2+ρ h4ϕ(k)272uxxxx(εm,t)+ρϕ2(k)2u)utt(x,τn)-(h212-ρh436u)uxxxx(εm,t)

where Γ=1-2ϕ(k)h2+ρϕ(k),R=ϕ(k)h2, and xm < εm < xm + 1 and tn < τn < tn + 1.

We define emn=umn-vmn≡emn+1=umn+1-vmn+1. It follows by considering symmetry condition Γ = R

umn+1-vmn+1={Γ(um+1n+umn+um-1n)1+ρϕ(k)3(um+1n+umn+um-1n)-Γ(vm+1n+vmn+vm-1n)1+ρϕ(k)3(vm+1n+vmn+vm-1n)+(ϕ(k)2+ρ h4ϕ(k)272uxxxx(εm,t)+ρϕ2(k)2u)utt(x,τn)-(h212-ρh436u)uxxxx(εm,t)}.

It follows

emn+1=Γ (em+1n+emn+em-1n)(1+ρϕ(k)3(um+1n+umn+um-1n))(1+ρϕ(k)3(vm+1n+vmn+vm-1n))+(ϕ(k)2+ρ h4ϕ(k)272uxxxx(εm,t)+ρϕ(k)22u)utt(x,τn)-(h212-ρh436u)uxxxx(εm,t)

Let M=max{(x,t)}∈Q{|uxxxx(x,t)|,|utt(x,t)|} and Emn=max0<m<N{|emn|}. We have

(1+ρϕ(k)3(um+1n+umn+um-1n))(1+ρϕ(k)3(vm+1n+vmn+vm-1n))>1,

∀uin,vin∈[0,1],i=m-1,m,m+1

and by using Theorem 3, we have

|emn+1|≤|Γ|(|em+1n|+|emn|+|em-1n|)+|(ϕ(k)2+ρ h4ϕ(k)272uxxxx(εm,t)+ρϕ(k)22u)||utt(x,τn)|+|-(h212-ρh436u)||uxxxx(εm,t)|≤3 Γ Emn+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2

Let

Emn+1=3 Γ Emn+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M               +ρ h4ϕ(k)272M2.

We have

|emn+1|≤Emn+1=3 Γ Emn+(ϕ(k)2+ρϕ(k)22+h212+ρh436)M+ρ h4ϕ(k)272M2

For n = 0, Em1=3ΓEm0+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρh4ϕ(k)272M2

For n = 1, we have

Em2=3 Γ Em1+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2=3Γ(3 Γ Em0+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2)+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2=(3 Γ)2Em0+(1+3 Γ)[(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2]

For n = 2, we have

Em3=3 Γ Em2+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2=3 Γ(32Γ2+(1+3 Γ)[(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2])+(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2=(3 Γ)3Em0+(1+3 Γ+(3 Γ)2)[(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2]

By recurrence for n, we have

Emn+1=(3 Γ)nEm0+(∑i=0i-1(3 Γ)i)    [(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2]    =(3 Γ)nEm0+(1-(3 Γ)n1-3 Γ)    [(ϕ(k)2+ρϕ(k)22+h212+ρh436) M+ρ h4ϕ(k)272M2] 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 (xm, tn) in which xm = mh, (m = 0, 1, 2, ..., N) and tn = nk, (n = 0, 1, 2, ...), h=Δx=b-aN is the spatial mesh size and k = Δt is the time step, umn denotes the EEFD approximation and u(x, t) denotes the exact solution.

Dividing Equation (6) by u gives

∂ln u∂t=1u(u(1-u)+∂2u∂x2).

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

umn+1=umnexp{k(1-umn)+R(um+1n-2 umn+um-1numn)}

where R=kh2. 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=d(u(1-u))du=1-2u:

a) If the linearization is performed with respect to the previous time level, one obtains umn+1-umnk=θum-1n+1-2 umn+1+um+1n+1h2                       +(1-θ)um-1n-2 umn+um+1nh2                       +umn(1-umn)                       +θJmn(umn+1-umn)

which yields a non-iterative time linearization method.

b) If the linearization is performed with respect to the previous iteration, one obtains umi+1-umnk=θum-1i+1-2 umi+1+um+1i+1h2+(1-θ)um-1n-2 umn+um+1nh2+(1-θ)umn(1-umn)+θumn+1(1-umi+1)+θJmi+1(umi+1-umn)

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

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

u(xm-1)=um-1,  u(xm)=um,  u(xm+1)=um+1

and yield exponential solutions in (xm−1, xm+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

(θJmn-1k) umn+1+θum-1n+1-2 umn+1+um+1n+1h2=(θJmn-1k) umn-(1-θ)um-1n-2 umn+um+1nh2-umn(1-umn).

Let (θJmn-1k)umn+1+θum-1n+1-2umn+1+um+1n+1h2=zmn+1. Equation (103) becomes

zmn+1=(θJmn-1k)  umn-(1-θ)um-1n-2 umn+um+1nh2-umn(1-umn)

which accounts for diffusion, reaction, and transient effects in the differential operator [25]. The solution of Equation (104) depends on the sign of Dmn≡θJmn-1k. We have the following three cases:

a) If Dmn=0, the solution of Equation (104) subject to the condition (102) gives the finite difference expression zmn+1=um-1n-2 umn+um+1nh2.

b) If Dmn=-θ(λm)2<0, the solution of Equation (104) subject to the condition (102) gives the three-point finite difference expression zmn+1=Dmn2(um-1n-2cosh(λmh)umn+um+1n1-cosh(λmh)).

c) If Dmn=θ(λm)2>0, the solution of Equation (104) subject to the condition (102) gives the three-point finite difference expression zmn+1=Dmn2(um-1n-2cos(λmh)umn+um+1n1-cos(λmh)).

Remark 1. We can make the following remarks

a) The values θ=12,1, corresponding to the time linearization methods.

b) The quasilinear full exponential corresponds to the iterative solution of zmi+1=-(1-θ)um-1n-2 umn+um+1nh2-(1-θ)umn(1-umn)-θ umi(1-umi)+θ Jmiumi-umnk

for i = 1, 2, ··, n − 1, have analogous solutions to those reported in time-linearized full exponential techniques section and the cases θ=12 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.

θJmnumn+1+θum-1n+1-2 umn+1+um+1n+1h2=θJmnumn-(1-θ)um-1n-2 umn+um+1nh2+um-umnk-umn(1-umn)

Let θJmnumn+1+θum-1n+1-2umn+1+um+1n+1h2=zmn+1. The Equation (109) becomes

zmn+1=θJmnumn-(1-θ)um-1n-2 umn+um+1nh2+um-umnk-umn(1-umn)

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

a) If Jmn=0 um-1-(2+h2k)um+um+1=-h2k(umn+kumn(1-umn)),  tn<t<tn+1.

b) If Jmn=-(λm)2<0                                   um-1-2+(kJmn-1)(e-λmh+eλmh)kJmnum                                  +um+1                                  =2-(e-λmh+eλmh)kJmn                                  [-um+k(Jmnum-umn(1-umn))],tn<t<tn+1.

c) If Jmn=(λm)2>0

                             um-1-21+(kJmn-1)cos(λmh)kJmnum+um+1                            =21-(cos(λmh)kJmn                            [-um+k(Jmnum-umn(1-umn))],tn<t<tn+1.

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

θJmiumn+1+θum-1n+1-2 umn+1+um+1n+1h2=θJmnumn-(1-θ)um-1n-2 umn+um+1nh2+um-umnk-(1-θ)umn(1-umn)-θumi(1-umi)

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≪1θJin 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 zmn+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

(emn)2=1N∑j=1N(uji+1-uji)2≤10-|α|

where i = 1, 2···n and |α| is an integer obtained from numerical computation, N denotes the number of grid points, and emn is the error which is defined by emn=umn-vmn.

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

 8. Numerical Results

The stability region of FTCS is k≤h22. For h = 0.01, we obtain k ≤ 5 × 10−5. In the case of NSFD, the condition for positivity gives ϕ(k)≤h22 where ϕ(k)=1-eλkλ. We tabulate L1 and L∞ errors at certain values of k using ρ = 1, h = 0.01, at time, Tmax=2.5×10-3 using FTCS, NSFD, and EEFD schemes at certain various values of time-step size as chosen from Tmax/52, Tmax/100, Tmax/200, ⋯ , Tmax/1,800, Tmax/1,900 and Tmax/2,000. The errors are shown in Table 1.

L1 and L∞ errors at some different values of time-step size, k with ρ = 1 at time, Tmax=2.5×10-3 with spatial mesh size, h = 0.01 using three methods.

Values of k FTCS NSFD EEFD
L1 error L∞error L1 error L∞error L1 error L∞error
Tmax/52 5.4192 × 10−9 6.3334 × 10−9 3.2264 × 10−8 3.7216 × 10−8 4.1091 × 10−9 5.1492 × 10−9
Tmax/100 2.8854 × 10−9 3.3694 × 10−9 1.7026 × 10−8 1.9624 × 10−8 2.0683 × 10−9 2.6091 × 10−9
Tmax/200 1.5135 × 10−9 1.7645 × 10−9 8.7757 × 10−9 1.0098 × 10−8 9.6251 × 10−10 1.2354 × 10−9
Tmax/300 1.0563 × 10−9 1.7645 × 10−9 6.0260 × 10−9 6.0260 × 10−9 5.9431 × 10−10 7.7761 × 10−10
Tmax/400 8.2780 × 10−10 9.6234 × 10−10 4.6512 × 10−9 5.3361 × 10−9 4.1024 × 10−10 5.4763 × 10−10
Tmax/500 6.9066 × 10−10 8.0194 × 10−10 3.8264 × 10−9 4.3837 × 10−9 2.9973 × 10−10 4.1032 × 10−10
Tmax/600 5.9924 × 10−10 6.9502 × 10−10 3.2766 × 10−9 3.7488 × 10−9 2.2612 × 10−10 3.1904 × 10−10
Tmax/700 5.3392 × 10−10 6.1863 × 10−10 2.8838 × 10−9 3.2954 × 10−9 1.7352 × 10−10 2.5366 × 10−10
Tmax/800 4.8495 × 10−10 5.6136 × 10−10 2.5894 × 10−9 2.9554 × 10−9 1.3402 × 10−10 2.0476 × 10−10
Tmax/900 4.4686 × 10−10 5.1681 × 10−10 2.3602 × 10−9 2.6910 × 10−9 1.0334 × 10−10 1.6694 × 10−10
Tmax/1,000 4.1639 × 10−10 4.8117 × 10−10 2.1769 × 10−9 2.4795 × 10−9 7.8784 × 10−11 1.3668 × 10−10
Tmax/1,100 3.9144 × 10−10 4.5199 × 10−10 2.0269 × 10−9 2.3064 × 10−9 5.8705 × 10−11 1.1192 × 10−10
Tmax/1,200 3.7066 × 10−10 4.2769 × 10−10 1.9019 × 10−9 2.1622 × 10−9 4.1956 × 10−11 9.1419 × 10−11
Tmax/1,300 3.5308 × 10−10 4.0714 × 10−10 1.7962 × 10−9 2.0402 × 10−9 3.0241 × 10−11 7.4118 × 10−11
Tmax/1,400 3.3801 × 10−10 3.8953 × 10−10 1.7056 × 10−9 1.9356 × 10−9 2.3908 × 10−11 5.9401 × 10−11
Tmax/1,500 3.2496 × 10−10 3.7427 × 10−10 1.6270 × 10−9 1.8449 × 10−9 2.1174 × 10−11 4.6738 × 10−11
Tmax/1,600 3.1353 × 10−10 3.6091 × 10−10 1.5583 × 10−9 1.7656 × 10−9 2.0904 × 10−11 3.5802 × 10−11
Tmax/1,700 3.0345 × 10−10 3.4913 × 10−10 1.4976 × 10−9 1.6956 × 10−9 2.2344 × 10−11 4.2194 × 10−11
Tmax/1,800 2.9449 × 10−10 3.3866 × 10−10 1.4438 × 10−9 1.6335 × 10−10 2.4984 × 10−11 4.8536 × 10−11
Tmax/2,000 2.9449 × 10−10 3.3866 × 10−10 1.3521 × 10−9 1.5279 × 10−10 3.2504 × 10−11 5.9412 × 10−11

Bold values indicate the lowest errors.

In the case of FTCS, NSFD, minimum L1, and L∞ errors occur at k = Tmax/2,000 while in the case of EEFD, the errors are least k = Tmax/1,600. The least error is of order 10−9 and 10−10 in the case of NSFD and FTCS respectively while the least error is of order 10−11 in the case of EEFD. EEFD is a better scheme than FTCS at all values of k used.

We obtain values for L1 and L∞ errors at certain values of k using ρ = 102, h = 0.01, at time, Tmax=2.5×10-3 using the three methods in Table 2. The schemes behave differently. In all the three methods FTCS, NSFD, and EEFD, the L1 and L∞ errors keep on decreasing as the values of k are decreased gradually from k = Tmax/52 to k = Tmax/2,000.

L1 and L∞ errors at some different values of time-step size, k with ρ = 102 at time, Tmax=2.5×10-3 with spatial mesh size, h = 0.01 using three methods.

Values of k FTCS NSFD EEFD
L1 error L∞error L1 error L∞error L1 error L∞error
Tmax/52 3.2613 × 10−5 6.9142 × 10−5 6.3178 × 10−5 1.3471 × 10−4 2.2654 × 10−5 7.2263 × 10−5
Tmax/100 1.7120 × 10−5 3.6652 × 10−5 3.3705 × 10−5 7.2133 × 10−5 1.1628 × 10−5 3.7401 × 10−5
Tmax/200 8.7183 × 10−6 1.9060 × 10−5 1.7709 × 10−5 3.8414 × 10−5 5.6602 × 10−6 1.8584 × 10−5
Tmax/300 5.9159 × 10−6 1.3196 × 10−5 1.2371 × 10−5 2.7303 × 10−5 3.6732 × 10−6 1.2324 × 10−5
Tmax/400 4.5154 × 10−6 1.0268 × 10−5 9.7020 × 10−6 2.1809 × 10−5 2.6794 × 10−6 9.1874 × 10−6
Tmax/500 3.6749 × 10−6 8.5169 × 10−6 8.1001 × 10−6 1.8545 × 10−5 2.0833 × 10−6 7.3251 × 10−6
Tmax/600 3.1146 × 10−6 7.3489 × 10−6 7.0319 × 10−6 1.6379 × 10−5 1.6864 × 10−6 6.0864 × 10−6
Tmax/700 2.7143 × 10−6 6.5146 × 10−6 6.2689 × 10−6 1.4842 × 10−5 1.4024 × 10−6 5.2014 × 10−6
Tmax/800 2.4141 × 10−6 5.8925 × 10−6 5.6966 × 10−6 1.3689 × 10−5 1.1903 × 10−6 4.5376 × 10−6
Tmax/900 2.1806 × 10−6 5.4093 × 10−6 5.2515 × 10−6 1.2802 × 10−5 1.0244 × 10−6 4.0219 × 10−6
Tmax/1,000 1.9941 × 10−6 5.0227 × 10−6 4.8954 × 10−6 1.2093 × 10−5 8.9174 × 10−7 3.6143 × 10−6
Tmax/1,100 1.8416 × 10−6 4.7064 × 10−6 4.6040 × 10−6 1.1513 × 10−5 7.8356 × 10−7 3.2846 × 10−6
Tmax/1,200 1.7144 × 10−6 4.4428 × 10−6 4.3612 × 10−6 1.1030 × 10−5 7.0876 × 10−7 3.0099 × 10−6
Tmax/1,300 1.6068 × 10−6 4.2197 × 10−6 4.1557 × 10−6 1.0621 × 10−5 6.5558 × 10−7 2.7774 × 10−6
Tmax/1,400 1.5153 × 10−6 4.0298 × 10−6 3.9797 × 10−6 1.0274 × 10−5 6.1557 × 10−7 2.5767 × 10−6
Tmax/1,500 1.4371 × 10−6 3.8658 × 10−6 3.8279 × 10−6 9.9737 × 10−6 5.8466 × 10−7 2.4043 × 10−6
Tmax/1,600 1.3698 × 10−6 3.7222 × 10−6 3.6960 × 10−6 9.7113 × 10−6 5.6048 × 10−7 2.2528 × 10−6
Tmax/1,700 1.3113 × 10−6 3.5956 × 10−6 3.5803 × 10−6 9.4798 × 10−6 5.4168 × 10−7 2.1186 × 10−6
Tmax/1,800 1.2601 × 10−6 3.4830 × 10−6 3.4781 × 10−6 9.2740 × 10−6 5.2654 × 10−7 2.0034 × 10−6
Tmax/2,000 1.1748 × 10−6 3.2916 × 10−6 3.3055 × 10−6 8.9242 × 10−6 5.0458 × 10−7 1.807 × 10−6

Bold values indicate the lowest errors.

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

L1 and L∞ errors at some different values of time-step size, k with ρ = 104 at time, Tmax=2.5×10-3 with spatial mesh size, h = 0.01 using three methods.

Values of k FTCS NSFD EEFD
L1 error L∞error L1 error L∞error L1 error L∞error
Tmax/52 3.6440 × 10−1 1.6012 1.1960 × 10−1 9.1526 × 10−1 over flow over flow
Tmax/100 6.5143 × 10−2 6.9877 × 10−1 7.0119 × 10−2 7.1253 × 10−1 1.3594 × 10−2 1.5978 × 10−1
Tmax/200 3.2980 × 10−2 3.9927 × 10−1 3.6170 × 10−2 4.1776 × 10−1 1.0084 × 10−2 1.1786 × 10−1
Tmax/300 2.0726 × 10−2 2.5613 × 10−1 2.3034 × 10−2 2.7854 × 10−1 8.9264 × 10−3 1.0476 × 10−1
Tmax/400 1.4262 × 10−2 1.7830 × 10−1 1.6098 × 10−2 1.9602 × 10−1 8.3519 × 10−3 9.8394 × 10−2
Tmax/500 1.0268 × 10−2 1.2943 × 10−1 1.1810 × 10−2 1.4691 × 10−1 8.0087 × 10−3 9.4591 × 10−2
Tmax/600 7.5556 × 10−3 9.5604 × 10−2 8.8965 × 10−3 1.1270 × 10−1 7.7810 × 10−3 9.2071 × 10−2
Tmax/700 5.5928 × 10−3 7.1037 × 10−2 6.7882 × 10−3 8.7558 × 10−2 7.6182 × 10−3 9.0294 × 10−2
Tmax/800 4.1067 × 10−3 5.2473 × 10−2 5.1924 × 10−3 6.8387 × 10−2 7.4974 × 10−3 8.8949 × 10−2
Tmax/900 2.9425 × 10−3 3.7988 × 10−2 3.9442 × 10−3 5.3326 × 10−2 7.4020 × 10−3 8.7921 × 10−2
Tmax/1, 000 2.0059 × 10−3 2.6389 × 10−2 2.9442 × 10−3 4.1202 × 10−2 7.3264 × 10−3 8.7094 × 10−2
Tmax/1,100 1.2403 × 10−3 1.6901 × 10−2 2.1334 × 10−3 3.1242 × 10−2 7.2643 × 10−3 8.6421 × 10−2
Tmax/1,200 6.2582 × 10−4 9.0012 × 10−3 1.4824 × 10−3 2.2920 × 10−2 7.2120 × 10−3 8.5864 × 10−2
Tmax/1,300 1.9347 × 10−4 2.3239 × 10−3 9.7801 × 10−4 1.5866 × 10−2 7.1694 × 10−3 8.5384 × 10−2
Tmax/1,400 4.2509 × 10−4 4.3457 × 10−3 6.3252 × 10−4 9.8126 × 10−3 7.1310 × 10−3 8.4982 × 10−2
Tmax/1,500 8.3336 × 10−4 8.7691 × 10−3 4.5235 × 10−4 4.5623 × 10−3 7.0992 × 10−3 8.4631 × 10−2
Tmax/1,600 1.1913 × 10−3 1.3138 × 10−2 5.0611 × 10−4 5.3708 × 10−3 7.0703 × 10−3 8.4321 × 10−2
Tmax/1,700 1.5077 × 10−3 1.7003 × 10−2 8.4385 × 10−4 7.6874 × 10−3 7.0451 × 10−3 8.4052 × 10−2
Tmax/1,800 1.7895 × 10−3 2.0446 × 10−2 1.1467 × 10−3 1.0154 × 10−2 7.0231 × 10−3 8.3812 × 10−2
Tmax/2,000 1.5538 × 10−2 1.8402 × 10−1 1.6628 × 10−3 1.5383 × 10−2 6.9852 × 10−3 8.3399 × 10−2

Bold values indicate the lowest errors.

We obtain plots of numerical solution vs. x at time, Tmax=2.5×10-3 using three methods FTCS, NSFD, and EEFD in Figure 1.

Plot of initial solution, numerical solution of (FTCS), NSFD, and EEFD respectively in (A–C) and exact solution against x at time, Tmax=2.5×10-3 using h = 0.01 and optimal value of time step size of k for three different values of ρ namely 1, 102, and 104.

 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, 102, and 104. 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 L1 and L∞ errors from the three methods is very small which lead to conclude that the best methods are FTCS and NSFD. Also for ρ = 102 and ρ = 104, 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.

Rate of convergence in time for FTCS using different values of ρ at time, Tmax=2.5×10-3 with h = 0.01.

Values of ρ Values of k Rate of convergence in time using L1 Rate of convergence in time using L∞error
1 Tmax/100
Tmax/200 0.9309 0.9332
Tmax/400 0.8705 0.8746
Tmax/800 0.7714 0.7776
102 Tmax/100
Tmax/200 0.9736 0.9433
Tmax/400 0.9492 0.8924
Tmax/800 0.9034 0.8012
104 Tmax/100
Tmax/200 0.9820 0.8075
Tmax/400 1.2094 1.1631

Rate of convergence in time for NSFD using different values of ρ at time, Tmax=2.5×10-3 with h = 0.01.

Values of ρ Values of k Rate of convergence in time using L1 Rate of convergence in time using L∞error
1 Tmax/100
Tmax/200 0.9561 0.9585
Tmax/400 0.9159 0.9202
102 Tmax/100
Tmax/200 0.9284 0.9090
Tmax/400 0.8681 0.8167
104 Tmax/100
Tmax/200 0.9769 0.7702

Rate of convergence in time for EEFD using different values of ρ at time, Tmax=2.5×10-3 with h = 0.01.

Values of ρ Values of k Rate of convergence in time using L1 Rate of convergence in time using L∞error
1 Tmax/100
Tmax/200 1.1030 1.0795
Tmax/400 1.2306 1.1727
102 Tmax/100
Tmax/200 1.0384 1.0091
Tmax/400 1.0789 1.0162
Tmax/800 1.1714 1.0176
104 Tmax/100
Tmax/200 0.4317 0.4388
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.

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 Gardner CS Greene JM Kruskal MD Miura RM. Method for solving the Korteweg-de Vries equation. Front Neurosci Phys Rev Lett. (1967) 19:1095. 10.1103/PhysRevLett.19.1095 Weiss J Tabor M Carnevale G. The Painlevé property for partial differential equations. J Math Phys. (1983) 155:522–6. 10.1063/1.525721 Kudryashov NA. On types of nonlinear nonintegrable equations with exact solutions. Phys Lett A. (1991) 24:269–75. 10.1016/0375-9601(91)90481-M Lou SY Huang GX Ruan HY. Exact solitary waves in a convecting fluid. J Phys A Math Gen. (1991) 24:L587. 10.1088/0305-4470/24/11/003 Kudryashov NA. Exact solutions of a family of Fisher equations. Theor Math Phys. (1993) 94:211–218. 10.1007/BF01019332 Kudryashov NA. Numerical convective schemes based on accurate computation of space derivatives. J Comput Phys. (1973) 13:100–13. 10.1016/0021-9991(73)90128-9 Cooley JW Lewis P Welch P. The finite Fourier transform. IEEE Trans Audio Electroacoust. (1969) 17:77–85. 10.1109/TAU.1969.1162036 Gazdag J Canosa J. Numerical solution of Fisher's equation. J Appl Probab. (1974) 11:445–57. 10.2307/3212689 Fisher RA. The wave of advance of advantageous genes. Ann Hum Genet. (1937) 17:355–69. 10.1111/j.1469-1809.1937.tb02153.x Moran PAP. The Statistical Processes of Evolutionary Theory. Oxford: Clarendon Press; Oxford University Press (1962). Kendall DG. A Form of Wave Propagation Associated With the Equation of Heat Conduction. Cambridge: Cambridge University Press (1948). 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. Canosa J. On a nonlinear diffusion equation describing population growth. IBM J Res Dev. (1937) 17:307–13. 10.1147/rd.174.0307 Jeffrey A Kakutani T. Weak nonlinear dispersive waves: a discussion centered around the Korteweg-de Vries equation. Siam Rev SIAM. (1972) 14:582–643. 10.1137/1014101 Kudryashov NA Zakharchenko AS. A note on solutions of the generalized Fisher equation. Appl Math Lett. (2014) 32:53–6. 10.1016/j.aml.2014.02.009 Newman TJ Kolomeisky EB Antonovics J. Population dynamics with global regulation: the conserved Fisher equation. Phys Rev Lett. (2004) 92:228103. 10.1103/PhysRevLett.92.22810315245262 Chen BM Kojouharov HV. Non-standard numerical methods applied to subsurface biobarrier formation models in porous media. Bull Math Biol. (1999) 61:779–98. 10.1006/bulm.1999.011317883224 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. 10.1016/j.ecocom.2017.06.004 Wang XY. Exact and explicit solitary wave solutions for the generalised Fisher equation. Phys Lett A. (1988) 131:277–9. 10.1016/0375-9601(88)90027-8 Chandraker V Awasthi A Jayaraj S. A numerical treatment of Fisher equation. Procedia Eng. (2015) 131:1256–62. 10.1016/j.proeng.2015.11.481 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. 10.1016/j.amc.2003.10.050 Moghimi M Hejazi FSA. Variational iteration method for solving generalized Burger-Fisher and Burger equations. Chaos Solitons Fractals. (2007) 33:1756–61. 10.1016/j.chaos.2006.03.031 Fahmy ES. Travelling wave solutions for some time-delayed equations through factorizations. Chaos Solitons Fractals. (2008) 38:1209–16. 10.1016/j.chaos.2007.02.007 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. 10.1002/(SICI)1098-2426(199701)13:1<51::AID-NUM4>3.0.CO;2-L Ramos JI. Exponential methods for one-dimensional reaction-diffusion equations. Math Methods Appl Sci. (2005) 170:380–98. 10.1016/j.amc.2004.12.003 Anguelov R Kama P Lubuma JMS. On non-standard finite difference models of reaction-diffusion equations. J Comput Appl Math. (2005) 175:11–29. 10.1016/j.cam.2004.06.002 Hagstrom T Keller HB. The numerical calculation of traveling wave solutions of nonlinear parabolic equations. SIAM J Scientific Stat Comput. (1986) 7:978–88. 10.1137/0907065 Li S Petzold L Ren Y. Stability of moving mesh systems of partial differential equations. SIAM J Scientific Comput. (1998) 20:719–38. 10.1137/S1064827596302011 Kudryashov NA. One method for finding exact solutions of nonlinear differential equations. Commun Nonlinear Sci Numer Simulat. (2012) 17:2248–53. 10.1016/j.cnsns.2011.10.01635730291 Polyanim AD Zaitsev ZF. Handbook of Nonlinear Partial Differential Equations. 2nd ed. Chapman and Hall/CRC Press. (2011). p. 1912. Qiu Y Sloan DM. Numerical solution of Fisher's equation using a moving mesh method. J Comput Phys. (1998) 146:726–746. 10.1006/jcph.1998.6081 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). Sutton OJ. Long time L∞(L2) a posteriori error estimates for fully discrete parabolic problems. arXiv[Preprint].arXiv:180303207. (2018). 10.1093/imanum/dry078 Chen-Charpentier BM Kojouharov HV. An unconditionally positivity preserving scheme for advection-diffusion reaction equations. Math Comput Model. (2013) 57:2177–85. 10.1016/j.mcm.2011.05.005 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. 10.1186/s13662-019-2080-x Durran DR. Numerical Methods for Fluid Dynamics (With Applications to Geophysics). Springer Science & Business Media (2010). 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. 10.1016/0021-9991(84)90004-4 Zabusky NJ Kruskal MD. Interaction of solitons in a collisionless plasma and the recurrence of initial states. Phys Rev Lett. (1965) 15:240. 10.1103/PhysRevLett.15.240 Mickens RE. Dynamic consistency: a fundamental principle for constructing nonstandard finite difference schemes for differential equations. J Diff Equat Appl. (2005) 11:645–53. 10.1080/10236190412331334527 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. 10.1080/10236198.2019.1689236 Bhattacharya MC. An explicit conditionally stable finite difference equation for heat conduction problems. Int J Numer Methods Eng. (1985) 21:239–65. 10.1002/nme.1620210205 Handschuh RF Keith TG. Applications of an exponential finite-difference technique. Numer Heat Transfer A Appl. (1992) 22:363–78. 10.1080/1040778920894477320715071 Ĩ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. Ĩnan B Bahadır AR. Finite difference methods for the generalized Huxley and Burgers-Huxley equations. Kuwait J Sci. (2017) 44:20–7. 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. 10.1080/00207160.2018.1478416 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. 10.1002/mma.6067 Lubuma J Roux A. An improved theta-method for systems of ordinary differential equations. J Diff Equat Appl. (2003) 9:1023–35. 10.1080/1023619031000146904