Effect of spatial configuration of an extended nonlinear Kierstead–Slobodkin reaction-transport model with adaptive numerical scheme

In this paper, we consider the numerical simulations of an extended nonlinear form of Kierstead–Slobodkin reaction-transport system in one and two dimensions. We employ the popular fourth-order exponential time differencing Runge–Kutta (ETDRK4) schemes proposed by Cox and Matthew (J Comput Phys 176:430–455, 2002), that was modified by Kassam and Trefethen (SIAM J Sci Comput 26:1214–1233, 2005), for the time integration of spatially discretized partial differential equations. We demonstrate the supremacy of ETDRK4 over the existing exponential time differencing integrators that are of standard approaches and provide timings and error comparison. Numerical results obtained in this paper have granted further insight to the question ‘What is the minimal size of the spatial domain so that the population persists?’ posed by Kierstead and Slobodkin (J Mar Res 12:141–147, 1953), with a conclusive remark that the population size increases with the size of the domain. In attempt to examine the biological wave phenomena of the solutions, we present the numerical results in both one- and two-dimensional space, which have interesting ecological implications. Initial data and parameter values were chosen to mimic some existing patterns.


Background
The study of nonlinear reaction-diffusion equation of the form (1) ∂u ∂t = D∇ 2 u + f (u) for 0 < x < l, t > 0 (2) u(x, 0) = d (x) (3) u(0, t) = a(t) (4) u(l, t) = b(t) has a long-standing history in mathematical modeling of propagation phenomena that mostly occurs in distributed dissipative dynamics, as well as diffusive transport of mass and heat, there is some reaction term representing, for instance, population growth or heat generation. Most realistic physical problems such as Allen-Chan, Burgers, Cahn-Hilliard, Fisher-KPP, Nagumo, Gray-Scott (or cubic autocatalytic), Kierstead, Slobodkin and Skellam (KiSS), Kuramoto-Sivashinsky and host of others, naturally exist in form of higher-order partial differential equations. In Holmes et al. (1994), PDEs of the class (1) have shown to provide a natural framework for investigating the influence patch size and geometry on the population dynamics of organisms living within an habitat. Many researchers have used equations of the form (1) in different forms especially in relation to three applications that model the behavior of biological systems in a spatial setting. The three major and popular applications of reaction-diffusion models relate to critical patch size (Kierstead and Slobodkin 1953), spread of advantageous genes (Fisher 1937), and pattern formation (Turing 1952). For instance, if the reaction or interaction term f(u) is replaced by κu(1 − u), where κ and D are positive parameters regarded as the carrying capacity and diffusion respectively in the context of biology, then (1) becomes the classic simplest case of a nonlinear reaction-diffusion equation popularly referred to as Fisher equation (1937) with history dated back to 1937, which has since becomes one of the most well-studied reaction-diffusion models in population biology to describe the spread of an advantageous allele.
In addition, A large number of nonlinear reaction-diffusion PDEs of the form (1), sometimes exist in the form of blow-up problems if f (u) = u α , α > 0. In PDEs systems, solution with the given initial data often lead to singularity (which could either be a point where a discontinuity occurs or the dependent variables tend to infinity) in finite time, such a phenomenon is widely referred to as the blow-up, and the time at which such occurs is known as the blow-up time, see for instance, (Owolabi 2015b;Roberts 2007) and references therein. A generalization of (1), which frequently occurs as a limiting case of a system, is when there is nonlocal dependency, often of the source term, upon u, we do not treat this problem here. We briefly present a review in "An extended KiSS model" section for special nonlinear KiSS model.
It is important to check for the uniqueness of solution of the class of diffusion Eq. (1).
Proof Let us adopt energy integral method for this proof, by assuming that there exist two solutions u α (x, t) and u β (x, t); let u = u α − u β . Based on the definition, u satisfies Furthermore, we obtain u t − Du xx = 0, for 0 < x < l and t > 0, u(x, 0) = 0, u(0, t) = 0 and u(l, t) = 0.
On integration, we get The second term to the right is zero from the boundary condition, thus we can take the t derivative as This implies that the integral l 0 is monotonically decreasing in t, meaning that in each case we expect to have We can also obtain the qualitative information based on the time of local existence of solution. By considering the system (1) in a single space variable u = (u 1 , u 2 ,…, u n ), We also in addition consider the initial data Next, solutions that are continuous function of time can be obtained in Banach spaces. We proceed by using definition to describe these space. Let B denote a Banach space of functions in R n → R, let ∥·∥ ∞ denotes L ∞ -norm, and ∥·∥ B .
Definition 2 B is acceptable if the following conditions hold: (i) B is a bounded subset and continuous functions on R, and if v ∈ R, then ∥v∥ B ≥ ∥v∥ ∞ . (ii) B is the translation-invariant. That is, voϱ ∈ B for every v ∈ B and every translate ϱ: we have The aims of this work are in folds. We first give an extension to the existing linear KiSS model in the form of a nonlinear type. Secondly, we formulate a viable numerical techniques in space and time for the numerical simulation of an extended KiSS model in one and two dimensional spaces. We finally justify the suitability and applicability of the present numerical techniques with the existing family of higher-order time-stepping schemes.
The rest of paper is structured as follows. In "An extended Kiss model" section, an extension is given to the KiSS model with underlying theory to back the patchy size selection. We discuss adaptive numerical methods in space and time as well as the stability analysis of the scheme in "Numerical method" section. Numerical experiments in one-and two dimensions are examined in "Numerical simulations" section. Finally, conclusion is drawn in "Conclusion" section.

An extended Kiss model
In the present paper, numerical solution of an exponential growth model of the form (1), where the reaction term f (u) is given as τ uα, so that Eq. (1) becomes where the diffusion coefficient D, the growth rate τ and the critical exponent constant α, are all positive parameters. This equation is the critical patch model popularly known as the KiSS model named after Skellam (1951) and Kierstead and Slobodkin (1953) which was originally developed to describe the spread of red tide outbreaks. Red tide is a name given to the discolored waters caused by the aggregation or blooming of microscopic organisms. A model for growth and spread of a population is used to determine the minimal size of the spatial domain needed for population to survive and this minimal size is is referred to as the critical patch size (Allen 2007).
In the classical paper (Kierstead and Slobodkin 1953), the critical patch size was determined for a simple reaction-diffusion equation with exponential growth, their model was applied to study phytoplankton plants living in the ocean. Determination of patch size of one-dimensional form of (6) have been considered (Allen 2007;Kierstead and Slobodkin 1953;Kot 2001;Murray 1989;Okubo 1978) on the spatial domain [0, l] via separation of variables method. In one-dimension with the choice α = 1, we have the KiSS model Subject to initial and homogeneous boundary conditions where D > 0 and τ > 0. The solution is given as By examining the solution reveals the condition that supports population growth and extinction.
For instance, if then u(x, t) will approach zero as time progresses, while if u(x, t) will increase indefinitely with time, thus leading to the bloom of the plankton. An attempt to have a better understanding of how the solution of the reaction-diffusion Eq. (6) behaves, we let τ = 0. Hence, Eq. (6) reduces to diffusion equation. We can now find the solution of the general initial value problem of solving (6) in spatial variable x, subject to with the aid of the Fourier transforms, as then u(x, t) will approach zero as time progresses, while if On using the initial condition as the localized source of the spread of species population, u0(x) = δ(x), then, (13) becomes As time increases the solution spreads out, having a typical width of O √ 4πDt and a maximum height of 1/ √ 4πDt. It is also noticeable that the diffusion transports the species within the interval of integration [−l, l], since u(x, t) > 0 for all x when t > 0. For |x| ≫ 1 and t ≪ 1, the corresponding species concentration are very small. If u 0 (x) = G(−x), then the solution takes the form Further, we consider the reaction-diffusion system (7) on (x, t) ∈ Ω × R+ , for Ω is defined as a bounded domain in R m , ∂Ω is smooth, u ∈ R n , and f (u) = τ (u) is smooth in U ⊂ R n → R n for each t ≥ 0. D is nonnegative diagonal matrix of size n × n. So, we assume that system (7) has a bounded invariant region In addition to (7), we have the initial conditions where u 0 (x) lies in Σ∀x ∈ Ω. We also assume that u is admissible of the homogeneous Neumann boundary conditions Proof Let us assume that a i = 0, i = 1, 2,…, n. For any u ∈ Σ, the set then ρ + (u, t) = sup{ρ(ξ, t): ξ ∈ Lu}. If we let u, z ∈ Σ, and since ρ is a continuous, there exists ξ ′ ∈ L u so that

and that
Known that (19) is admissible, and by suppressing the t's we obtain If the roles of u and z are interchanged, the proof is completed. □ In spite of considerable progress so far made in the field of population dynamics some years back, there are still many open problems. In particular, the numerical exploration of (6) for α > 1 has received little or no attention when the domain of interaction is considered wide enough to contain the population spread. Put together all these findings, we are motivated to seek for an appropriate and efficient numerical solution of (6) in one and two-dimensional space which we consider on an infinite domain truncated at some large, but finite value of l. We proceed in the next section to describe these methods.

Numerical method
We discuss briefly the spatial discretization methods used in in this paper. When a time dependent partial differential equation is discretized in space especially with either a finite difference or spectral approximations, it results to system of coupled ordinary differential equations in time, the resulting ODEs coming from the notion of method of lines (MOL) (Owolabi and Patidar 2014a) is stiff, such a system requires the use of higher-order approximation scheme in both space and time since naturally some of these time-dependent problems are found of combining lower-order nonlinear terms with higher-order linear terms. In one-dimension, we consider the semi-linear partial differential equation with D > 0, τ > 0 and α > 0.

Spatial discretization method
We discretize in space with step size h = x/(N − 1) and approximate the second-order spatial derivative by the fourth order central difference operator, we obtain a system of nonlinear ordinary differential equations and u = [u 1 , u 2 ,…, u l ]T, for 1 ≤ i, j ≤ l. Again, the two-dimensional form of system (6) can be written as now, we discretize the spatial domain by mesh ( where h x = (l 2 − l 1 )/(N x + 1), h y = (l 2 − l 1 )/(N y + 1) and 0 ≤ i ≤ N x + 1 and 0 ≤ j ≤ N y + 1. Using fourth order central difference discretization on the linear term, we obtain a system of nonlinear ODEs of the form where Spatial discretisation of Eq. (6) can also be done using Fourier spectral method with periodic boundary conditions (Boyd 2001;Craster and Sassi 2006;de la Hoz and Vadilo 2008;Kassam and Trefethen 2005;Trefethen 2000;Weideman and Reddy 2001). We adapt the Fourier spectral method from (Trefethen 2000) and applied it to (6), leaving all the time stepping in Fourier space gives the following system of ordinary differential equations so that the linear term of (6) now becomes a diagonal. Next, the systems (21), (23) and (25) will now be integrated using a time integration method as explained below.
At this junction, we have discretised the system (6) in spatial variables with both finite difference approximations and spectral approximations, we can now present the system of ODEs obtained in the form Let us use the idea of the so-called abstract ODEs in Hilbert space (H) to nonlinear problem (26). Let us take an assumption that the linear operator L:D(L) ⊂ H → H is defined by where ψ j , j = 1, 2,… is the complete orthonormal system in (H), it follows that Following this assumption, we require powers (−L) δ for 0 ≤ δ ≤ 1. In that Known, (−L) δ as (L) is regarded as closed operator, and (H δ ) is a Banach space w.r.t graph norm ∥u∥ + ∥u∥ δ . We have (24) u =      u 1,1 u 1,2 · · · u 1,N y u 1,N y +1 u 2,1 u 2,2 · · · u 2,N y u 2,N y +1 . . . . . . . . . . . .

Lemma 4 Let
where where for instance, if (27) is used in conjunction with ODE (25), L is the linear diffusion term, −Dk 2 and F represents the term τ uα which could be linear or nonlinear.

Stability analysis
We follow the general stability idea as discussed in (Beylkin et al. 1998;Cox and Matthews 2002;Fornberg and Driscoll 1999;Owolabi and Patidar 2014b) for the numerical scheme that incorporates the use of different methods for the treatment of both the linear and nonlinear parts of Eq. (26). To examine the stability of ETDRK4 method (27), we linearize the nonlinear autonomous ODE about a fixed point u 0 , such that gu 0 + F (u 0 ) = 0. As a result of linearizing, we obtain where u is the perturbation to u 0 and λ = F′(u, t) at u(t) = u 0 is either a diagonal or block matrix that contains the eigenvalues of F. For the fixed point u 0 to be stable, it is required that Re(g + λ) < 0. On applying the ETDRK4 method to (29), a recurrence relation involving u n and u n + 1 is obtained (de la Hoz and Vadilo 2008). By introducing the notation x = λh, y = gh, the amplification factor is given as where obviously, as y → 0, approximation in (31) reduces to (27) u n+1 = u n e Lh + F (u n , t n ) −4 − Lh + e L�t 4 − 3Lh + L 2 h 2 + 2(F (a n , t n + h/2) + F (b n , t n + h/2)) 2 + Lh + e Lh (−2 + Lh) a n = u n e Lh/2 + e Lh/2 − I F (u n , t n )/L, b n = u n e Lh/2 + e Lh/2 − I F (a n , t n + h/2)/L, c n = u n e Lh/2 + e Lh/2 − I (2F (b n , t n + h/2) − F (u n , t n ))/L, which corresponds to the amplification factor of fourth-order Runge-Kutta method. Hence, we continue with our analysis by taking the real negative values of y in the complex x plane where |r| < 1, by setting r = e iθ , with θ ∈ [0, 2π]. The curves in Fig. 1 correspond to y = 0, −3.5, −5, −7, −9, −11 from the inner curve to the outer curve. It is noticeable that the stability region of the inner curve obtained at y = 0, coincides with the stability region of classical fourth-order Runge-Kutta method.

Numerical simulations
To examine the efficiency and accuracy of our approach for ETDRK4 methods, we consider the numerical simulations of system (6) in one and two dimensions. We further justify the supremacy of ETDRK4 in comparison with the existing standard schemes of higher-orders by reporting the root mean square norm error of the solution defined by respectively, where e j and c j are the exact and computed values of the solution u at point j, and n is the number interior points.

Test 1: one-dimensional nonlinear KiSS model
In one-dimension, we consider the KiSS model of Kierstead and Slobodkin (1953) and Skellam (1951), subject to initial and zero-flux boundary conditions where u(x, t) is the density of the organisms at spatial domain x and time t, τ and α are both positive parameters, and D is the diffusion coefficient that measures the rate of dispersal. The particular choice of boundary conditions indicates that the organisms cannot boom or live beyond the domain. This assumption is taken to ensure that the experiment is not influenced by any external factor. All simulations here run for N = 200.
In Fig. 2 It is clear from the result presented in Fig. 4, that ETDRk4 has the best convergence when compared to other exponential time differencing schemes, such as ETDM4, ETDM5, ETDM6 and ETDADAMS4 methods (for details of these schemes, see Hochbruck and Ostermann 2011;Hochbruck et al. 2010). In Table 1, we illustrate the tradeoff between the computational [CPU (s)] time and the accuracy as time step k is refined for each of the methods with parameter values T = 1, D = 0.5, τ = 0.5 and α = 2 on interval [−1, 1] for N = 200. We display accuracy as a function of CPU time respectively for each of the schemes, to add a competing factor in differentiating between the methods.

Test 2: two-dimensional nonlinear KiSS model
Our major aim in this paper is to examine the behavior of system (6) numerically in two-dimensional space, that is, when the Laplacian operator ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 . Onedimensional form of KiSS equations are relatively simple to undertake using method of lines coupled with spatial adaptive schemes. In-fact, solutions of the form (6) have been sought theoretically (Allen 2007;Kot 2001;Okubo 1978). Unfortunately, in two space dimensions, numerical solutions of KiSS model (22) still requires some attentions, since simulations based upon the more conventional ideas become more time consuming. In the spirit of (Owolabi 2015a;Patidar 2015, 2016), we consider the twodimensional case (35) ∂u ∂t = D ∂ 2 u ∂x 2 + ∂ 2 u ∂y 2 + τ u α , (x, y) ∈ Ω = (l 1 ≤ x, y ≤ l 2 ), t > 0, u(x, y, 0) = u 0 (x, y), l 1 ≤ x, y ≤ l 2 , u(0, t) = u(l 2 , t) = 0, t > 0, where u(x, y, t) is the density of organisms at spatial coordinates x, y and time t. D > 0 remains the diffusion coefficient, while τ > 0 and α ≥ 1 are the respective growth rate and critical exponent. The initial data and parameter values were carefully chosen to make the Figures replicate some of the existing patterns. In all cases, the space step h was kept equal to l, that is, h x = h y = l in the spatial domain −l ≤ x, y ≤ l. Figures 5, 6, 7 show various patterns formation process that could emerge as a result of variation of the initial data.
The two-dimensional results presented in Figs. 5, 6, 7 are much more meaningful in the contexts of mathematical biology and ecology. It is clear from the numerical results that the increase in the spatial domain l is the factor responsible for the spread of the phytoplankton plants living in the ocean.

Conclusion
In this paper, we have further justified the assertion made by Kassam and Trefethen (2005) on the efficiency and suitability of ETDRK4 schemes in conjunction with spatial discretization methods by comparing it with exponential time differencing method (ETDADAMS4) of Adams type and exponential time differencing multistep (ETDM4,