Research Article  Open Access
Sangkwon Kim, Chaeyoung Lee, Hyun Geun Lee, Hyundong Kim, Soobin Kwak, Youngjin Hwang, Seungyoon Kang, Seokjun Ham, Junseok Kim, "An Unconditionally Stable PositivityPreserving Scheme for the OneDimensional Fisher–Kolmogorov–Petrovsky–Piskunov Equation", Discrete Dynamics in Nature and Society, vol. 2021, Article ID 7300471, 11 pages, 2021. https://doi.org/10.1155/2021/7300471
An Unconditionally Stable PositivityPreserving Scheme for the OneDimensional Fisher–Kolmogorov–Petrovsky–Piskunov Equation
Abstract
In this study, we present an unconditionally stable positivitypreserving numerical method for the Fisher–Kolmogorov–Petrovsky–Piskunov (Fisher–KPP) equation in the onedimensional space. The Fisher–KPP equation is a reactiondiffusion system that can be used to model population growth and wave propagation. The proposed method is based on the operator splitting method and an interpolation method. We perform several characteristic numerical experiments. The computational results demonstrate the unconditional stability, boundedness, and positivitypreserving properties of the proposed scheme.
1. Introduction
In an infinitely long domain, the behavior of virile mutant propagation can be modeled by Fisher’s equation [1]. It was also independently studied by Kolmogorov, Petrovsky, and Piskunov [2]. There are some numerical methods to solve the Fisher–Kolmogorov–Petrovsky–Piskunov (Fisher–KPP) or fractional Fisher–KPP equations such as the spectral collocationtype method [3, 4], Adams–Bashforth–Moulton method [5], and MittagLeffler kernel [6]. In addition, various finite difference methods [7–10] and the radial basis functions [11] have been developed.
In this paper, we present an unconditionally stable positivitypreserving numerical method for the Fisher–KPP equation:where is the population density at position and time . Here, is the diffusion coefficient, is the coefficient of nonlinearity, and is a positive integer. There exist several other unique approaches to solve the Fisher–KPP equation. Saeed and Mustafa [12] used the reduced differential transformation method to solve the equation and showed that their result coincides with the exact solution. From a biological point of view, the population density should be nonnegative and bounded from above by one if the density is scaled by the carrying capacity. ElHachem et al. [13] studied approaches to modelling biological invasion and recession using the Fisher–KPP model. In addition, they utilized the Fisher–Stefan model, another generalization of the Fisher–KPPtype model that simulates biological recession. Xu et al. [14] proposed a nonlinear finite volume scheme for twodimensional diffusion equations which preserves the maximum principle. Fourthorder Fisher–KPP equation is one of the widely studied and extended variations, and various methods have been proposed to solve it. Kadri and Omrani [15] utilized the extrapolation technique, and Başhan et al. [16] applied the differential quadrature method. An explicit positivitypreserving numerical method for the Fisher–KPP equation has been proposed [7]. Furthermore, the positivitypreserving schemes for the Fisher–KPP equation were developed using continuous finite element discretizations [17], and discontinuous Galerkin schemes were proved to be entropystable [18]. Ilati and Dehghan [19] proposed the local boundary integral equation method to solve the extended Fisher–KPP equation. The presented method is based on generalized moving least squares, which can be used to directly approximate local weak forms. Oruç [20] developed an efficient wavelet collocation method for the nonlinear 2D Fisher–KPP equation and extended Fisher–Kolmogorov equation. Castillo and Gómez [21] used the local discontinuous Galerkin method with the symmetric Strang operator splitting scheme to solve Fisher–KPP. Kabir [22] considered a model with travelling wave solutions with a qualitatively similar structure to that observed in the Fisher–KPP equation and performed numerical simulations. McCue et al. [23] studied exact sharpfronted travelling wave solutions of the Fisher–KPP equation to the Fisher–Stefantype moving boundary problem. They demonstrated large time evolutionary dynamics of the exact travelling wave solution and considered these travelling waves for cell migration and proliferation mathematical models. Recently, numerical research studies for the nonlocal Fisher–KPP equation have been actively conducted to validate the asymptotic behavior [24], positivity of the solution [25, 26], and stability of the travelling wave solution [27, 28].
However, if in equation (1) is large, then the Fisher–KPP equation becomes a very stiff problem to be numerically solved. To resolve this stiff problem, we propose the operator splitting method (OSM) with a recently developed interpolation method [29]. The computational results demonstrate the unconditional stability and positivitypreserving properties of the proposed scheme.
The contents of this article are as follows. In Section 2, we present the proposed numerical scheme for the Fisher–KPP equation. In Section 3, we conduct computational experiments to verify the robustness of the proposed scheme. In Section 4, conclusions are drawn.
2. Numerical Solution Algorithm
To efficiently solve equation (1), we apply a recently developed interpolation method [29]. Now, we present the numerical solution scheme of the Fisher–KPP equation in onedimensional space . Let be an integer and be the uniform spatial step size. Let be the numerical approximation of , where for and . Using the OSM, we divide equation (1) into two parts: diffusion and nonlinear equations. First, we solve the diffusion equation using the fully implicit Euler finite difference method (FDM):
Equation (2) can be discretized using the FDM as
Equation (3) is reorganized as
For the boundary condition, we use the zero Neumann boundary condition: and for all values of . The Dirichlet boundary condition can also be used. Then, we solve equation (4) using the Thomas algorithm [30]. Second, we solve the nonlinear equation using an interpolation method:
For simplicity of exposition, let ; then, equation (5) becomes
To solve equation (6), we use a recently developed interpolation method [29]. Let us assume that the time step is given. We uniformly discretize the interval , i.e., . Note that is a positive integer. We solve equation (6) on with a smaller time step until . Here, will be defined later. Figure 1 shows a simulation result with the initial condition for . That is, for ,
Once we compute all for , then we will use these values when we numerically solve equation (6) by using an interpolation.
Next, we find a value of appropriate to stably integrate equation (7) for given . If or , then from equation (7), we have
Next, when and , we want for . Here, to show , by mathematical induction, it is sufficient to show . Then, let us suppose that , , and ; from (7), we have
As shown in Figure 2, for ; therefore, we have from equation (9).
Because is always satisfied, we only need to find satisfying the following condition:
The righthand side of inequality (10) is decreasing with respect to ; see Figure 3.
Therefore, we consider the limit of the term as approaches 1. By using L’Hopital’s rule, we have
From the definition of and stability condition (11), should be satisfied. Therefore, unless otherwise specified, we set , where is the greatest integer less than or equal to . Given the intermediate solution from the first step, , the second step is to solve equation (6) by using the precomputed values and interpolation. We can find an interval which satisfies for some and . Then, we definewhich is schematically illustrated in Figure 4.
We should note that, in the solution algorithm, the first step is a diffusion equation solver, which is unconditionally stable and satisfies the discrete maximum principle. The second step is an interpolation, which is unconditionally stable, monotonically increasing, and bounded by one. Because two steps are both unconditionally stable and positive preserving, therefore, the full step is unconditionally stable, positive preserving, and bounded by one.
3. Numerical Experiments
In this section, we conduct various computational experiments such as the effect of the model and numerical parameters of the proposed algorithm. In order to demonstrate the performance of the proposed algorithm, we conduct comparative experiments with the previous work.
3.1. Convergence Tests
Let us consider an initial condition for equation (1) on the computational domain with and :
The exact solution is given in [31]:
We run the simulation up to terminal time with different time steps: , 0.2, 0.1, and 0.05. The other numerical parameter values used are and . Figure 5 shows the numerical results and the exact solution (14) with . We can observe that the numerical solutions converge to the exact solution as we refine the time step.
Table 1 shows the norm error and convergence rate of numerical solutions with the same conditions above except for the final time and the fixed , where is the number of time iterations of the simulation for each . We can observe that the numerical scheme is firstorder accurate in time.

Next, we investigate the accuracy of numerical results on the value of . Figure 6 shows the convergence of numerical solutions with respect to values on subdomain . Here, we used the same parameter values as before except values and .
Last, we study the convergence of numerical results with respect to . The applied conditions are the same as the time regarding simulation, except for . With all the other parameters fixed, the spatial step size is refined to calculate the convergence rate. norm error and convergence rate are shown in Table 2, and it can be observed that the numerical scheme is secondorder accurate in space.

Time convergence rate and spatial convergence rate are clearly different: first order and second order. This is due to the fact that the interpolation method is second order in both time and space, yet the fully implicit Euler method is first order in time and second order in space. Therefore, after combining the two methods using the OSM, time convergence rate is one, while space convergence rate is two.
3.2. Stability Test
We demonstrate the stability of the proposed scheme with a numerical simulation with initial condition (13) for (1) on the computational domain with , , , and . We run the simulation up to terminal time with different time steps: , 0.3, and 3. Figure 7 shows the numerical results and the exact solution (14) with . It can be observed that the numerical solutions do not blow up with large time step. Therefore, the proposed scheme is stable. However, when a large time system is used, the accuracy of the numerical solution is low.
3.3. Effect of
Now, let us consider the effect of on the evolution dynamics. Let us consider equation (6). By using the explicit Euler’s method, we present the evolution of with the initial condition on various values of . Here, we use the time step and the terminal time :
Figure 8 shows the evolution of with , and 100. As the parameter value increases, the solution evolves faster.
However, as shown in Figure 8, there is little difference between and . Next, we investigate the numerical solution of equation (9) on as , , and ; see Figure 9.
To clearly identify this effect, numerical experiments are conducted using initial condition . We set the computational domain , , , terminal time , , , and . Figure 10 shows the numerical solutions with various values of , and 100.
The evolution with large is faster than that with small , especially near . However, evolution with in the vicinity of and has a small difference.
Next, let us set an initial condition equation (13) on the computational domain with and . We run the simulation up to time with , , and . Figure 11 shows the numerical results with different values on subdomain .
(a)
(b)
(c)
Figure 12 shows the numerical solution with an initial condition if ; otherwise, for different values. Here, we used the same parameter values as before. We can observe that the numerical solutions become stiffer as the value increases.
3.4. Positive Preserving
We conduct numerical experiment using the initial conditionwhere to confirm that the algorithm satisfies positive preserving and boundedness by one. The parameters used in the numerical experiment are as follows: , , , , , and . Figures 13(a) and 13(b) show the evolution of numerical solutions for and at . Figures 13(c) and 13(d) show the temporal evolution of the maximum and minimum of over time when and , respectively. We can confirm that is satisfied regardless of the value. The proposed algorithm is a positivepreserving scheme.
(a)
(b)
(c)
(d)
3.5. Comparison with the Previous Work
To confirm the performance of our proposed method, we compare it with the recently developed method. Izadi [10] proposed a secondorder accurate and unconditionally stable FDM for the classical Fisher–KPP equation:where and with and . Let us consider an initial condition on : if ; otherwise, . We run the simulation up to time with , , and . Figure 14 shows the numerical solutions with various values of , and 800. When , the numerical solutions for both methods converge to the exact solution. For the large value , the proposed method obtained a stable numerical solution, but the Izadi method [10] obtained an unstable numerical solution; see Figures 14(b) and 14(c).
(a)
(b)
(c)
3.6. SecondOrder Accurate Scheme
In this section, we present a secondorder accurate scheme in time. Previously, the proposed scheme, which is a firstorder accurate scheme in time, separates the governing equation into two parts over a time step: the diffusion equation and nonlinear equation . That is, the numerical solution is with time step size . For a secondorder accurate scheme, meanwhile, . Note that all numerical schemes used for and should be at least secondorder accurate in time. Please refer to [29] for more details. Here, we solve the diffusion equation using the Crank–Nicolson scheme and use a sufficiently large for the nonlinear equation to guarantee the secondorder accuracy in time. To verify the convergence of the numerical scheme, we set the same initial condition (13) and parameters as in Section 3.1. Figure 15 illustrates that the numerical solutions converge to the exact solution as the time step is refined. In addition, Table 3 indicates that the entire scheme is secondorder accurate in time.

We tested the spatial convergence for the aforementioned scheme. It can be assumed that it is secondorder accurate in space since both Crank–Nicolson method and interpolation method are second order. To confirm the convergence, the same condition as the spatial convergence test in Section 3.1 is used. The derived results are given in Table 4, and it shows that the method is secondorder accurate in space.

4. Conclusions
In this study, we presented an unconditionally stable positivitypreserving numerical method for the Fisher–KPP equation. The proposed method is based on the OSM and an interpolation method. The computational results demonstrated the unconditional stability and positivitypreserving properties of the proposed scheme. In the future work, we will extend the proposed scheme to multidimensional spaces [20].
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The corresponding author (J. S. Kim) was supported by the Brain Korea 21 FOUR from the Ministry of Education of Korea.
References
 R. A. Fisher, “The wave of advance of advantageous genes,” Annals of eugenics, vol. 7, pp. 355–369, 1936. View at: Publisher Site  Google Scholar
 A. Kolmogorov, N. Petrovsky, and S. Piscounov, “Étude de l’équations de la diffusion avec croissance de la quantitaté de matiére et son application a un probléme biologique,” Bull Univ Moskow, Ser Internat, Sec A., vol. 1, pp. 1–25, 1937. View at: Google Scholar
 M. M. Khader and K. M. Saad, “A numerical approach for solving the fractional Fisher equation using Chebyshev spectral collocation method,” Chaos, Solitons & Fractals, vol. 110, pp. 169–177, 2018. View at: Publisher Site  Google Scholar
 M. M. Khader and M. Adel, “Numerical and theoretical treatment based on the compact finite difference and spectral collocation algorithms of the space fractionalorder Fisher’s equation,” International Journal of Modern Physics C, vol. 31, no. 9, pp. 1–13, 2020. View at: Publisher Site  Google Scholar
 S. Kumar, A. Kumar, B. Samet, and H. Dutta, “A study on fractional hostparasitoid population dynamical model to describe insect species,” Numerical Methods for Partial Differential Equations, vol. 37, no. 2, pp. 1673–1692, 2021. View at: Publisher Site  Google Scholar
 P. Veeresha, D. G. Prakasha, J. Singh, I. Khan, and D. Kumar, “Analytical approach for fractional extended Fisher–Kolmogorov equation with Mittag–Leffler kernel,” Advances in Difference Equations, vol. 174, no. 1, pp. 1–17, 2020. View at: Google Scholar
 J. E. MacíasDíaz and A. Puri, “An explicit positivitypreserving finitedifference scheme for the classical FisherKolmogorovPetrovskyPiscounov equation,” Applied Mathematics and Computation, vol. 218, no. 9, pp. 5829–5837, 2012. View at: Publisher Site  Google Scholar
 V. Chandraker, A. Awasthi, and S. Jayaraj, “Implicit numerical techniques for Fisher equation,” Journal of Information and Optimization Sciences, vol. 39, no. 1, pp. 1–13, 2018. View at: Publisher Site  Google Scholar
 M. Izadi, “Splitstep finite difference schemes for solving the nonlinear Fisher equation,” Journal of Mahani Mathematical Research Center, vol. 7, no. 1, pp. 37–55, 2018. View at: Google Scholar
 M. Izadi, “A secondorder accurate finitedifference scheme for the classical FisherKolmogorovPetrovskyPiscounov equation,” Journal of Information and Optimization Sciences, vol. 42, no. 2, pp. 431–448, 2021. View at: Publisher Site  Google Scholar
 X. Zhang, L. Yao, and J. Liu, “Numerical study of Fisher’s equation by the RBFFD method,” Applied Mathematics Letters, vol. 120, Article ID 107195, 2021. View at: Google Scholar
 R. K. Saeed and A. A. Mustafa, “Numerical solution of Fisher—KPP equation by using reduced differential transform method,” In AIP Conference Proceedings, vol. 1888, Article ID 020045, 2017. View at: Google Scholar
 M. ElHachem, S. W. McCue, and M. J. Simpson, “Invading and receding sharpfronted travelling waves,” Bulletin of Mathematical Biology, vol. 83, no. 4, pp. 1–25, 2021. View at: Publisher Site  Google Scholar
 J. Xu, F. Zhao, Z. Sheng, and G. Yuan, “A nonlinear finite volume scheme preserving maximum principle for diffusion equations,” Communications in Computational Physics, vol. 29, no. 3, pp. 747–766, 2021. View at: Publisher Site  Google Scholar
 T. Kadri and K. Omrani, “A fourthorder accurate finite difference scheme for the extendedFisher–Kolmogorov equation,” Bulletin of the Korean Mathematical Society, vol. 55, no. 1, pp. 297–310, 2020. View at: Google Scholar
 A. Başhan, Y. Ucar, N. M. Yamurlu, and A. Esen, “Numerical solutions for the fourth order extended Fisher–Kolmogorov equation with high accuracy by differential quadrature method,” Sigma Journal of Engineering and Natural Sciences, vol. 9, no. 3, pp. 273–284, 2018. View at: Google Scholar
 O. P. Yadav and R. Jiwari, “Finite element analysis and approximation of Burgers’Fisher equation,” Numerical Methods for Partial Differential Equations, vol. 33, no. 5, pp. 1652–1677, 2017. View at: Publisher Site  Google Scholar
 F. Bonizzoni, M. Braukhoff, A. Jüngel, and I. Perugia, “A structurepreserving discontinuous Galerkin scheme for the FisherKPP equation,” Numerische Mathematik, vol. 146, no. 1, pp. 119–157, 2020. View at: Publisher Site  Google Scholar
 M. Ilati and M. Dehghan, “Direct local boundary integral equation method for numerical solution of extended FisherKolmogorov equation,” Engineering with Computers, vol. 34, no. 1, pp. 203–213, 2018. View at: Publisher Site  Google Scholar
 Ö. Oruç, “An efficient wavelet collocation method for nonlinear twospace dimensional FisherKolmogorovPetrovskyPiscounov equation and twospace dimensional extended FisherKolmogorov equation,” Engineering with Computers, vol. 36, no. 3, pp. 839–856, 2020. View at: Publisher Site  Google Scholar
 P. Castillo and S. Gómez, “An interpolatory directional splittinglocal discontinuous Galerkin method with application to pattern formation in 2D/3D,” Applied Mathematics and Computation, vol. 397, Article ID 125984, 2021. View at: Publisher Site  Google Scholar
 M. H. Kabir, “Reactiondiffusion modeling of the spread of spruce budworm in boreal ecosystem,” Journal of Applied Mathematics and Computing, vol. 66, no. 1, pp. 203–219, 2021. View at: Publisher Site  Google Scholar
 S. W. McCue, M. ElHachem, and M. J. Simpson, “Exact sharpfronted travelling wave solutions of the FisherKPP equation,” Applied Mathematics Letters, vol. 114, Article ID 106918, 2021. View at: Publisher Site  Google Scholar
 J. Cai and H. Gu, “Asymptotic behavior of solutions of free boundary problems for Fisherkpp equation,” Journal of Dynamics and Differential Equations, vol. 33, no. 2, pp. 913–940, 2021. View at: Publisher Site  Google Scholar
 J. Brasseur, “The role of the range of dispersal in a nonlocal Fisher–KPP equation: an asymptotic analysis,” Communications in Contemporary Mathematics, vol. 23, no. 3, Article ID 2050032, 2021. View at: Google Scholar
 R. B. Salako and W. Shen, “Long time behavior of random and nonautonomous Fisherkpp equations: Part Istability of equilibria and spreading speeds,” Journal of Dynamics and Differential Equations, vol. 33, no. 2, pp. 1035–1070, 2021. View at: Publisher Site  Google Scholar
 J. F. Leyva and R. G. Plaza, “Spectral stability of traveling fronts for reaction diffusiondegenerate Fisherkpp equations,” Journal of Dynamics and Differential Equations, vol. 32, no. 3, pp. 1311–1342, 2020. View at: Publisher Site  Google Scholar
 G. Tian, Z.C. Wang, and G.B. Zhang, “Stability of traveling waves of the nonlocal FisherKPP equation,” Nonlinear Analysis, vol. 211, Article ID 112399, 2021. View at: Publisher Site  Google Scholar
 C. Lee, H. Kim, S. Yoon et al., “An unconditionally stable scheme for the AllenCahn equation with highorder polynomial free energy,” Communications in Nonlinear Science and Numerical Simulation, vol. 95, Article ID 105658, 2021. View at: Publisher Site  Google Scholar
 S. Kim, D. Jeong, C. Lee, and J. Kim, “Finite difference method for the multiasset Black–Scholes equations,” Mathematics, vol. 8, no. 3, p. 391, 2020. View at: Publisher Site  Google Scholar
 X. Y. Wang, “Exact and explicit solitary wave solutions for the generalised Fisher equation,” Physics Letters A, vol. 131, no. 45, pp. 277–279, 1988. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2021 Sangkwon Kim et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.