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 formhas 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 . 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).

Proposition 1

The Dirichlet problem (1)–(4) has at most one solution.

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

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 is monotonically decreasing in t, meaning that in each case we expect to have

So, it follows that u(x, t) = 0, u = u for all t ≥ 0.□

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 = (u1, u2,…, u), D = diag(d1, d2, …, d), with d ≥ 0 for all i. 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 Rn → R, let ∥·∥∞ denotes L-norm, and ∥·∥.

Definition 2

B is acceptable if the following conditions hold:

B is a bounded subset and continuous functions on R, and if v ∈ R, then ∥v∥ ≥ ∥v∥∞.

B is the translation-invariant. That is, voϱ ∈ B for every v ∈ B and every translate ϱ:R → R. Also ∥voϱ∥ = ∥v∥.

If R → R is smooth, and f (0) = 0, then f (v) ∈ B for every v ∈ B, and for any G > 0, there exists a constant k (G) such that

for all v, v′ in B with ∥v∥∞, ∥v′∥∞ ≤ G.

If ϱh:R → R represents translation by h, that is ϱh (x) = x + h, then for each v ∈ B, 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) becomeswhere 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 modelSubject to initial and homogeneous boundary conditionswhere D > 0 and τ > 0. The solution is given aswith

By examining the solution reveals the condition that supports population growth and extinction.

For instance, ifthen u(x, t) will approach zero as time progresses, while ifu(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 towith 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 and a maximum height of . 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 u0(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, ∂Ω is smooth, u ∈ R, and f (u) = τ (u) is smooth in U ⊂ R → R for each t ≥ 0. D is nonnegative diagonal matrix of size n × n. So, we assume that system (7) has a bounded invariant regionwhere −∞ < a < b < ∞, i = 1, 2,…, n. In addition to (7), we have the initial conditionswhere u0(x) lies in Σ∀x ∈ Ω. We also assume that u is admissible of the homogeneous Neumann boundary conditions

Lemma 3

Let Σ be defined by (16) and ρ = ρ(u, t) is Lipschitz continuous in Σ, for each t ≥ 0. If

Then ρ+is Lipschitz continuous in Σ.

Let us assume that a = 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 so that

If ξ0 = max(0, + z − u), that is, () = max(0,  + z − u), j = 1, 2,…, n. It is obvious that  ∈ L, and that

Known thatand if  + z − u ≥ 0,  − max(0,  + z − u) = u − z, while  + z − u < 0, then  − max(0,  + z − u) =  < u − z. Thus (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 equationwith 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 equationsand u = [u1, u2,…, u]T, for 1 ≤ i, j ≤ l.

Again, the two-dimensional form of system (6) can be written asnow, we discretize the spatial domain by mesh (x, y) = (l1 + i × h, l1 + j × h) where h = (l2 − l1)/(N + 1), h = (l2 − l1)/(N + 1) and 0 ≤ i ≤ N + 1 and 0 ≤ j ≤ N + 1. Using fourth order central difference discretization on the linear term, we obtain a system of nonlinear ODEs of the formwhere

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 equationsso 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 bywhere ψ, j = 1, 2,… is the complete orthonormal system in (H), it follows that

Following this assumption, we require powers (−L) for 0 ≤ δ ≤ 1. In thatKnown, (−L) as (L) is regarded as closed operator, and (H) is a Banach space w.r.t graph norm ∥u∥ + ∥u∥. We have

Due to Parseval’s equation, we have from the first inequality due to 0 > λ1 > λ2 >···, limλ = −∞ that (−λ)2/(−λ1)2 ≥ 1, which means that the graph norm ≡ to the norm ∥u∥ which can be used as a norm on H.

Lemma 4

Let 0 ≤ δ ≤ 1, it follows that

Proof

For all x ∈ R, we choose ex ≤ e, x = −λt/δ, we obtain

From the definition,and by using the Parseval equation, we havewhere s1 = (δ/e) t− and s2 = ∥u∥2. Taken together, we have

Time-stepping method

By adopting the description in (Cox and Matthews 2002; Du and Zhu 2005; Kassam and Trefethen 2005; Krogstad 2005; Owolabi and Patidar 2016), Cox and Matthews fourth order exponential time differencing Runge–Kutta formula (Cox and Matthews 2002) was used to advance the resulting ODEs are:wherewhere for instance, if (27) is used in conjunction with ODE (25), L is the linear diffusion term, −Dk2 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 ODEabout a fixed point u0, such that gu0 + F (u0) = 0. As a result of linearizing, we obtainwhere u is the perturbation to u0 and λ = F′(u, t) at u(t) = u0 is either a diagonal or block matrix that contains the eigenvalues of F. For the fixed point u0 to be stable, it is required that Re(g + λ) < 0. On applying the ETDRK4 method to (29), a recurrence relation involving u and u + 1 is obtained (de la Hoz and Vadilo 2008). By introducing the notation x = λh, y = gh, the amplification factor is given aswhereobviously, as y → 0, approximation in (31) reduces towhich 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, 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.

Fig. 1

The boundary of stability regions of the ETDRK4 method in the complex plane x for some negative values of y

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 byrespectively, where e and c 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 conditionswhere 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, the successive profile in (a) is obtained at T = 0.01(0.005)0.05. Panel (b) and the contour plots (d) are obtained for [0, 20], T = 0.05. Surface plot (c) is obtained at parameter values T = 0.05, in the interval [−2, 2]. Panels (e) is obtained at parameter value T = 0.1 in the interval [0, 1] while (f) obtained at final time T = 0.05 on domain [−1, 1]. The results presented have shown various behavioral patterns that could evolved when the patch size of one-dimensional KISS model (34) is varied in spatial domain.

Fig. 2

Numerical solutions of one-dimensional KiSS model (34). Time dependent density profiles u(x, t) versus position x on interval −l ≤ x ≤ l for D = τ = α = 1

The plots in Fig. 3 indicate the results from initial time (t0) to final time T showing the density profiles u(x, t) versus position x on a closed interval −l ≤ x ≤ l for the choice of growth rate τ = 0.5 and critical exponent α = 2. The successive profile in (a) is obtained at T = 0.01 with D = 0.5 on [−1, 1]. For panel (b), large D = 2, T = 0.05 on [0, 3]. For (c), D = 0.1, T = 0.01 on [0, 1]; plots (d) is obtained on [−4, 4], T = 0.02 for large diffusion coefficient D = 1.5. Panel (e) is obtained on the spatial domain [0, 5] with T = 0.1 and D = 0.05. Contour plot (f) is obtained with parameter values D = 0.2, T = 0.1 on domain of size [0, 4]. The results presented here have equally revealed some of the dispersal-driven patterns that arise as a result of diffusion.

Fig. 3

Diffusion-driven spatial patterning in one-dimensional KISS model (34) as it changes with both spatial domain x with varying time t

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.

Fig. 4

Performance of ETDRK4, ETDM4, ETDM5, ETDM6 and ETDADAMS4 methods for solving the KISS Eq. (34) with parameter values T = 1, D = 0.5, τ = 0.5 and α = 2 on interval [−1, 1] for N = 200

Table 1

Error norm ||L 2|| at some selected time steps for solving KiSS Eq. (34)

Method	Time step (k)	||L 2||	CPU time (s)
ETDM4	1/32	0.2150	1.5086
	1/64	0.0301	1.5467
	1/256	3.3910e−004	1.7189
	1/1024	3.4784e−006	3.8817
ETDM5	1/32	0.2170	1.9896
	1/64	0.0195	2.0277
	1/256	1.7793e−005	2.3538
	1/1024	1.7923e−006	4.6010
ETDM6	1/32	0.0664	2.0405
	1/64	8.5299e−004	2.0515
	1/256	4.2840e−004	2.2500
	1/1024	1.8248e−007	8.2387
ETDADAMS4	1/32	0.0050	1.4850
	1/64	0.0028	1.5511
	1/256	1.3334e−004	1.8339
	1/1024	2.3963e−006	4.4020
ETDRK4	1/32	7.2013e−004	1.1146
	1/64	5.0921e−005	1.0236
	1/256	3.5038e−007	1.3198
	1/1024	4.4086e−009	3.8726

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/∂x2 + ∂2/∂y2. One-dimensional 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; Owolabi and Patidar 2015, 2016), we consider the two-dimensional casewhere 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 = h = 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.

Fig. 5

Solutions of two-dimensional KISS model (35) on different spatial domain. The initial and parameter values: u0(x, y) = sinc[√(x/π)2 +√ (y/π)2], with D = 1, τ = 0.01, T = 0.5 and α = 2 on a l = 10 and b l = 20. Plots c T = 1.2, l = 10 and d T = 2, l = 20 are obtain with initial data u0(x, y) = cos(x) cos(y) exp(−√(x 2 + y /4))

Fig. 6

Solutions of two-dimensional KISS model (35) on different spatial domain. The initial and parameter values: u0(x, y) = exp(1/10)(cos x + sin y), with T = 0.05 D = 1, τ = 0.01 and α = 2 on a l = 8 and b l = 14

Fig. 7

Solutions of two-dimensional KISS model (35) on different spatial domain. The initial and parameter values are: u 0 = sin(√(x/π)2 +√(y/π)2), with D = 1, τ = 0.01, α = 2 and T = 1 for l = 10

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, ETDM5, ETDM6) methods all are of the higher orders. This approach was tested with the reaction–diffusion equation, a nonlinear form of KiSS model that was named after Kierstead and Slobodkin (1953) and Skellam (1951), which was originally developed to investigate the size of nutrient patches needed to sustain phytoplankton blooms. We carried out numerical simulations in both one- and two-dimensional space on spatial domain x ∈ [−l, l], that are chosen large enough to support the boom. Our numerical results revealed that the population size increases if the domain size l also increased. Some initial data and parameter values are chosen to mimic some existing patterns. The methodology presented in this paper can be extended to higher-order fractional derivative (Atangana and Nieto 2015; Atangana 2015) and time-dependent parabolic problems.