Analysis of an efficient integrator for a size-structured population model with a dynamical resource

Computers and Mathematics with Applications xx (xxxx) xxx–xxx Contents lists available at ScienceDirect Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa Analysis of an efficient integrator for a size-structured population model with a dynamical resource Q1 ∧ O. ∧ Anguloa,∗,∧ J.C. ∧ López-Marcosb,∧ M.A. ∧ López-Marcosb aDepartamento de Matemática Aplicada, ETSIT, Universidad de Valladolid, Pso. Belén 5, 47011 Valladolid, Spain bDepartamento de Matemática Aplicada, Facultad de Ciencias, Universidad de Valladolid, Valladolid, Spain article info Article history: Available online xxxx Keywords: Structured population models Numerical methods Convergence Daphnia magna abstract In this paper, an efficient numerical method for the approximation of a nonlinear sizestructured population model is presented. The nonlinearity of the model is given by dependency on the environment through the consumption of a dynamical resource. We ∧ analyse the properties of the numerical scheme and optimal second-order convergence is derived. We report experiments with academical tests to demonstrate numerically the predicted accuracy of the scheme. The model is applied to solve a biological problem: the dynamics of an ectothermic population (the water flea, Daphnia magna). We ∧ analyse its long time evolution and describe the asymptotically stable steady states, both equilibria and limit cycles. ©2014 Elsevier Ltd. All rights reserved. 1. Introduction1 In nature we observe that some physiological characteristics, like age, level of satiation (of a predator), energy reserves2 or the body size of the individual, play an important role in its behaviour. Physiologically structured population models ex-3 press the dynamics of the population in terms of the processes taking place at the basic unit (individual level) considering4 physiological differences. They describe the changes in the number of individuals of a population due to growth, death and5 reproduction and reflect the effect of the physiological state of the individuals on the population dynamics.6 Unlimited population growth does not exist either. A population influences its environment and therefore its own be-7 haviour. In addition, many biological feedback loops can only be described properly in terms of the interaction of the8 physiological processes within the individuals (e.g. the availability of food). Consequently, the use of nonlinear structured ∧ 9 population models is an ideal tool to give a realistic mathematical formulation of density dependence. They allow us to take10 into account the effect of competition for natural resources in the structured-specific growth, mortality and fertility rates.11 In this paper, we study a situation in which the vital rates are influenced by the availability of food in the environment, for12 which the individuals in the population compete. The dynamics of food density is determined by two processes: the regen-13 eration of the food within the environment (which models the changes in the food density in the absence of any consumers)14 and the feeding by the individuals of the physiological structured population. It could be explained as a predator–prey model15 in which the predator is considered physiologically structured and the prey is unstructured. A similar theoretical framework16 could be found when we try to model the dynamics of a cellular population for a continuously stirred batch culture in a tank17 reactor [1]. In this case, the model includes a different integral term which takes into account how cells divide. With respect18 to individual growth rate, a fraction of the ingested food is ∧ channelled to maintenance and growth. The necessity of main-19 tenance could make this quantity negative. Therefore, the animal shrinks in fact (sea-anemones, flatworms, water fleas,20 ∧ etc.) [2–4].21 ∗Corresponding author. Tel.: +34 983 423000x5835; fax: +34 983 423661. E-mail addresses: [email protected] (O. Angulo), [email protected] (J.C. López-Marcos), [email protected] (M.A. López-Marcos). http://dx.doi.org/10.1016/j.camwa.2014.04.009 0898-1221/©2014 Elsevier Ltd. All rights reserved. 2O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx Theoretical analysis of these kinds of problems is highly difficult. With respect to models with a nonlinear dependence 1 on weighted populations, the analysis of existence and uniqueness was proven in [5]. Calsina et al. [6] made an in depth 2 theoretical study in two settings: a finite and infinite size range, although it was made only in the case of the nonlinear au3 tonomous case. It was pointed out in [4], that for an analysis of models that involves a dynamical resource we had to perform 4 as carried out in [6]. The study of this theoretical framework appeared first in a Thieme’s presentation http://math.la.asu. 5 edu/∼dieter/workshop/schedule.html. Nevertheless, a problem with only positive growth was analysed in [7]. They showed 6 stability properties and bifurcation phenomena within a study of the renewal equation for the consumer population birth 7 rate coupled to a delay differential equation for the resource concentration. We can find an extensive study of physiolog8 ically structured population models, with analytical studies of aspects such as existence and uniqueness, smoothness and 9 the asymptotic behaviour of solutions in [8,2,9,10,4]. 10 Often, models such as those discussed above cannot be solved analytically and require numerical integration to obtain an 11 approximation to the solution. Nevertheless, standard numerical methods cannot be indiscriminately applied, because they 12 could lead to inaccurate results and, therefore, wrong conclusions. We could introduce the following examples. We showed 13 that the application of a non A-stable numerical method for the simulation to this problem in [11] was not suitable for the 14 approximation of its asymptotically stable steady state. Some spurious oscillations occur if we use a wrong choice of the 15 numerical integrator [12]. The choice of the selection procedure makes the approximation to singular asymptotically stable 16 steady states better [13]. Also, the lack of regularity of the solution plays an essential role in using an ad hoc method [14]. 17 The numerical solution to the model, due to its obvious mathematical complexity, entails a serious challenge. De Roos [3]18 introduced a semidiscrete method, the ‘‘escalator boxcar train’’, in terms of momenta of the original density function. How19 ever, it did not provide a direct approximation to the density function and its convergence has not yet been considered yet. 20 In [15], we consider a direct integration of the system by means of a version of a numerical method, successfully employed 21 in [11], to obtain this missing information. Nevertheless, it was shown that some numerical methods were not appropri22 ate for a long-time integration. In that work, we presented a new suitable numerical method based on the modified Euler 23 method and the mid-point quadrature rule. On the other hand, a modification of such a numerical procedure was introduced 24 in [13] in order to approximate singular asymptotic states for these kinds of models. Details about the numerical integration 25 of physiologically-structured population models can be found in [16] and the references therein. 26 On the other hand, numerical methods for the solution to these kinds of models have been successfully applied to struc27 tured models to replicate the available field and/or laboratory data for a variety of different systems: rotifers [17,18], where 28 we showed the existence of an asymptotic stable equilibrium state and the existence of a stable periodic solution with ad hoc 29 schemes due to the lack of regularity of the problem; slugs [19], in which we studied equilibrium and oscillatory solutions 30 of a general mass structured system with a boundary delay: the numerical calculations revealed oscillations, pulse solutions 31 and irregular dynamics; marine invertebrates [20,21], where we approximated accurately the steady states and ∧ analysed 32 the asymptotic behaviour of the solutions to the linear model and provided original knowledge about the mechanisms that 33 govern the stability of a nonlinear system with a dynamical larvae behaviour; forest dynamics [22,23], in which we described 34 coexistence mechanisms in a size-structured model in terms of competitive differences at the regeneration state, etc. 35 Finally, when we have to face a numerical simulation in a problem, we must carry out the analysis of the following 36 numerical properties: consistency, stability and convergence. These properties guarantee the goodness of the method to 37 approach the solution. 38 In this paper, we have developed a new, more efficient second-order numerical method for the problem and performed its 39 complete convergence analysis. We have developed numerical simulations for an academic problem to confirm numerically 40 the convergence order. Also, we have applied it to a significant biological example: the dynamics of a Daphnia magna popu41 lation. It has been studied numerically, but the convergence analysis for the numerical integration has not yet finished [3]. 42 Nevertheless, an analysis of a different scheme that utilized an intermediate value to perform the numerical integration was 43 initiated in [24]. Here, we place emphasis on the approximation to the ∧ asymptotically stable states of the model. In Section 244 we proceed to present the general model. Section 3is devoted to introducing the numerical method employed to approx45 imate the solution to the model. The analysis of the convergence properties is shown in Sections 4–7 and the numerical 46 results performed complete ∧ the last section.47 2. The model 48 We consider the following nonlinear size-structured population model where the population feedback on the individuals 49 life history is given by an integro-ordinary differential equation, 50 ut+(g(x,S(t), t)u)x= −µ(x,S(t), t)u,0<x<xM(t), t>0,(2.1) 51 g(0,S(t), t)u(0,t)=xM(t) 0 α(x,S(t), t)u(x,t)dx,t>0,(2.2) 52 u(x,0)=u0(x), 0≤x≤xM(0), (2.3) 53 S′(t)=f(S(t), I(t), t), t>0,(2.4) 54 S(0)=S0.(2.5) 55 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 3 The model involves a nonlinear hyperbolic partial differential equation which represents a balance law, a nonlinear and a1 nonlocal boundary condition which reflects the reproduction process and an initial size-distribution for the population. This2 size-structured problem is coupled with an initial value one.3 The independent variables xand trepresent size and time, respectively. The dependent variable u(x,t)is the size-specific4 density of individuals with size xat time t. The influence of the environment is given by a function S(t)which represents a5 physiological resource. Its dynamics is driven with the initial value problem (2.4)–(2.5), which is coupled with (2.1)–(2.3).6 The evolution of such a resource also depends on the population, which is performed by means of the nonlocal term I(t)7 defined by8 I(t)=xM(t) 0 γ (x,S(t), t)u(x,t)dx,t≥0.(2.6)9 Functions g,αand µrepresent the growth, fertility and mortality rates, respectively. These are usually called the vital10 functions and define the life history of an individual. Functions αand µare nonnegative. Note that all the vital functions11 (g,µand α) depend on the size x(the structuring internal variable), on time tand on the value of the resource at time t,12 which can reflect the influence of the environmental changes on the vital functions. The function fon the right-hand side13 of (2.4) depends on the value of the resource at time t, on the total amount of individuals in the population by means of the14 weighted functional I(t)(which represents the way of weighting the size distribution density in order to model the different15 influence of the individuals of different sizes on such dynamics) and on time t. In general, the size of any individual varies16 according to the following ordinary differential equation17 dx dt =g(x,S(t), t). (2.7)18 Therefore, as the growth rate does not have a determined sign, the maximum size xM(t)that an individual could have at19 time t, changes with time and its dynamics is described by20 d dt xM(t)=g(xM(t), S(t), t), t>0, xM(0)=xM. (2.8)21 Model (2.1)–(2.5) represents a way of describing the dynamics of a Daphnia magna population into an algal food source and22 was introduced first in [3].23 With regard to the theoretical and numerical ∧ analyses of the model, they were mentioned in the Introduction. We note24 that the mathematical analysis of the solutions to the problem (2.1)–(2.5) in its general setting is beyond the scope of this25 paper and we will concentrate on the numerical aspects of the problem. Thus, we assume the existence and uniqueness of26 sufficiently smooth solutions under suitable hypotheses.27 In the present paper, we introduce a new, more efficient second-order numerical method for this nonlinear model. We28 also validate it with its convergence analysis. Throughout the paper we assume the following regularity conditions on the29 data functions and the solution to the problem (2.1)–(2.6):30 (H1) u∈C2([0,xM(t)]×[0,T]),u(x,t)≥0, x∈ [0,xM(t)],t≥0.31 (H2) S∈C2([0,T]),S(t)≥0, t≥0.32 (H3) γ∈C2([0,xM(t)] × D× [0,T]), where Dis a compact neighbourhood of33 {S(t),0≤t≤T}.34 (H4) µ∈C2([0,xM(t)] × D× [0,T]), is nonnegative and Dis a compact neighbourhood of35 {S(t),0≤t≤T}.36 (H5) α∈C2([0,xM(t)] × D× [0,T]), is nonnegative and Dis a compact neighbourhood of37 {S(t),0≤t≤T}.38 (H6) f∈C2(D×DI× [0,T]), is nonnegative, Dis a compact neighbourhood of39 {S(t),0≤t≤T},40 and DIis a compact neighbourhood of41 xM(t) 0 γ (x,S(t), t)u(x,t)dx,0≤t≤T.42 (H7) g∈C3([0,xM(t)] × D× [0,T]), where Dis a compact neighbourhood of43 {S(t),0≤t≤T}.44 and g(0,S,t)≥C>0, t≥0, S∈R. In addition, the characteristic curves x(t;t∗,x∗)defined in (2.7) are continuous45 and differentiable with respect to the initial values (t∗,x∗)∈ [0,T]×[0,xM(t)].46 The above hypotheses may be based on three possible reasons. First, biological assumptions such as the nonnegativity of47 some of the vital functions or, in (H7), to reflect that any individual in the studied population could shrink [2]. Second, the 4O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx mathematical requirements to obtain the existence and uniqueness of solutions for the problem (2.1)–(2.6) [2]. Finally, the 1 regularity properties needed in the numerical analysis to derive optimal rates of convergence [15]. 2 3. Numerical integration 3 The numerical method we employ to approximate the solution to (2.1)–(2.6) is based on the discretization of a repre4 sentation of the solution along the characteristic curves [11]. First of all we rewrite the partial differential equation (2.1) in 5 a more suitable form for its numerical treatment. So we define 6 µ∗(x,z,t)=µ(x,z,t)+gx(x,z,t). 7 Thus, Eq. (2.1) has the form 8 ut+g(x,S(t), t)ux= −µ∗(x,S(t), t)u,0<x<xM(t), t>0.(3.9) 9 We denote by x(t;t∗,x∗)the characteristic curve of Eq. (3.9) that takes the value x∗at time t∗. Such a characteristic curve 10 is the solution to the initial value problem 11 d dt x(t;t∗,x∗)=g(x(t;t∗,x∗), S(t), t), t≥t∗, x(t∗;t∗,x∗)=x∗. (3.10) 12 Now, we consider the function that represents the solution to (3.9) along the characteristic curves 13 w(t;t∗,x∗)=u(x(t;t∗,x∗), t), t≥t∗,14 which satisfies the initial value problem 15 d dt w(t;t∗,x∗)= −µ∗xt;t∗,x∗,S(t), tw(t;t∗,x∗), t≥t∗,16 w(t∗;t∗,x∗)=u(x∗,t∗), 17 and, therefore, it can be represented in the following integral form 18 w(t;t∗,x∗)=u(x∗,t∗)exp −t t∗ µ∗xτ;t∗,x∗,S(τ), τdτ,t≥t∗.(3.11) 19 Given a constant step k>0, we introduce the discrete time levels tn=n k,n=0,1,2, . . .. We also take Ja positive integer, 20 as a parameter related to the size variable which describes the number of points in the uniform initial grid. The diameter of 21 such a mesh grid is h=xM/J, and the initial grid nodes are X0 j=j h, 0 ≤j≤J. In order to start the integration, we consider 22 as an approximation to the density at the initial time (t0), the grid restriction of the initial condition in (2.3),U0 j=u0(X0 j),23 0≤j≤J. Also, we use S0in (2.5) as the initial value of the resource. Then, the numerical method provides, at each dis24 crete time level, a mesh grid on the size interval in Xn, the approximation to the density on such a mesh grid in Unand the 25 approximation of the value of the resource Sn, from the approximations we computed at the previous time level, by using 26 discretizations of Eqs. (3.10),(3.11),(2.2),(2.5),(2.6). Thus, for n=0,1,2,...,the numerical solution at time tn+1=tn+k,27 is obtained from the known values of the numerical solution at time tnas follows, 28 Xn+1 0=0,(3.12) 29 Xn+1 j+1=Xn j+k 2g(Xn j,Sn,tn)+g(Xn+1,∗ j+1,Sn+1,∗,tn+1),0≤j≤J,(3.13) 30 Sn+1=Sn+k 2fSn,Q(Xn,γn·Un), tn+fSn+1,∗,Q(Xn+1,∗,γn+1,∗·Un+1,∗,tn+1),(3.14) 31 Un+1 j+1=Un jexp −k 2µ∗Xn j,Sn,tn+µ∗Xn+1,∗ j+1,Sn+1,∗,tn+1,(3.15) 32 0≤j≤J,33 Un+1 0=Q(Xn+1,αn+1·Un+1) g(0,Sn+1,tn+1),(3.16) 34 where we have to compute approximations at time level tn+1,35 Xn+1,∗ 0=0,Xn+1,∗ j+1=Xn j+k g(Xn j,Sn,tn), 0≤j≤J,(3.17) 36 Sn+1,∗=Sn+k f Sn,Qn(Xn,γn·Un), tn,(3.18) 37 Un+1,∗ j+1=Un jexp −kµ∗Xn j,Sn,tn,0≤j≤J,(3.19) 38 Un+1,∗ 0=Q(Xn+1,∗,αn+1,∗·Un+1,∗) g(0,Sn+1,∗,tn+1).(3.20) 39 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 5 Note that the general step of the method increases the number of grid points and also the dimension of the vector with the1 numerical densities: at time tn, we have J+1 grid nodes in Xnand the (J+1)-dimensional vector Un, and at time tn+1we2 obtain J+2 grid nodes in Xn+1and the (J+2)-dimensional vector Un+1. In order to maintain the number of grid points3 suitable to perform the next step, we eliminate at time tn+1the first grid node Xn+1 lwhich satisfies4 |Xn+1 l+1−Xn+1 l−1| = min 1≤j≤J+1|Xn+1 j+1−Xn+1 j−1|.(3.21)5 We reproduce the same reduction in the corresponding vector Un+1. However, we do not recompute the approximations to6 the nonlocal terms at such a time level.7 In the description of the method, we use the following notation; vectors αpand γpcontain the evaluations of the functions8 αand γin (2.2) and (2.6), respectively, at the grid points in Xp, at the resource value Spand at time tp. Products γp.Up 9 and αp.Upmust be considered componentwise. In order to approximate integrals over the interval [0,xM(tp)], we use the10 composite trapezoidal quadrature rule based on the grid points Xp= [Xp 0,Xp 1,...,Xp J], that is11 Q(Xp,Vp)= J  j=1 Xp j−Xp j−1 2Vp j−1+Vp j.(3.22)12 Note that the method is implicit: all the expressions provide explicit equations for the numerical values at the highest time13 level, except those which involve the numerical density Up 0at the first grid point, but it is easy to implement the method in14 an explicit form.15 4. Convergence analysis: preliminaries16 Below, we will ∧ analyse numerical methods based on integration along characteristics that use a general quadrature rule17 with suitable properties to approximate the integral terms. Convergence will be obtained by means of consistency and18 nonlinear stability. We use the discretization framework developed by López-Marcos et al. [25].19 We assume that the spatial discretization parameter, h, takes values in the set H= {h>0:h=xM/J,J∈N}. Now, we20 suppose that the time step, k, satisfies k=r h, where ris an arbitrary and positive constant fixed throughout the analysis.21 In addition, we set N= [T/k]. For each h∈H, we define the space22 Ah= N  n=0RJ+n×RJ+n+1×RN+1,23 where RJ+nis used for the approximations to the interior grid nodes and RJ+n+1for the approaches to the theoretical solution24 on them and on the left boundary node, at time level tn, 0 ≤n≤N; and RN+1is employed for the approximations to the25 theoretical solution to the initial value problem. We also consider the space26 Bh=RJ×RJ+1×R×RN× N  n=1RJ+n×RJ+n×RN,27 where RJ×RJ+1×Ris employed to compare with the initial approximations; RNconsiders the residuals that take place28 on the boundary node for every time step; and N n=1RJ+n×RJ+n, is used for the residuals which arise in the formulae29 that define the grid nodes and the solution values; and RNconsiders the residuals that are computed for the resource. Both30 spaces have the same dimension.31 In order to measure the size of the errors, we define32 ∥η∥∞=max 1≤j≤p|ηj|,η=(η1, η2, . . . , ηp)∈Rp,33 34 ∥Vn∥1= J+n  j=0 h|Vn j|,Vn∈RJ+n+1.35 Thus, we endow the spaces Ahand Bhwith the following norms.36 If y0,V0,...,yN,VN,a∈Ah, then37 ∥y0,V0,...,yN,VN,a∥Ah=max max 0≤n≤N∥yn∥∞,max 0≤n≤N∥Vn∥∞,∥a∥∞.38 On the other hand, if Y0,Z0,A0,Z0,Y1,Z1,...,YN,ZN,A∈Bh, thus39 ∥Y0,Z0,A0,Z0,Y1,Z1,...,YN,ZN,A∥Bh= ∥Y0∥∞+ ∥Z0∥∞+ |A0|+∥Z0∥∞ 40 + N  n=1 k∥Zn∥∞+ N  n=1 k∥Yn∥∞+ N  n=1 k|An|.41 6O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx Now, for each h∈H, we define 1 xh=(x0,x1,x2,...,xN), xn=(xn 1,...,xn J+n)∈RJ+n,2 x0 j=j h,1≤j≤J,3 xn j=x(tn;tn−1,xn−1 j−1), 1≤j≤J+n,1≤n≤N;(4.23) 4 and we denote xn 0=0, n≥0. Recall that x(t;t∗,x∗)represents the theoretical solution to problem (2.7),t∗∈ [0,T],x∗∈5 [0,xM(t)]. In addition, if urepresents the theoretical solution to (2.1)–(2.6) we define 6 uh=(u0,u1,u2,...,uN), un=(un 0,un 1,...,un J+n)∈RJ+n+1,7 un j=u(xn j,tn), 0≤j≤J+n,0≤n≤N.(4.24) 8 Finally, if Sis the theoretical solution to (2.5) then we define 9 sh=(s0,s1,s2,...,sN), sn=S(tn), 0≤n≤N.(4.25) 10 Therefore, ˜ uh=(x0,u0,x1,u1,...,xN,uN,sh)∈Ah.11 Next, we introduce the discretization operator. Let Rbe a positive constant and we denote by BAh(˜ uh,R hp)⊂Ahthe 12 open ball with centre ˜ uhand radius R hp, 1 <p<2, 13 Φh:BAh(˜ uh,R hp)→Bh,14 Φhy0,V0,...,yN,VN,a=Y0,P0,A0,P0,Y1,P1,...,YN,PN,A,(4.26) 15 defined by the following equations: 16 Y0=y0−X0∈RJ,(4.27) 17 P0=V0−U0∈RJ+1,(4.28) 18 A0=a0−S0∈R.(4.29) 19 Vectors X0,U0and value S0represent approximations at t=0, respectively, to the initial grid nodes, to the theoretical 20 solution at these nodes and to the initial resource. Also, 21 Pn+1 0=Vn+1 0− Qyn+1,αn+1·Vn+1 g0,an+1,tn+1,(4.30) 22 Yn+1 j+1=1 kyn+1 j+1−yn j−k 2g(yn j,an,tn)+g(yn+1,∗ j+1,an+1,∗,tn+1),(4.31) 23 Pn+1 j+1=1 kVn+1 j+1−Vn jexp −k 2µ∗yn j,an,tn+µ∗yn+1,∗ j+1,an+1,∗,tn+1,(4.32) 24 0≤j≤J+n−1, 25 An+1=1 kan+1−an−k 2fan,Q(yn,γn·Vn), tn+fan+1,∗,Q(yn+1,∗,γn+1,∗·Vn+1,∗), tn+1,(4.33) 26 0≤n≤N−1 ∧ , where, with the notation introduced in Section 3,27 yn+1,∗ j+1=yn j+k g(yn j,an,tn), (4.34) 28 Vn+1,∗ j+1=Vn jexp −kµ∗yn j,an,tn,(4.35) 29 0≤j≤J+n−1, 30 Vn+1,∗ 0=Q(yn+1,∗,αn+1,∗·Vn+1,∗) g(0,an+1,∗,tn+1),(4.36) 31 an+1,∗=an+k f an,Q(yn,γn·Vn), tn,(4.37) 32 0≤n≤N−1. We denote by Q(X,V)=M l=0ql(X)Vl, the general quadrature rule employed in (4.30)–(4.37). We have 33 to highlight that the number of nodes considered at each time level is J+n+1, which does not coincide with the number 34 of nodes of Xn, 0 ≤n≤N. This is due to the fact that the quadrature rules use the fixed values of the nodes Xn 0=xn 0=0, 35 0≤n≤N. This notation is also valid if we consider quadrature rules whose nodes are, at each time level, a subgrid of Xn,36 0≤n≤N.37 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 7 Note that, Φhtakes into account all the possible nodes and their corresponding solution values at each time level, and it1 employs quadrature rules possibly based on a subgrid. If ˜ Uh=(X0,U0,X1,U1,...,XN,UN,sh)∈BAh(˜ uh,R hp), satisfies2 Φh(˜ Uh)=0∈Bh,(4.38)3 the nodes and the corresponding values of the solution at such nodes of ˜ Uhare a numerical solution to the scheme defined4 by (3.13)–(3.16) when the composite trapezoidal quadrature rule is given. On the other hand, the numerical solution to the5 scheme defined by (3.13)–(3.16) satisfies (4.38).6 Henceforth, Cwill denote a positive constant, independent of h,k(k=r h), j(0 ≤j≤J+n) and n(0 ≤n≤N); Cmay7 have different values in different places.8 We assume that the quadrature rule satisfies the following properties:9 (P1) |I(tn)−Q(xn,γn·un)|≤C h2, when h→0, 0 ≤n≤N.10 (P2) xM(t) 0α(x,S(tn), tn)u(x,tn)dx −Q(xn,αn·un)≤C h2, when h→0, 0 ≤n≤N.11 (P3) |qj(xn)| ≤ q h, where qis a positive constant independent of h,k,j(0 ≤j≤J+n) and n(0 ≤n≤N), for 0 ≤j≤J+n,12 0≤n≤N.13 (P4) Let Rand pbe positive constants with 1 <p<2. The quadrature weights qjare Lipschitz continuous functions on14 B∞(xn,R hp), 0 ≤j≤J+n, 0 ≤n≤N.15 (P5) Let Rand pbe positive constants with 1 <p<2. If yn,zn∈B∞(xn,R hp),Vn∈B∞(un,R hp)and an∈B∞(Sn,R hp),16 then17  J+n  i=0qi(yn)−qi(zn)γ (zn i,an,tn)Vn i≤C∥yn−zn∥∞,18 when h→0, 0 ≤n≤N.19 (P6) Let Rand pbe positive constants with 1 <p<2. If yn,zn∈B∞(xn,R hp),Vn∈B∞(un,R hp)and an∈B∞(Sn,R hp),20 then21  J+n  i=0qi(yn)−qi(zn)αzn i,an,tnVn i≤C∥yn−zn∥∞,22 when h→0, 0 ≤n≤N.23 The following result establishes that the composite trapezoidal rule used in our experiments satisfies these properties [15].24 Theorem 1. Assume that the hypotheses (H1)–(H7) hold. If the quadrature rules are the composite trapezoidal quadrature on25 subgrids xn jn lM(n) l=0 ,0≤n≤N with the property26 (SR) There exists a positive constant C such that, for h sufficiently small, xn jn l+1 −xn jn l ≤C h, 0≤l≤M(n)−1, xn jn 0 =0, xn jn M(n) =xJ+n,27 with xn jn lM(n)−1 l=1 contained in xn,0≤n≤N.28 Then, properties (P1)–(P6) hold.29 The following result shows that operator (4.26) is well defined.30 Proposition 1. Assume that hypotheses (H1)–(H7) hold and that the quadrature rules used in (4.30)–(4.37) satisfy the properties31 (P1)–(P6). If32 X0,V0,...,XN,VN,S∈BAh(˜ uh,R hp),33 where R is a fixed positive constant and 1<p<2, then, for h sufficiently small,34 Q(Xn,γn·Vn)∈DI,(4.39)35 0≤n≤N. Furthermore, as h →0,Xn,∗∈B∞(xn,R′hp), Sn,∗∈B∞(sn,R′hp)and Vn,∗∈B∞(un,R′hp), and36 Q(Xn,∗,γn,∗·Vn,∗)∈DI,(4.40)37 1≤n≤N.38 Proof. The definition of Q, the hypotheses (H1)–(H7), the properties (P1)–(P6) and that Vnis bounded, allow us to obtain39 Q(Xn,γn·Vn)−I(tn)≤Q(Xn,γ(Xn,Sn)·Vn)−Q(xn,γ(Xn,Sn)·Vn) 40 +Q(xn,γ(Xn,Sn)·Vn)−Q(xn,γ(Xn,S(tn)) ·Vn) 41 +Q(xn,γ(Xn,S(tn)) ·Vn)−Q(xn,γn·un)+Q(xn,γn·un)−I(tn) 42 8O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx ≤ J+n  j=0qj(Xn)−qj(xn)γ (Xn j,Sn)Vn j+ J+n  j=0 qj(xn)γ (Xn j,Sn)−γ (xn j,sn)Vn j1 + J+n  j=0 qj(xn) γ (xn j,sn)Vn j−un j+o(1)2 ≤C R hp(J+n+1)∥γ∥∞∥Vn∥∞+C q R hp+1(J+n+1)∥Vn∥∞3 +R q hp+1(J+n+1)∥γ∥∞+o(1), (4.41) 4 0≤n≤N,h→0. Therefore (4.39) holds, for hsufficiently small. On the other hand, (4.40) is derived following the same 5 arguments. 6 5. Consistency 7 We define the local discretization error as 8 lh=Φh(˜ uh)∈Bh,9 and we say that the discretization (4.26) is consistent if, as h→0, 10 lim ∥Φh(˜ uh)∥Bh=lim ∥lh∥Bh=0.11 The following theorem establishes the consistency of the numerical scheme defined by Eqs. (4.27)–(4.32).12 Theorem 2. Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Then, 13 as h →0, the local discretization error satisfies, 14 ∥Φh(˜ uh)∥Bh= ∥u0−U0∥∞+ ∥x0−X0∥∞+ |s0−S0| + O(h2+k2). (5.1) 15 Proof. We denote Φh(˜ uh)=(Z0,L0, σ0,L0,Z1,L1,...,ZN,LN,σ).16 First, we set the next bounds for the auxiliary values. Then, by means of regularity hypotheses (H1)–(H7), properties (P1), 17 (P2) of the quadrature rule and error bounds for the explicit Euler method and the rectangular quadrature rule, we obtain 18 |xn j−xn,∗ j| = xtn;tn−1,xn−1 j−1−xn−1 j−1−k g(xn−1 j−1,sn−1,tn−1)19 ≤C k2,(5.2) 20 1≤j≤J+n,21 |sn−sn,∗| ≤ S(tn)−sn−1−k f sn−1,I(tn−1), tn−122 +C k fsn−1,I(tn−1), tn−1−ftn−1,sn−1,Q(xn−1,γn−1·un−1)23 ≤Ck2+kI(tn−1)−Q(xn−1,γn−1·un−1)24 ≤C(k2+h2), (5.3) 25 and 26 |un j−un,∗ j| = |un−1 j−1|exp −tn tn−1 µ∗x(τ;tn−1,xn−1 j−1), S(τ), τdτ27 −exp −kµ∗xn−1 j−1,sn−1,tn−1≤C k2,(5.4) 28 1≤j≤J+n, 1 ≤n≤N. By means of inequality (5.3) we arrive at 29 |un 0−un,∗ 0| ≤ Cg0,sn,∗,tnun 0−Qxn,∗,αn,∗un,∗30 ≤Cg0,sn,∗,tn−g0,sn,tn|un 0|31 +xM(tn) 0 α(x,sn,tn)u(x,tn)dx −Qxn,αnun32 +Qxn,αnun−Qxn,∗,αn,∗un,∗33 ≤C(h2+k2)+Qxn,αnun−Qxn,∗,αn,∗un,∗.(5.5) 34 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 9 Now, by means of inequalities (5.2)–(5.5), hypotheses (H1)–(H7) and the properties of the quadrature rule (P3) and (P6), we1 obtain2 Qxn,αnun−Qxn,∗,αn,∗un,∗≤ J+n  j=0qj(xn)−qj(xn,∗)αxn,∗ j,sn,∗,tnun,∗ j 3 + J+n  j=0qj(xn)αxn,∗ j,sn,∗,tn−αxn j,sn,tnun,∗ j 4 + J+n  j=0qj(xn) α xn j,sn,tnun j−un,∗ j 5 ≤C(h2+k2)+C h Qxn,αnun−Qxn,∗,αn,∗un,∗,6 thus, for hsufficiently small, we arrive at7 Qxn,αnun−Qxn,∗,αn,∗un,∗≤C(h2+k2). (5.6)8 Therefore, by means of (5.5) and (5.6), it follows that9 |un 0−un,∗ 0| ≤ C(h2+k2). (5.7)10 Finally, using the inequalities (5.2)–(5.4) and (5.6)–(5.7), the hypotheses (H1)–(H7), the properties of the quadrature11 rule (P3) and (P5) in (5.6), we obtain12 Qxn,γnun−Qxn,∗,γn,∗un,∗≤C(h2+k2). (5.8)13 Now, we set the bounds for Zn, 1 ≤n≤N. By means of (3.10) and (4.31), the regularity hypotheses (H1)–(H7), inequal-14 ities (5.2)–(5.3) and the error bound of ∧ R–Kschemes employed, we have15 |Zn j| ≤ 1 kxn j−xn−1 j−1−k 2gxn−1 j−1,sn−1,tn−1+gxn j,sn,tn 16 +k 2gxn j,sn,tn−gxn,∗ j,sn,tn 17 +k 2gxn,∗ j,sn,tn−gxn,∗ j,sn,∗,tn 18 ≤Ck2+xn j−xn,∗ j+ |sn−sn,∗| 19 ≤C(k2+h2), (5.9)20 1≤j≤J+n, 1 ≤n≤N.21 Next, arguments analogous to those used to derive (5.9) lead us to establish the bound for the truncation errors produced22 by the resource solution. By means of (4.33), the regularity hypotheses (H1)–(H7), the property of the quadrature rule (P1),23 inequalities (5.3) and (5.8) and the error bound of ∧ R–Kschemes employed, we arrive at24 |σn| ≤ 1 ksn−sn−1−k 2fsn−1,I(tn−1), tn−1+fsn,I(tn), tn 25 +k 2fsn−1,I(tn−1), tn−1−fsn−1,Q(xn−1,γn−1·un−1), tn−1 26 +k 2fsn,I(tn), tn−fsn,∗,I(tn), tn 27 +k 2fsn,∗,I(tn), tn−fsn,∗,Q(xn,∗,γn,∗·un,∗), tn, 28 ≤Ck2+sn−sn,∗+I(tn−1)−Q(xn−1,γn−1·un−1) 29 +I(tn)−Q(xn,γn·un)+Q(xn,γn·un)−Q(xn,∗,γn,∗·un,∗) 30 ≤Ck2+h2,(5.10)31 1≤n≤N. Finally, we establish the bounds for the truncation errors produced by the solution to the PDE with a similar32 reasoning. By means of (3.11) and (4.32), the regularity hypotheses (H1)–(H7), the property of the quadrature rule (P1),33 16 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 1≤n≤N. And, finally, we substitute (6.38)–(6.40) in (6.27),(6.29) and (6.31) to arrive at 1 ∥En∥∞≤C∥∆0∥∞+ ∥E0∥1+ |σ0|+∥P0−R0∥∞+ n−1  m=1 k|Am−Bm|2 + n  m=1 k∥Pm−Rm∥∞+ n  m=1 k∥Ym−Zm∥∞,(6.41) 3 1≤n≤N. So, due to (6.38),(6.39) and (6.41) we have 4 ∥∆0,E0,...,∆N,EN,σ∥Ah≤C∥(∆0,E0, σ0,P0−R0,Y1−Z1,P1−R1,...,YN−ZN,PN−RN,A−B)∥Bh.5 7. Convergence 6 The global discretization error is defined as 7 ˜ eh=˜ uh−˜ Uh∈Ah.8 We say that the discretization (4.26) is convergent if there exists h0>0 such that, for each h∈Hwith h≤h0,(4.38) has a 9 solution ˜ Uhfor which, as h→0, 10 lim ∥˜ uh−˜ Uh∥Ah=lim ∥˜ eh∥Ah=0.11 In our analysis, we shall use the following result of the general discretization framework introduced by López-Marcos 12 et al. [25]. 13 Theorem 4. Assume that (4.26) is consistent and stable with thresholds Rh. If Φhis continuous in B(˜ uh,Rh)and ∥lh∥Bh=o(Rh)14 as h →0, then: 15 (i) For h sufficiently small, the discrete equations (4.38) possess a unique solution in B(˜ uh,Rh).16 (ii) As h →0, the solutions converge and ∥˜ eh∥Ah=O(∥lh∥Bh).17 Finally, we propose the next theorem which establishes the convergence of the numerical method defined by Eqs. (4.27)– 18 (4.32).19 Theorem 5. Assume that hypotheses (H1)–(H7) hold and that the considered quadrature rules satisfy properties (P1)–(P6). Then, 20 for h sufficiently small, the numerical method defined by Eqs. (4.27)–(4.32) has a unique solution ˜ Uh∈B(˜ uh,Rh)and 21 ∥˜ Uh−˜ uh∥Ah≤C∥x0−X0∥∞+ ∥u0−U0∥∞+ |s0−S0| + O(h2+k2).(7.1) 22 The proof of Theorem 5 is immediately derived by means of consistency (Theorem 2), stability (Theorem 3) and Theorem 4.23 Specifically, if X0=x0,U0=u0and S0=s0, the proposed numerical scheme is second-order accurate. 24 At this moment, we have obtained convergence of the numerical method (4.27)–(4.33) which does not employ selection 25 at each time level. Also, we have proven the convergence of numerical methods which employ a selection criterion, when26 ever the positions, which are determined by the criterion we have chosen, lead us to subgrids which satisfy property (SR). 27 For the criterion presented in this paper, this property may be shown in two stages. First, as proven in [15], it leads us to 28 subgrids with such a property when we applied it over nodes which are in a neighbourhood of the theoretical ones with 29 radius R hp. In a second stage, it is proven that the nodes, which in fact the numerical method computes, are in such neigh30 bourhoods. In order to do this, it is enough to realize that such nodes could be seen, up to each time level, as the solutions 31 obtained by a discrete operator which has the form of that defined in (4.26).32 Remark. It is well known that regularity hypotheses are necessary, in numerical analysis, to derive optimal rates of 33 convergence for numerical quadrature rules and numerical methods for differential equations. This is the meaning of our 34 smoothness assumptions. However, an interesting question is to consider how the numerical scheme analysed in this paper 35 should be used to carry out the numerical integration of problems with non-smooth biological data functions. In such case, 36 sometimes it is possible to locate these singular points, thus we would obtain the convergence result taking into account 37 that the method is based on the approximation along the characteristics curves by means of quadrature rules. For example, 38 with respect to the discontinuities caused by the lack of compatibility among the initial and boundary conditions, the 39 discussion introduced in [15], for the model without dynamical resource, is also valid in this case. In other physiologically 40 structured models, we can observe finite jump discontinuities caused by the problem itself (as in [14]), or by non-smooth 41 coefficients (as in [26,27]) which propagate along the characteristic curves. While these points are located, a ∧ proper choice 42 of the mesh grid and a suitable simple adaptation of the quadrature rules, keep the order of convergence of the method. The 43 same studies could be made for our model if they would be necessary. Also we have found situations, as in [13], providing 44 unbounded solutions close to the maximum size. Again, we can perform a different selection procedure and a modification 45 of the quadrature rule to describe the approximation to the solution. Finally, we refer to [28,29], and the references therein, 46 for numerical approximation to weak solutions for a similar kind of models. 47 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 17 Table 1 Error and experimental order of convergence. k\h2.5e−2 1.25e−2 6.25e−2 3.13e−2 1.56e−2 2.5e−2 1.565536e−4 2.674563e−4 2.852092e−4 2.888894e−4 2.897345e−4 1.25e−26.434510e−4 3.741278e−5 6.469671e−5 6.920614e−5 7.013337e−5 2.07 2.05 2.04 2.04 6.25e−27.647515e−4 1.656411e−4 9.157412e−6 1.590974e−5 1.704878e−5 1.96 2.03 2.02 2.02 3.13e−28.238862e−4 1.916386e−4 4.209129e−5 2.265792e−6 3.944845e−6 2.00 1.98 2.01 2.01 1.56e−28.438404e−4 2.051189e−4 4.790993e−5 1.060032e−5 5.635547e−7 2.01 2.00 1.99 2.01 8. Numerical results1 We have carried out different numerical experiments with the scheme defined in Section 3.2 First, we have considered an academical problem. It consists of a theoretical test problem with meaningful nonlinearities3 (both from a mathematical and abiological point of view). In this case, the numerical integration for the experiment was4 carried out on the time interval [0,10]. The size interval was taken as [0,1]. Below, we describe the functions involved in5 the experiment. The size-specific growth, fertility and mortality moduli are chosen as6 g(x,z,t)=λ 2 1+z zz 1+z2 −x2+xr 1+z29 30 −z k,7 α(x,z,t)=3 2λ 1+z C29 30 −z k−29λ 30r 1+2z C29 30 −z k−29λ 30r ,8 µ(x,z,t)=λ1+z zz 1+z+2x−3r 1+z29 30 −z k.9 The weight function is taken as γ (x,z,t)=x2and, finally,10 f(z,i,t)=rz 1−z k−rzi (1+z)5 z5(1+4e−λt).11 With this choice of data functions, the problem (2.1)–(2.6) has the following solution12 u(x,t)=S(t) 1+S(t)−x2 −e−λtS(t) 1+S(t)2 −x2,13 S(t)=29 30 Ce29rt/30 1+Ce29rt/30/k,14 r=0.1, C=24, k=5, λ=0.3. Since we know the exact solution to the problem, we can show numerically that our15 method is second-order accurate by means of an error table. In Table 1, each entry in columns two to seven of the upper16 value represents the global error17 eh,k=max max 0≤j≤J|u(X0 j,t0)−U0 j|,|S(t0)−S0|,max 1≤n≤Nmax 0≤j≤J|u(Xn j,tn)−Un j|,|S(tn)−Sn| 18 and the lower number is the experimental order sof the method as computed from19 s=log(e2h,2k/eh,k) log 2 .20 Each column and each row of the ∧ table correspond to different values of the spatial and time discretization parameter, re-21 spectively. The results in the ∧ table clearly confirm the expected second-order convergence. The property of convergence22 18 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx in finite time interval is important to carry out experiments in which the long-time behaviour of the population is investi1 gated. Therefore, the numerical integration of the model with an efficient method allows us to consider a more realistic test 2 problem. 3 Second, the numerical method has been employed to describe the dynamics of a population of ectothermic invertebrates: 4 the water flea, Daphnia magna. In this particular case, the data functions employed in [3] are given by 5 g(x,z,t)=gz 1+z−x,6 µ(x,z,t)=µ, 7 8 α(x,z,t)=αz 1+zx2,9 f(z,i,t)=r z 1−z K−iz 1+z,10 γ (x,z,t)=x2.11 Simple calculations provide the following steady state of the model,Q2 12 u∗(x)=αr g(1+S∗)1−S∗ K1−x x∗ Mµ g−1 ,0≤x≤x∗ M,(8.1) 13 with the following maximum size and resource 14 x∗ M=3 µ(µ +g)(µ +2g) 2αg2,S∗=x∗ M 1−x∗ M .(8.2) 15 This model, with the following set of parameter values µ=0.1, α=0.75, r=3, K=8.3 and xM=1, was ∧ analysed 16 with different numerical techniques in [3,13]. This set is beyond the theoretical analysis settings because the steady state is 17 unbounded (although integrable in whatever case g> µ). 18 We have performed simulations with g=0.0075, µ=0.1, α=0.75, r=3 and xM=1. This set of parameters 19 makes the solution be bounded and therefore, within the theoretical analysis settings. In a first experiment, we take the 20 value K=8.3. We employ initial conditions (2.3) and (2.5) which ensure compatibility between u0,S0and the problem, 21 u0(0)=α gxM 0 x2u0(x)dx,S0 1+S0=3g u0(0) 2αxM 0x u0(x)dx +g u′ 0(0),22 and integrate numerically until T=1000. We employ small perturbations of (8.1) and (8.2) as initial conditions: xM=0.875, 23 S0=7, and 24 u0(x)=     αr g(1+S0)1−S0 K1−x  xMβ if 0 ≤x≤ xM, 0 if  xM<x≤xM. (8.3) 25 The value of βis taken in order to ensure the compatibility conditions 26 (β +3)(β +2)(β +1)=2α x3 M g,27 which results in β≈0.5153057420. 28 The computations performed with different values of the parameters of the discretization, kand J, show that the numer29 ical solution is attracted to a stationary state. Therefore, we can state that equilibrium is asymptotically stable. For a fixed 30 number of grid points on the size interval, J, the numerical steady state converges to the theoretical one as the time step k31 decreases. However, we observe that a critical value of the time step k, depending on J, appears. Below this critical value, 32 the numerical method does not produce better approximations to the theoretical values of the equilibrium state. Our expe33 rience allows us to estimate that an optimal ratio of the discretization parameters is about k J ≈250. In Table 2, we present 34 the errors produced by the numerical approximations to the theoretical equilibrium state with different values of the dis35 cretization parameters (Jand k). Note that the convergence is second-order. For this choice of the final time integration, the 36 approach to the theoretical steady state is good enough as shown in Table 2.37 In Fig. 1, we show the long time evolution of the numerical solution obtained with the values of the discretization pa38 rameters k=0.0625, J=4000. In the ∧ figure, we draw the evolution of the density function, the total population computed 39 from the population density by the composite trapezoidal quadrature rule, the maximum individual size and the dynamical 40 resource. As we can observe, all of them evolve towards the equilibrium state. 41 Taking into account the good behaviour of the numerical solution in a predictable situation, we now consider a more 42 complicated situation. As the value of the parameter Kincreases, the equilibrium state becomes unstable. In a second O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 19 Table 2 Errors in the approximations with respect to the theoretical values: N∗ 0=x∗ M 0u∗(x)dx ≈46.6763867098 (total population), S∗≈4.0859721214 (resource). k J Total population Resource 1 250 4.946335e−3 4.112717e−2 5e−1 500 1.004377e−3 9.992129e−3 2.5e−1 1 000 2.270716e−4 2.461613e−4 1.25e−1 2 000 5.403338e−5 6.108879e−4 6.25e−2 4 000 1.318184e−5 1.521058e−4 3.125e−2 8 000 3.255560e−6 3.788916e−5 1.5625e−2 16 000 8.089576e−7 9.394259e−6 7.8125e−3 32 000 2.016265e−7 2.277898e−6 Fig. 1. Evolution of the numerical solution. Case of a stable steady state. experiment, we present an example of this situation by taking K=9.64. As an initial condition, we take the unstable1 equilibrium state (8.1)–(8.2) (we extend this function to the whole interval [0,xM]by taking u0(x)=0 if x∗ M<x≤xM).2 As in the previous experiment, the values of the discretization parameters are k=0.0625, J=4000. We can observe3 the instability of the equilibrium and the solution evolving towards a cycled situation (Fig. 2). Taking into account that the4 numerical solution is attracted to a limit cycle, considering a sufficiently large time, the numerical solution obtained after5 this long time integration lies practically on such a cycle. In this way, the numerical method provides an approximation to6 the limit cycle. In Fig. 3, the representation of such a cycle in the tridimensional space defined by the total population, the7 maximum individual size and the dynamical resource, is drawn. From the numerical results obtained for the total population,8 the maximum size and the dynamical resource, we can obtain an in depth analysis of these quantities throughout a period9 of the limit cycle. For example, in Table 3, we present the behaviour of some numerical quantities of the solution ∧ throughout10 a period of the limit cycle. Also, we can estimate such a period by interpolation: it is about 64.6824.11 9. Conclusions12 We have ∧ analysed a second-order numerical method for a problem that describes a population with a possible shrinking13 size and with a dependency on the environment managed by the evolution of a vital resource. The method involves only one14 level of time each step. The second-order convergence has been theoretically proven by means of an argument of consistency15 and stability of the scheme. The academical test problem allows us to report numerical experiments which demonstrate the16 predicted accuracy of the scheme.17 This knowledge leads us to make a long time integration of the model. The biological problem we consider is well known.18 It describes the dynamics of a Daphnia magna population. With its long time integration, we observe the good stability19 20 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx Fig. 2. Evolution of the numerical solution. Case of an unstable steady state. Fig. 3. Limit cycle appearing in the case of an unstable steady state. properties that our numerical scheme possesses. In the settings of the theoretical analysis, our numerical scheme makes it 1 possible to determine the rich dynamics of the model. It presents an equilibrium that is stable until the parameter Kreaches 2 the bifurcation value. After that, it becomes unstable. In the stable case, we show how the numerical steady state approaches 3 the theoretical steady state. When the equilibrium is unstable, a stable limit cycle appears, the characteristics of which are 4 also described. 5 Acknowledgments 6 We thank the reviewers for their constructive and helpful suggestions. This work was supported in part by the project 7 of the Ministerio de Ciencia e Innovación MTM2011-25238 and the project of the Consejería de Educación, JCyL VA191U13. 8 O. Angulo et al. / Computers and Mathematics with Applications xx (xxxx) xxx–xxx 21 Table 3 Behaviour of some numerical quantities of the solution throughout a period of the limit cycle. Time Resource Total population Maximum size 0 Abs. maximum (5.91) ↗ ↗ 14.22 ↘ ↗ Abs. maximum(0.82) 19.57 ↘Rel. maximum(53.41) ↘ 27.01 Rel. minimum(3.42) ↘ ↘ 30.90 ↗ ↘ Rel. minimum(0.80) 32.15 ↗Rel. minimum(51.37) ↗ 36.21 Rel. maximum(4.92) ↗ ↗ 43.09 ↘ ↗ Rel. maximum(0.81) 45.92 ↘Abs. maximum(53.52) ↘ 57.50 Abs. minimum(2.55) ↘ ↘ 60.20 ↗ ↘ Abs. minimum(0.77) 61.79 ↗Abs. minimum(47.35) ↗ 64.68 Abs. maximum(5.92) ↗ ↗ References1 [1] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, A semi-lagrangian method for a cell population model in a dynamical environment, Math. Comput. Modelling 57 (2013) 1860–1866. 2 [2] J. A. J. Metz, E. O. Diekmann (Eds.), The Dynamics of Physiologically Structured Populations, in: Springer Lecture Notes in Biomathematics, vol. 68, Springer, Heildelberg, 1986. 3 [3] A. M. de Roos, Numerical methods for structured population models: the escalator boxcar train, Numer. Methods Partial Differential Equations 4 (1988) 173–195. 4 [4] B. Perthame, Transport Equations in Biology, Birkhäuser Verlag, Basel, 2007.5 [5] N. Kato, A general model of size-dependent population dynamics with nonlinear growth rate, J. Math. Anal. Appl. 297 (2004) 234–256.6 [6] A. Calsina, J. Saldaña, A model of physiologically structured population dynamics with a nonlinear individual growth rate, J. Math. Biol. 33 (4) (1995) 335–364. 7 [7] O. Diekmann, M. Gyllenberg, J.A.J. Metz, S. Nakaoka, A.M. de Roos, Daphnia revisited: local stability and bifurcation theory for physiologically structured population models explained by way of an example, J. Math. Biol. 61 (2010) 277–318. 8 [8] G.F. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985.9 [9] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, in: Applied Mathematics Monographs, C.N.R., Giardini Editori e Stampatori, Pisa, 1995. 10 [10] J.M. Cushing, An Introduction to Structured Populations Dynamics, in: CMB-NSF Regional Conference Series in Applied Mathematics, SIAM, 1998.11 [11] L. Abia, O. Angulo, J. C. López-Marcos, M. A. López-Marcos, Long-time simulation of a size-structured population model with a dynamical resource, Math. Model. Nat. Phenom. 5 (6) (2010) 1–21. 12 [12] O. Angulo, A. Durán, J.C. López-Marcos, Numerical study for size-structured population models: A case of Gambussia affinis, C. R. Biol. 328 (2005) (2005) 387–402. 13 [13] O. Angulo, J. C. López-Marcos, M. A. López-Marcos, Numerical approximation of singular asymptotic states for a size-structured population model with a dynamical resource, Math. Comput. Modelling 54 (2011) 1693–1698. 14 [14] L. Abia, O. Angulo, J. C. López-Marcos, M. A. López-Marcos, Numerical schemes for a size-structured cell population model with equal fission, Math. Comput. Modelling 50 (2009) 653–664. 15 [15] O. Angulo, J.C. López-Marcos, Numerical integration of fully nonlinear size-structured population models, Appl. Numer. Math. 50 (2004) 291–327.16 [16] L.M. Abia, O. Angulo, J.C. López-Marcos, Size-structured population models and their numerical solutions, Discrete Contin. Dyn. Syst. Ser. B 4 (2004) 1203–1222. 17 [17] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, A numerical simulation for the dynamics of the sexual phase of monogonont rotifera, C. R. Biol. 327 (2004) 293–303. 18 [18] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, A numerical integrator for a model with a discontinuous sink term: the dynamics of the sexual phase of monogonont rotifera, Nonlinear Anal. RWA 6 (2005) 935–954. 19 [19] O. Angulo, J.C. López-Marcos, M.A. Bees, Mass structured systems with boundary delay: oscillations and the effect of selective predation, J. Nonlinear Sci. 22 (2012) 961–984. 20 [20] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, J. Martínez-Rodríguez, Numerical investigation of the recruitment process in open marine population21 models, J. Stat. Mech. (2011) P01003. http://dx.doi.org/10.1088/1742-5468/2011/01/P01003.22 [21] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, J. Martínez-Rodríguez, Numerical analysis of a population model of marine invertebrates with different life stages, Commun. Nonlinear Sci. Numer. Simul. 18 (2013) 2153–2163. 23 [22] M.A. Zavala, O. Angulo, R. Bravo de la Parra, J.C. López-Marcos, An analytical model of stand dynamics as a function of tree growth, mortality and recruitment: the shade tolerance-stand structure hypothesis revisited, J. Theoret. Biol. 244 (2007) 440–450. 24 [23] O. Angulo, R. Bravo de la Parra, J.C. López-Marcos, M.A. Zavala, Stand dynamics and tree coexistence in an analytical structured model: the role of recruitment, J. Theoret. Biol. 333 (2013) 91–101. 25 [24] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, Numerical analysis of a size-structured population model with a dynamical resource, Biomath (2013) submitted. Q326 [25] J. C. López-Marcos, J. M. Sanz-Serna, Stability and convergence in numerical analysis III: linear investigation of nonlinear stability, IMA J. Numer. Anal. 8 (1988) 71–84. 27 [26] O. Angulo, J.C. López-Marcos, M.A. López-Marcos, A numerical integrator for a model with a discontinuous sink term: the dynamics of the sexual phase of monogonont rotifera, Nonlinear Anal. RWA 6 (2005) 935–954. 28 [27] M.A. Bees, O. Angulo, J.C. López-Marcos, D. Schley, Dynamics of a structured slug population model in the abscence of seasonal variation, Math. Models Methods Appl. Sci. 16 (2006) 1961–1985. 29 [28] J. Shen, C. Shu, M. Zhang, High resolution schemes for a hierarchical size-structured model, SIAM J. Numer. Anal. 45 (2007) 352–370.30 [29] A.S. Ackleh, B. Ma, J.J. Thibodeaux, A second-order high resolution finite difference scheme for a structured erythropoiesis model subject to malaria infection, Math. Biosci. 245 (2013) 2–11. 31