# An Unconditionally Stable Positivity-Preserving Scheme for the One-Dimensional Fisher–Kolmogorov–Petrovsky–Piskunov Equation

Authors: Sangkwon Kim; Chaeyoung Lee; Hyun Geun Lee; Hyundong Kim; Soobin Kwak; Youngjin Hwang; Seungyoon Kang; Seokjun Ham; Junseok Kim
Journal: Discrete Dynamics in Nature and Society
 (2021)

---

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


---

## Body

## 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:(1)∂u∂tx,t=D∂2u∂2xx,t+Kux,t1−upx,t,x∈Ω,t>0where ux,t is the population density at position x and time t. Here, D is the diffusion coefficient, K is the coefficient of nonlinearity, and p 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, ifp 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 Section2, 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 Ω=Lx,Rx. Let Nx be an integer and h=Rx−Lx/Nx be the uniform spatial step size. Let uin be the numerical approximation of uxi,tn, where xi=Lx+i−0.5h for i=1,2,…,Nx and tn=nΔt. 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):(2)∂u∂tx,t=D∂2u∂x2x,t.Equation (2) can be discretized using the FDM as(3)uin+1/2−uinΔt=Dui−1n+1/2−2uin+1/2+ui+1n+1/2h2,fori=1,…,Nx.Equation (3) is reorganized as(4)−Dh2ui−1n+1/2+2Dh2+1Δtuin+1/2−Dh2ui+1n+1/2=uinΔt.For the boundary condition, we use the zero Neumann boundary condition: anduNx+1n+1/2=uNxn+1/2 for all values of n+1/2. 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:(5)∂u∂tx,t=Kux,t1−upx,t.For simplicity of exposition, letK=1; then, equation (5) becomes(6)∂u∂tx,t=ux,t1−upx,t.To solve equation (6), we use a recently developed interpolation method [29]. Let us assume that the time step Δt is given. We uniformly discretize the interval 0,1, i.e., I=Is|Is=s−1/M−1,fors=1,…,M. Note that M is a positive integer. We solve equation (6) on I with a smaller time step Δτ=Δt/Nτ until t=Δt. Here, Nτ will be defined later. Figure 1 shows a simulation result with the initial condition Φs0=Is for s=1,…,M. That is, for m=0,…,Nτ−1,(7)Φsm+1Δτ=ΦsmΔτ+ΔτΦsmΔτ1−ΦspmΔτ.Figure 1
Numerical solutionΦΔt of equation (6) on the uniform grid I.Once we compute allΦsΔt for s=1,…,M, then we will use these values when we numerically solve equation (6) by using an interpolation.Next, we find a value ofNτ appropriate to stably integrate equation (7) for given Δt. If Is=0 or Is=1, then from equation (7), we have(8)Φsm+1Δτ=Is,form=0,…,Nτ−1.Next, when0<Ιs<1 and Φs0=Is, we want 0≤Φsm+1Δτ≤1 for m=0,…,Nτ−1. Here, to show 0≤ΦsΔτ≤1, by mathematical induction, it is sufficient to show 0≤Φsm+1Δτ≤1. Then, let us suppose that 0<Is<1, Φs0=Is, and m=0; from (7), we have(9)ΦsΔτ=Is+ΔτIs1−Isp.As shown in Figure2, Is1−Isp>0 for 0<Is<1; therefore, we have ΦsΔτ>0 from equation (9).Figure 2
Nonlinear partu1−up of equation (1) for p=1,3,7,20, and 100.Because0≤Is+ΔτIs1−Isp is always satisfied, we only need to find Δτ satisfying the following condition:(10)Δτ≤1−IsIs1−Isp.The right-hand side of inequality (10) is decreasing with respect to Is; see Figure 3.Figure 3
Semilogy plot of1−Is/Is1−Isp with respect to Is. Here, p=10 is used.Therefore, we consider the limit of the term asIs approaches 1. By using L’Hopital’s rule, we have(11)Δτ≤limIs⟶11−IsIs1−Isp=limIs⟶1−11−p+1Isp=1p.From the definition ofΔτ=Δt/Nτ and stability condition (11), Nτ≥pΔt should be satisfied. Therefore, unless otherwise specified, we set Nτ=pΔt+1, where x is the greatest integer less than or equal to x. Given the intermediate solution from the first step, uin+1/2, the second step is to solve equation (6) by using the precomputed values and interpolation. We can find an interval Is,Is+1 which satisfies Is≤uin+1/2≤Is+1 for some s and 1≤s≤M−1. Then, we define(12)uin+1=Is+1−uin+1/2Is+1−IsΦsΔt+uin+1/2−IsIs+1−IsΦs+1Δt,which is schematically illustrated in Figure 4.Figure 4
Schematic diagram of the interpolation.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 Ω=−30,120 with D=1 and K=1:(13)ux,0=0.5−0.5tanhpx22p+42/p.The exact solution is given in [31]:(14)ux,t=0.5−0.5tanhp22p+4x−p+4t2p+42/p.We run the simulation up to terminal timet=46 with different time steps: Δt=0.4, 0.2, 0.1, and 0.05. The other numerical parameter values used are M=800 and Nx=300. Figure 5 shows the numerical results and the exact solution (14) with p=1. We can observe that the numerical solutions converge to the exact solution as we refine the time step.Figure 5
Convergence of the numerical solutions with respect toΔt values at t=46.Table1 shows the l2-norm error and convergence rate of numerical solutions with the same conditions above except for the final time t=4 and the fixed Nx=1024, where Nt is the number of time iterations of the simulation for each Δt=t/Nt. We can observe that the numerical scheme is first-order accurate in time.Table 1
l2-norm error and temporal convergence rate of numerical solutions.
Nt=16RateNt=32RateNt=64RateNt=128Error8.8861e − 30.93284.6551e − 30.96362.3870e − 30.98121.2091e − 3Next, we investigate the accuracy of numerical results on the value ofM. Figure 6 shows the convergence of numerical solutions with respect to M values on subdomain 75,105. Here, we used the same parameter values as before except M values and Δt=0.01.Figure 6
Convergence of the numerical solutions with respect toM values.Last, we study the convergence of numerical results with respect toh. The applied conditions are the same as the time regarding simulation, except for Δt=0.0001. With all the other parameters fixed, the spatial step size Nx is refined to calculate the convergence rate. l2-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.Table 2
l2-norm error and temporal convergence rate of numerical solutions.
Nx=51RateNx=101RateNx=201RateNx=401Error2.1493e − 31.92655.6542e − 41.98811.4253e − 42.04233.4576e − 5Time 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 Ω=−30,120 with D=1, K=1, M=800, and Nx=600. We run the simulation up to terminal time t=30 with different time steps: Δt=0.03, 0.3, and 3. Figure 7 shows the numerical results and the exact solution (14) with p=1. 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.Figure 7
Numerical solutions with respect to three different temporal step sizesΔt=0.03,0.3,and3 at t=30.


### 3.3. Effect ofp

Now, let us consider the effect ofp on the evolution dynamics. Let us consider equation (6). By using the explicit Euler’s method, we present the evolution of u with the initial condition u0=0.1 on various values of p. Here, we use the time step Δt=0.0001 and the terminal time t=7:(15)un+1=un+Δtun1−unp.Figure8 shows the evolution of u with p=1,3,7,20, and 100. As the parameter p value increases, the solution u evolves faster.Figure 8
Evolution ofu with p=1,3,7,20, and 100, Δt=0.0001, and t=7.However, as shown in Figure8, there is little difference between p=20 and p=100. Next, we investigate the numerical solution of equation (9) on I as uM=400, Δt=0.0001, and t=6; see Figure 9.Figure 9
Numerical solutions ofΦ with p=1,3,7,20, and 100, M=400, Δt=0.0001, and t=6.To clearly identify this effect, numerical experiments are conducted using initial conditionux,0=0.3cosπx/5+0.3. We set the computational domain Ω=10, D=1.0, K=0.1, terminal time t=0.5, Δt=0.001, M=800, and Nx=400. Figure 10 shows the numerical solutions with various values of p=1,3,7,20, and 100.Figure 10
Snapshots of the numerical solutions of the Fisher–KPP equation with the initial conditionux,0=0.3cosπx/5+0.3. Here, t=0.5 with Δt=0.001 and p=1,3,7,20, and 100.The evolution with largep is faster than that with small p, especially near u=0.5. However, evolution with p in the vicinity of u=0 and u=1 has a small difference.Next, let us set an initial condition equation (13) on the computational domain Ω=−30,120 with D=1 and K=1. We run the simulation up to time t=10 with Δt=0.005, M=800, and Nx=300. Figure 11 shows the numerical results with different p values on subdomain −10,40.Figure 11
Numerical solutions with initial condition equation (13) for different p values. (a–c) Solutions with p=1,3, and 6, respectively.
(a)(b)(c)Figure12 shows the numerical solution with an initial condition ux,0=1 if x≤0; otherwise, ux,0=0 for different p values. Here, we used the same parameter values as before. We can observe that the numerical solutions become stiffer as the p value increases.Figure 12
Numerical solutions with differentp values.


### 3.4. Positive Preserving

We conduct numerical experiment using the initial condition(16)ux,0=0.5+0.4cosπx10+0.09randx.where −1<randx<1 to confirm that the algorithm satisfies positive preserving and boundedness by one. The parameters used in the numerical experiment are as follows: Ω=0,10, D=1, K=1, Δt=0.0005, M=800, and Nx=101. Figures 13(a) and 13(b) show the evolution of numerical solutions u for p=1 and p=20 at t=1.2. Figures 13(c) and 13(d) show the temporal evolution of the maximum and minimum of u over time when p=1 and p=20, respectively. We can confirm that 0<u<1 is satisfied regardless of the p value. The proposed algorithm is a positive-preserving scheme.Figure 13
(a, b) Numerical solution ofu with (a) p=1 and (b) p=20. (c, d) Temporal evolution of the minimum and maximum of u when (c) p=1 and (d) p=20.
(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:(17)−αRui−1n+1+ui+1n+1+1+2αR−βΔt1−uinp+1−βΔtuinuin+1=1−2R1−α+Δt1−βuin+R1−αui−1n+ui+1n,where α=β=1/2 and R=DΔt/h2 with D=1 and K=2. Let us consider an initial condition on Ω=−30,40: ux,0=1 if x≤0; otherwise, ux,0=0. We run the simulation up to time t=5 with Δt=0.05, M=800, and Nx=60. Figure 14 shows the numerical solutions with various values of p=1,80, and 800. When p=1, the numerical solutions for both methods converge to the exact solution. For the large value p, 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).Figure 14
Numerical solution comparison of the proposed method and recently developed method [10] for various values of p=1,80, and 800. (a–c) Numerical solutions with p=1,80, and 800, respectively.
(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 N. That is, the numerical solution is ux,t+Δt=ℒΔt°NΔtux,t+OΔt2 with time step size Δt. For a second-order accurate scheme, meanwhile, ux,t+Δt=ℒΔt/2°NΔt°ℒΔt/2ux,t+OΔt3. Note that all numerical schemes used for ℒ and N 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 Nτ 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.Figure 15
Convergence of the numerical solutions with respect toΔt values.Table 3
l2-norm error and temporal convergence rates of the second-order accurate scheme.
Nt=16RateNt=32RateNt=64RateNt=128Error2.7444e− 31.96967.0073e− 41.99361.7595e− 42.00064.3970e− 5We 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 Section3.1 is used. The derived results are given in Table 4, and it shows that the method is second-order accurate in space.Table 4
l2-norm error and temporal convergence rate of numerical solutions.
Nx=51RateNx=101RateNx=201RateNx=401Error2.1491e− 31.92685.6525e− 41.98921.4238e− 42.04793.4432e − 5


## 3.1. Convergence Tests

Let us consider an initial condition for equation (1) on the computational domain Ω=−30,120 with D=1 and K=1:(13)ux,0=0.5−0.5tanhpx22p+42/p.The exact solution is given in [31]:(14)ux,t=0.5−0.5tanhp22p+4x−p+4t2p+42/p.We run the simulation up to terminal timet=46 with different time steps: Δt=0.4, 0.2, 0.1, and 0.05. The other numerical parameter values used are M=800 and Nx=300. Figure 5 shows the numerical results and the exact solution (14) with p=1. We can observe that the numerical solutions converge to the exact solution as we refine the time step.Figure 5
Convergence of the numerical solutions with respect toΔt values at t=46.Table1 shows the l2-norm error and convergence rate of numerical solutions with the same conditions above except for the final time t=4 and the fixed Nx=1024, where Nt is the number of time iterations of the simulation for each Δt=t/Nt. We can observe that the numerical scheme is first-order accurate in time.Table 1
l2-norm error and temporal convergence rate of numerical solutions.
Nt=16RateNt=32RateNt=64RateNt=128Error8.8861e − 30.93284.6551e − 30.96362.3870e − 30.98121.2091e − 3Next, we investigate the accuracy of numerical results on the value ofM. Figure 6 shows the convergence of numerical solutions with respect to M values on subdomain 75,105. Here, we used the same parameter values as before except M values and Δt=0.01.Figure 6
Convergence of the numerical solutions with respect toM values.Last, we study the convergence of numerical results with respect toh. The applied conditions are the same as the time regarding simulation, except for Δt=0.0001. With all the other parameters fixed, the spatial step size Nx is refined to calculate the convergence rate. l2-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.Table 2
l2-norm error and temporal convergence rate of numerical solutions.
Nx=51RateNx=101RateNx=201RateNx=401Error2.1493e − 31.92655.6542e − 41.98811.4253e − 42.04233.4576e − 5Time 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 Ω=−30,120 with D=1, K=1, M=800, and Nx=600. We run the simulation up to terminal time t=30 with different time steps: Δt=0.03, 0.3, and 3. Figure 7 shows the numerical results and the exact solution (14) with p=1. 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.Figure 7
Numerical solutions with respect to three different temporal step sizesΔt=0.03,0.3,and3 at t=30.


## 3.3. Effect ofp

Now, let us consider the effect ofp on the evolution dynamics. Let us consider equation (6). By using the explicit Euler’s method, we present the evolution of u with the initial condition u0=0.1 on various values of p. Here, we use the time step Δt=0.0001 and the terminal time t=7:(15)un+1=un+Δtun1−unp.Figure8 shows the evolution of u with p=1,3,7,20, and 100. As the parameter p value increases, the solution u evolves faster.Figure 8
Evolution ofu with p=1,3,7,20, and 100, Δt=0.0001, and t=7.However, as shown in Figure8, there is little difference between p=20 and p=100. Next, we investigate the numerical solution of equation (9) on I as uM=400, Δt=0.0001, and t=6; see Figure 9.Figure 9
Numerical solutions ofΦ with p=1,3,7,20, and 100, M=400, Δt=0.0001, and t=6.To clearly identify this effect, numerical experiments are conducted using initial conditionux,0=0.3cosπx/5+0.3. We set the computational domain Ω=10, D=1.0, K=0.1, terminal time t=0.5, Δt=0.001, M=800, and Nx=400. Figure 10 shows the numerical solutions with various values of p=1,3,7,20, and 100.Figure 10
Snapshots of the numerical solutions of the Fisher–KPP equation with the initial conditionux,0=0.3cosπx/5+0.3. Here, t=0.5 with Δt=0.001 and p=1,3,7,20, and 100.The evolution with largep is faster than that with small p, especially near u=0.5. However, evolution with p in the vicinity of u=0 and u=1 has a small difference.Next, let us set an initial condition equation (13) on the computational domain Ω=−30,120 with D=1 and K=1. We run the simulation up to time t=10 with Δt=0.005, M=800, and Nx=300. Figure 11 shows the numerical results with different p values on subdomain −10,40.Figure 11
Numerical solutions with initial condition equation (13) for different p values. (a–c) Solutions with p=1,3, and 6, respectively.
(a)(b)(c)Figure12 shows the numerical solution with an initial condition ux,0=1 if x≤0; otherwise, ux,0=0 for different p values. Here, we used the same parameter values as before. We can observe that the numerical solutions become stiffer as the p value increases.Figure 12
Numerical solutions with differentp values.


## 3.4. Positive Preserving

We conduct numerical experiment using the initial condition(16)ux,0=0.5+0.4cosπx10+0.09randx.where −1<randx<1 to confirm that the algorithm satisfies positive preserving and boundedness by one. The parameters used in the numerical experiment are as follows: Ω=0,10, D=1, K=1, Δt=0.0005, M=800, and Nx=101. Figures 13(a) and 13(b) show the evolution of numerical solutions u for p=1 and p=20 at t=1.2. Figures 13(c) and 13(d) show the temporal evolution of the maximum and minimum of u over time when p=1 and p=20, respectively. We can confirm that 0<u<1 is satisfied regardless of the p value. The proposed algorithm is a positive-preserving scheme.Figure 13
(a, b) Numerical solution ofu with (a) p=1 and (b) p=20. (c, d) Temporal evolution of the minimum and maximum of u when (c) p=1 and (d) p=20.
(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:(17)−αRui−1n+1+ui+1n+1+1+2αR−βΔt1−uinp+1−βΔtuinuin+1=1−2R1−α+Δt1−βuin+R1−αui−1n+ui+1n,where α=β=1/2 and R=DΔt/h2 with D=1 and K=2. Let us consider an initial condition on Ω=−30,40: ux,0=1 if x≤0; otherwise, ux,0=0. We run the simulation up to time t=5 with Δt=0.05, M=800, and Nx=60. Figure 14 shows the numerical solutions with various values of p=1,80, and 800. When p=1, the numerical solutions for both methods converge to the exact solution. For the large value p, 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).Figure 14
Numerical solution comparison of the proposed method and recently developed method [10] for various values of p=1,80, and 800. (a–c) Numerical solutions with p=1,80, and 800, respectively.
(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 N. That is, the numerical solution is ux,t+Δt=ℒΔt°NΔtux,t+OΔt2 with time step size Δt. For a second-order accurate scheme, meanwhile, ux,t+Δt=ℒΔt/2°NΔt°ℒΔt/2ux,t+OΔt3. Note that all numerical schemes used for ℒ and N 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 Nτ 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.Figure 15
Convergence of the numerical solutions with respect toΔt values.Table 3
l2-norm error and temporal convergence rates of the second-order accurate scheme.
Nt=16RateNt=32RateNt=64RateNt=128Error2.7444e− 31.96967.0073e− 41.99361.7595e− 42.00064.3970e− 5We 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 Section3.1 is used. The derived results are given in Table 4, and it shows that the method is second-order accurate in space.Table 4
l2-norm error and temporal convergence rate of numerical solutions.
Nx=51RateNx=101RateNx=201RateNx=401Error2.1491e− 31.92685.6525e− 41.98921.4238e− 42.04793.4432e − 5


## 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].



---

*Source: 7300471-2021-10-19.xml*