Abstract
In this study, we present an unconditionally stable positivity-preserving numerical method for the Fisher–Kolmogorov–Petrovsky–Piskunov (Fisher–KPP) equation in the one-dimensional space. The Fisher–KPP equation is a reaction-diffusion 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 positivity-preserving 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 collocation-type method [3, 4], Adams–Bashforth–Moulton method [5], and Mittag-Leffler 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 positivity-preserving 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. El-Hachem 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–KPP-type model that simulates biological recession. Xu et al. [14] proposed a nonlinear finite volume scheme for two-dimensional diffusion equations which preserves the maximum principle. Fourth-order 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 positivity-preserving numerical method for the Fisher–KPP equation has been proposed [7]. Furthermore, the positivity-preserving schemes for the Fisher–KPP equation were developed using continuous finite element discretizations [17], and discontinuous Galerkin schemes were proved to be entropy-stable [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 sharp-fronted travelling wave solutions of the Fisher–KPP equation to the Fisher–Stefan-type 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 positivity-preserving 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 one-dimensional 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 right-hand 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 first-order 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 second-order 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 positive-preserving 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 second-order 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. Second-Order Accurate Scheme
In this section, we present a second-order accurate scheme in time. Previously, the proposed scheme, which is a first-order 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 second-order accurate scheme, meanwhile, . Note that all numerical schemes used for and should be at least second-order 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 second-order 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 second-order accurate in time.
We tested the spatial convergence for the aforementioned scheme. It can be assumed that it is second-order 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 second-order accurate in space.
4. Conclusions
In this study, we presented an unconditionally stable positivity-preserving 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 positivity-preserving 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.