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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

We have the following solutions

We finally get as exact solution

where V₁ and V₂ stand for velocity.

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

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

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

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

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

The exact solution is given by

The boundary conditions are as follows

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

The rate of convergence with respect to time is defined by

5. Forward in Time Central Space

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

which leads to

or

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

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

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

which can be rewritten as

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

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

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

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

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

which can be written as

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

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

Taylor series expansion about (n, m) gives

which can be rewritten as

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

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

It follows that

We have

which can be rewritten as

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

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

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

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

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

For n = 2, we have

For n, we have

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

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

where

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]

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

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

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

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

and

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

for Equation (22) is

where

The Equation (66) gives the following single expression

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

It follows that

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

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

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

We also have from Equation (71)

Based on the symmetric condition, we can take

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

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

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

From Equation (78), we have

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

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

Hence,

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

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

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

and

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

Proof.

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

This gives

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

It follows that

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

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

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

It follows

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

$\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

Let

We have

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

For n = 2, we have

By recurrence for n, we have

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

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

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

which yields a non-iterative time linearization method.

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

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:

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

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

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

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

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

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

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.

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

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$

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

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

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

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

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.

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.

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⁻².

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.

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.

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