Numerical integration of a hierarchically size-structured population model with contest competition

Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx Contents lists available at ScienceDirect Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam Numerical integration of a hierarchically size-structured population model with contest competition Q1 ∧ L.M. ∧ Abiaa,∧ O. ∧ Angulob,∗,∧ J.C. ∧ López-Marcosa,∧ M.A. ∧ López-Marcosa aDepartamento de Matemática Aplicada. Facultad de Ciencias. Universidad de Valladolid. Valladolid, Spain bDepartamento de Matemática Aplicada. ETS de Ingenieros de Telecomunicación. Universidad de Valladolid. Paseo de Belén, 15 - Campus Miguel Delibes, 47011 Valladolid, Spain article info Article history: Received 21 June 2011 Received in revised form 9 July 2013 Keywords: Hierarchically size-structured population models Numerical integration Convergence abstract We formulate schemes for the numerical solution to a hierarchically size-structured population model. The schemes are analysed and optimal rates of convergence are derived. Some numerical experiments are also reported to demonstrate the predicted accuracy of the schemes and to show their behaviour to approaching stable steady states. ©2013 Published by Elsevier B.V. 1. Introduction1 The study of population dynamics has evolved from simple unstructured models to fairly sophisticated structured ones.2 In these, there exists a variable which structures the population into different classes. The most widely used is age, but it is3 difficult to measure in a large number of species. Physiological characteristics, grouped generically under the name of size,4 are the alternative to age. These structured population models involve a first-order hyperbolic partial differential equation5 which represents a conservation law for the population, a boundary condition which reflects the reproduction process and6 an initial condition given by the initial size-distribution of the population. The evolution of the population is determined by7 so-called ‘‘vital functions’’ (fertility, mortality and growth rates). Reciprocally, the population also has an influence on such8 evolution through the vital functions which depend on functionals of the population density, so that the model becomes a9 nonlinear one. We can find more details about physiologically structured models in [1–4], about their numerical integration10 we refer to [5,6] and about the computational study of real populations. See, for example, [7–11].11 The model we study in this work assumes that the influence of the population on the life-history of ∧ individuals is of a12 contest kind. This means that the individuals in the population do not have the same opportunity in the competition for13 resources. This kind of model is known as a hierarchical size-structured population one. Here, we study a standard problem.14 It consists of a nonlinear partial differential equation (the population balance law),15 ut+(g(x,B(x,t), t)u)x= −µ(x,B(x,t), t)u,xm<x<xM,t>0,(1.1)16 a nonlinear and nonlocal boundary condition which represents the birth law,17 g(xm,B(xm,t), t)u(xm,t)=Γ(t)+xM xm α(x,B(x,t), t)u(x,t)dx,t>0,(1.2)18 ∗Corresponding author. Tel.: +34 983 423000X5680; fax: +34 983 423 661. E-mail addresses: [email protected] (L.M. Abia), [email protected] (O. Angulo), [email protected] (J.C. López-Marcos), [email protected] (M.A. López-Marcos). 0377-0427/$ – see front matter ©2013 Published by Elsevier B.V. http://dx.doi.org/10.1016/j.cam.2013.09.005 2L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx and an initial condition, 1 u(x,0)=φ(x), xm≤x≤xM,(1.3) 2 where the nonlocal term B(x,t), is defined by 3 B(x,t)=xM x γ (σ) u(σ , t)dσ , xm<x<xM,t>0.(1.4) 4 The independent variables xand trepresent size and time, where xmand xMare, respectively, the minimum and maximum 5 value reached by an individual in a given population. The function u(x,t)is the size-specific population density of individuals 6 with size xat time t. We also assume that the size of any individual varies according to the following ordinary differential 7 equation 8 d dt x(t)=g(x(t), B(x(t), t), t). (1.5) 9 Population dynamics ∧ are determined by the growth rate g, the mortality rate µ, the reproduction rate αand the external 10 inflow Γ. These nonnegative vital functions (growth, mortality and reproduction rates) depend on the structuring variable 11 and on the functional B(x,t)used to describe the competition among individuals. In this case, no individual in a class of a 12 smaller size can affect the amount of available resources of an individual larger in size, which is a kind of contest competition. 13 Also, all the functions depend on the time ∧ variable in order to take into ∧ account environmental changes. We could consider 14 different functionals in the vital functions dependency (g, µ and β). It would increase the complexity of the numerical 15 method description. However, this more general case would not increase the difficulty of the analysis. 16 This model can be used to describe the dynamics of a size-structured tree population which takes into account the effect 17 of competition for light [12–14]. In this case, the size is given by the diameter at breast height (d.b.h.) and B(x,t)represents 18 the cumulative basal area of trees with a size greater than xand expresses the shading effect under light competition. 19 From a theoretical point of view, Cushing [15] was the first to study a model like this. It was structured by age and included 20 a more general type of nonlinearity in which the vital rates depended on the number of older and/or younger individuals 21 than itself. However, the approach was limited to vital rates with no explicit age dependence. Later, in [16], Calsina and 22 Saldaña considered a hierarchically size-structured model, but they did not include explicit dependence on size in the vital 23 rates. They obtained the existence, uniqueness, asymptotic behaviour of the solutions and global asymptotical stability of 24 a nontrivial equilibrium. The existence and uniqueness of solutions for an autonomous model which included the explicit 25 dependence on the structured variable were studied by Kraev [17]. And Ackleh and Ito, in [18], proved the existence of 26 measure-valued solutions in a more general setting with less restrictions for the growth rate. 27 We should point out here that, without other restrictive assumptions, this kind of model cannot be solved analytically. 28 Moreover, when such models include nonlinearities and environmental dependence on the different vital rates, the use of 29 efficient methods that provide a numerical approach is the most suitable mathematical tool for studying the problem and, 30 indeed, is often the only one available. Nevertheless, the numerical approach to these equations has important drawbacks 31 because they are usually nonlinear equations, and the nonlinearities of the PDE and nonlocal boundary condition are caused 32 by nonlocal terms. The model (1.1)–(1.4) was studied in [19], where an implicit first-order finite difference scheme was 33 ∧ analysed and its stability and convergence, as well as the existence, uniqueness, and well-posedness (in the L1norm) of 34 bounded variation weak solutions for (1.1), were proved. In [20], for the autonomous version of (1.1)–(1.4), Shen et al. 35 proved the convergence of an upwind scheme and another method which were only first-order in time and second-order 36 in size discretization. 37 In the present paper, we introduce second-order numerical methods for the solution to this nonlinear and 38 nonautonomous model and we also develop the convergence analysis of a general formulation. We carry out the numerical 39 integration of Eqs. (1.1)–(1.4) by means of methods which integrate along the characteristic curves in different settings: the 40 use of all the grid nodes, the use of a grid with a constant number of nodes in every time step by means of a suitable selection, 41 and the use of a more specific selection of grid nodes which allows us to change their number (increasing or diminishing). 42 The integral terms are approximated by means of second-order quadrature rules. We analyse the consistency, stability and 43 convergence properties of a general numerical scheme which covers the situation where the quadrature rule used at a time 44 step is based on nodes obtained by integration along the characteristics with the nodes employed at the previous time step. 45 We also provide numerical experiments to show the accuracy of the schemes and their behaviour in long time integration. 46 These methods could be clearly generalized for other situations like the ones studied in [21]. Also, we would like to point 47 out that the nonautonomous case is important from the biological point of view, as we can see in [22] and the references 48 therein. 49 Throughout the paper we assume the following regularity conditions on the data functions and the solution to the 50 problem (1.1)–(1.4) which will be integrated in a fixed time interval [0,T]:51 (H1) u∈C2([xm,xM]×[0,T]), u(x,t)≥0,x∈ [xm,xM],u(xM,t)=0,t≥0. 52 (H2) γ∈C2([xm,xM]).53 (H3) Γ∈C1([0,T]).54 (H4) µ∈C2([xm,xM]×DB×[0,T]), is nonnegative. 55 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 3 Table 1 Notation. Model xmMinimum size xMMaximum size gGrowth rate BCompetition functional µMortality rate γWeight function αFertility rate φInitial density Γexternal inflow DBCompact neighbourhood [0,T]Fixed time interval Numerical integration x(t;t∗,x∗)Characteristic curve mAuxiliary rate w(t;t∗,x∗)Solution along characteristic curve NNumber of time steps J+1 Number of initial grid nodes hSpatial discretization parameter kTime discretization parameter tnDiscrete time level Jn+1 Number of Grid nodes at tnXn jGrid node Un jNumerical Approximation of u(Xn j,tn)In jNumerical approximation to B(Xn j,tn) QjQuadrature rule Xn,∗ jAuxiliary node Un,∗ jAuxiliary solution approximation In.∗ jAuxiliary functional approximation βFixed value in selection strategy (2.20) λnnumber of eliminated nodes Convergence HSet of hvalues rFixed positive constant AhSpace of approaches BhSpace of residuals xhTheoretical grid values uhTheoretical solution values R,pPositive constants 8hDiscretization operator yn,znRelated to grid values qn,j iQuadrature weight Vn,WnRelated to solution values λn l(j)First subgrid subindex greater than j Yn,ZnRelated to residuals in grid approach lhLocal discretization error Pn,LnRelated to residuals in solution approach ˜ ehGlobal discretization error (H5) α∈C2([xm,xM]×DB×[0,T]), is nonnegative.1 (H6) g∈C3([xm,xM] × DB× [0,T]), is nonnegative, g(xm,z,t)≥C>0,g(xM,z,t)=0, and gz(x,z,t)≤0,x∈2 [xm,xM],z∈R,t>0.3 In hypotheses (H4)–(H6), DBis a compact ∧ neighbourhood of4 {B(x,t),xm≤x≤xM,0≤t≤T}.5 The above hypotheses can be based on three possible reasons. First, biological assumptions such as the nonnegativity of6 the vital functions or, in (H6), to reflect that the size of any individual in the population studied [23,1,24] is strictly increasing7 during its lifetime and never reaches the maximum value. Second, the mathematical requirements to obtain the existence8 and uniqueness of solutions to the problem (1.1)–(1.4), [21]. Last, the regularity properties needed in the numerical analysis9 to derive optimal rates of convergence [25].10 The paper is structured as follows. In Section 2, we introduce the numerical methods we employ to make the integration11 of the problem (1.1)–(1.4). From Sections 3 to 6we analyse the numerical schemes: first we introduce some preliminary12 results, then we prove the consistency and stability properties in order to obtain their convergence. In the final section we13 present numerical results which confirm the theoretical ones and the good behaviour of each numerical scheme when it14 approaches the stable steady state of a particular problem. The notation employed throughout the paper is given in Table 1.15 16 2. Numerical integration17 In this section, we describe the numerical methods employed for the integration of Eqs. (1.1)–(1.4). First, we rewrite the18 partial integro-differential equation in a more suitable form for its numerical treatment. Thus, we define19 m(x,z,u,t)=µ(x,z,t)+gx(x,z,t)−gz(x,z,t) γ (x)u.20 Thus, (1.1) has the form21 ut(x,t)+g(x,B(x,t), t)ux(x,t)= −m(x,B(x,t), u(x,t), t)u(x,t), (2.6)22 xm<x<xM,t>0. Next, we denote by x(t;t∗,x∗)the characteristic curve of Eq. (2.6) which takes the value x∗at time t∗.23 This is the solution to the following initial value problem,24 x′(t;t∗,x∗)=g(x(t;t∗,x∗), B(x(t;t∗,x∗), t), t), t≥t∗, x(t∗;t∗,x∗)=x∗.(2.7)25 4L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx Note that x(t;0,xM)=xM,t≥0; because we assume that g(xM,·,·)=0. Then, we define the function 1 w(t;t∗,x∗)=u(x(t;t∗,x∗), t), t≥t∗,(2.8) 2 which satisfies the following initial value problem, 3    d dt w(t;t∗,x∗)= −mxt;t∗,x∗,B(x(t;t∗,x∗), t), w(t;t∗,x∗), tw(t;t∗,x∗), t≥t∗, w(t∗;t∗,x∗)=u(x∗,t∗), (2.9) 4 and, therefore, it can be represented by the following formula 5 w(t;t∗,x∗)=u(x∗,t∗)exp −t t∗ mxτ;t∗,x∗,B(xτ;t∗,x∗, τ), w(τ;t∗,x∗), τdτ.(2.10) 6 We suppose that u(xM,0)=0 and then u(xM,t)=0,t≥0. We shall use this property in our numerical method. However, 7 it can be easily modified to cover other cases. 8 Our numerical methods are based on the numerical integration along characteristics curves because these kind of 9 schemes have shown its efficiency in other cases [25]. Also, other types of methods (like the total variation diminishing 10 schemes) might not be suitable in simple situations in which the total variation of the solution ∧ increases. This situation 11 could ∧ happen when the growth rate has no dependence on the population and gx+µ < 0 and, therefore, the solution could 12 be unbounded. Also, when gzis positive and the conditions of existence of Kraev [17] do not satisfy. 13 Now, we will introduce the numerical schemes we are going to use. The numerical integration along with characteristics 14 means that the number of grid nodes increases ∧ by one at every time step because there is one node which fluxes from the 15 minimum size, as we can see in [25]. Therefore, first we will introduce a numerical method which employs all the grid nodes 16 obtained in the computation. 17 Let Jand Nbe positive integers. We define the spatial and time discretization parameters as h=xM−xm Jand k=T N,18 respectively. We define the discrete time levels as tn=n k,0≤n≤N, and the initial grid nodes as X0 j=xm+j h,0≤j≤J.19 We suppose that the approximations to the theoretical solution at the initial time (φin (1.3)) in such nodes are known, 20 U0 j,0≤j≤J. Thus, we denote 21 X0=X0 0=xm,X0 1,...,X0 J−1,X0 J=xM,U0=U0 0,U0 1,...,U0 J−1,U0 J=0,22 the initial size grid and the initial approximation to the solution. Also, we introduce I0=I0 0,I0 1,...,I0 J−1,I0 J=0, where 23 I0 j=Qj(X0,γ(X0)U0), (2.11) 24 is defined as the discrete version of (1.4) at x=X0 j,j=0,1,...,J, by using approximations to the integral term. In (2.11) 25 and henceforth, the product of the vectors γ(X)Umust be interpreted componentwise. In this case, we propose the following 26 second-order composite quadrature rules throughout the paper; for (p+1)-dimensional vectors X={X0,X1,...,XP}and 27 V={V0,V1,...,VP}, representing the vector of nodes of the spatial grid and the vector of values of the function at the 28 spatial grid, respectively, 29 Q0(X,V)=V1(X1−xm)+ P−1  l=1 Xl+1−Xl 2(Vl+Vl+1),(2.12) 30 Qj(X,V)= P−1  l=j Xl+1−Xl 2(Vl+Vl+1),1≤j≤P−1.(2.13) 31 Expressions γ(X)and α(X,U)denote the vectors obtained from the evaluation of function γand αat the corresponding 32 values. 33 We obtain the numerical approximations at time level t1as follows. First, we compute the following auxiliary values 34 X1,∗ 0=xm;X1,∗ j+1=X0 j+k g(X0 j,I0 j,t0), 0≤j≤J−1;X1,∗ J+1=xM;35 U1,∗ j+1=U0 jexp −k m X0 j,I0 j,U0 j,t0,0≤j≤J−1;U1,∗ J+1=0;36 and I1,∗ j=Qj(X1,∗,γ(X1,∗)U1,∗), 1≤j≤J+1. These formulae are valid approximations for the solution to the problem, 37 but they only contribute with the first convergence order and we want to build a second-order convergence scheme. 38 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 5 Then, we obtain the grid nodes X1=X1 0=xm,X1 1,...,X1 J,X1 J+1=xM, by the numerical integration of (2.7) with the1 improved Euler method,2 X1 j+1=X0 j+k 2g(X0 j,I0 j,t0)+g(X1,∗ j+1,I1,∗ j+1,t1),(2.14)3 0≤j≤J−1. Now, we calculate U1=U1 0,U1 1,...,U1 J,U1 J+1=0, the corresponding approximations to the theoretical4 solution. For the approximation to the corresponding interior grid points, we use the following discretization of (2.10), based5 on the trapezoidal quadrature rule,6 U1 j+1=U0 jexp −k 2mX0 j,I0 j,U0 j,t0+mX1,∗ j+1,I1,∗ j+1,U1,∗ j+1,t1,(2.15)7 0≤j≤J−1. Then, we compute I1 j=Qj(X1,γ(X1)U1), 0≤j≤J+1. Finally, we derive the approximation U1 0to u(xm,t1)8 from a discrete version of the boundary condition (1.2) Q29 U1 0=Γ(t1)+Q0(X1,α(X1,U1)U1) g(xm,I1 0,t1).(2.16)10 As in (2.11), we denote the componentwise product of the vectors α(X,U)and Uby α(X,U)U. Note that the method11 introduces a new grid node which fluxes from the boundary. Thus, we set J0=Jand introduce the number of grid nodes at12 level t1as J1+1, with J1=J0+1. We also want to indicate that such a scheme is explicit.13 Next, we describe the general time step tn+1,0≤n≤N−1. Now, we suppose that the numerical approximations at14 the previous time level tnare known,15 Xn 0=xm,Xn 1,...,Xn Jn−1,Xn Jn=xM,Un 0,Un 1,...,Un Jn−1,Un Jn=0,16 and In 0,In 1,...,In Jn−1,In Jn, where Jnis the number of grid nodes at tn. We recall that Xn jand Xn+1 j+1,0≤j≤Jn−1, are17 (numerically) in the same characteristic curve. First, we compute the auxiliary values,18 Xn+1,∗ 0=xm;Xn+1,∗ j+1=Xn j+k g(Xn j,In j,tn), 0≤j≤Jn−1;19 Un+1,∗ j+1=Un jexp −k m Xn j,In j,Un j,tn,0≤j≤Jn−1;20 Xn+1,∗ Jn+1=xM;Un+1,∗ Jn+1=0;21 and In+1,∗ j=Qj(Xn+1,∗,γ(Xn+1,∗)Un+1,∗), 1≤j≤Jn+1. We correct these values to obtain second-order approximations.22 The grid values at the time level tn+1,23 Xn+1=Xn+1 0=xm,Xn+1 1,...,Xn+1 Jn,Xn+1 Jn+1=xM,24 by means of the numerical integration of (2.7),25 Xn+1 j+1=Xn j+k 2g(Xn j,In j,tn)+g(Xn+1,∗ j+1,In+1,∗ j+1,tn+1),(2.17)26 0≤j≤Jn−1, and the approximations to the theoretical solution in these nodes at such a time level,27 Un+1=Un+1 0,Un+1 1,...,Un+1 Jn,Un+1 Jn+1=0,28 using the discretization of (2.10) for the approximation at the interior grid points,29 Un+1 j+1=Un jexp −k 2mXn j,In j,Un j,tn+mXn+1,∗ j+1,In+1,∗ j+1,Un+1,∗ j+1,tn+1,(2.18)30 0≤j≤Jn−1, and the approximation Un+1 0to u(xm,tn+1), using a discretization of the boundary condition (1.2),31 Un+1 0=Γ(tn+1)+Q0(Xn+1,α(Xn+1,Un+1)Un+1) g(xm,In+1 0,tn+1),(2.19)32 where αj(Xn+1,Un+1)=α(Xn+1 j,In+1 j,tn+1), 0≤j≤Jn+1, and33 In+1 j=Qj(Xn+1,γ(Xn+1)Un+1), 0≤j≤Jn+1.34 We define Jn+1=Jn+1. We observe that, at consecutive time levels, we work with a different number of nodes because35 we have introduced a new node which fluxes through the boundary. So, at time level tn, we have (Jn+1)grid nodes and36 at time level tn+1we have (Jn+2). This fact increases the memory requirements and computational cost. We could think37 of another kind of method which, at each time step, selects the nodes used in the numerical integration with different38 efficiency motivations. For example, we could propose a method which keeps the number of nodes constant at every time39 6L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx step. Another possibility consists of selecting the grid nodes to obtain a suitable dynamic of the grid points in large time 1 integration. These kind of schemes need a selection procedure after each time step. 2 Next, we describe the general setting when we introduce such a selection. The first time step integration is given by 3 the Eqs. (2.14)–(2.16) and J1=J+1. The equations of the general time integration tn+1,0≤n≤N−1 are given by 4 (2.17)–(2.19). At this moment, the treatment of the number of nodes in the subsequent levels of time has to be decided. As 5 we pointed out, we could think of several strategies for the number of grid nodes. Therefore, we take the choice in the grid 6 nodes which would continue the numerical integration and we consider Jn+1related to the number of grid nodes we set at 7 time level tn+1. For example, we could select some characteristic curves and we do not compute the approximations at such 8 curves. We define λnas the first l∈Nwhich satisfies 9 |Xn+1 Jn+1−Xn+1 Jn−(l+1)|> β h,(2.20) 10 with βa fixed value. We eliminate the grid nodes Xn+1 Jn−l,l=0, . . . , λn−1. Also, we do not consider the corresponding values 11 in the vectors Un+1and In+1. Then, we define the value Jn+1as Jn−λn+1, where λn=0 means that we do not eliminate 12 any node. Another possible choice involves eliminating the computation of the characteristic curve which begins at Xn+1 l,13 so that 14 |Xn+1 l+1−Xn+1 l−1| = min 1≤j≤Jn|Xn+1 j+1−Xn+1 j−1|(2.21) 15 and the corresponding value in the vectors Un+1and In+1. Then, we define the value Jn+1as Jn. Again, we should point out 16 that the scheme is explicit. 17 3. Convergence analysis: preliminaries 18 In this section, we begin the analysis of a general class of numerical methods. More precisely, we consider a formulation 19 in which, for each time level, all the nodes are considered, as were done in the method described by (2.14)–(2.19). However, 20 at each time level, we have a quadrature rule which could be based on a subgrid. In this way, the numerical schemes with a 21 selection procedure could be viewed as methods which employ quadrature rules based on the position of the nodes selected 22 at each time step. When a node is eliminated, the quadrature rules considered for the following time steps do not take into 23 account the position of this node as moving along the characteristic curve. Actually, the numerical methods as presented 24 in the previous section do not compute the ∧ characteristic curves relative to those nodes to improve their efficiency, but to 25 carry out our convergence analysis it is necessary to take into account all of them. 26 The convergence result will be obtained by means of consistency and nonlinear stability. In order to carry out this analysis, 27 we have to rewrite it into the discretization framework developed by López-Marcos et al. [26]. 28 We assume that the spatial discretization parameter, h, takes values in the set H= {h>0:h=(xM−xm)/J,J∈N}.29 Now, we suppose that the time step, k, satisfies k=r h, where ris an arbitrary and positive constant fixed throughout the 30 analysis. In addition, we set N= [T/k]. For each h∈H, we define the space 31 Ah= N  n=0RJ+n−1×RJ+n,32 where, for each time level n,0≤n≤N,RJ+n−1is used for the approximations to the interior grid nodes, and RJ+nfor those 33 to the theoretical solution on them and on the left boundary node. We also consider the space 34 Bh=RJ−1×RJ×RN× N  n=1RJ+n−1×RJ+n−1,35 where RJ−1×RJis employed to compare with the initial approximations; RNconsiders the residuals which take place in 36 the approximation to the solution at the boundary node for every time step; and N n=1RJ+n−1×RJ+n−1, is used for the 37 residuals which arise in the formulae which define the grid nodes and the solution values. We note that in spaces Ahand 38 Bh, we do not consider the first and last grid nodes and the value of the solution at the last grid node because they are fixed 39 values. Thus, both spaces have the same dimension. 40 In order to measure the size of the errors, we define 41 ∥η∥∞=max 1≤j≤p|ηj|,η=(η1, η2, . . . , ηp)∈Rp,42 (B∞(η, ρ) represents the corresponding open ball with centre ηand radius ρ > 0), and 43 ∥Vn∥1= J+n−1  j=0 h|Vn j|,Vn∈RJ+n.44 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 7 Now, we endow spaces Ahand Bhwith the following norms. If y0,V0,...,yN,VN∈Ah, then1 ∥y0,V0,...,yN,VN∥Ah=max max 0≤n≤N∥yn∥∞,max 0≤n≤N∥Vn∥∞.2 On the other hand, if Y0,Z0,Z0,Y1,Z1,...,YN,ZN∈Bh, then3 Y0,Z0,Z0,Y1,Z1,...,YN,ZN∥Bh= ∥Y0∥∞+∥Z0∥∞+∥Z0∥∞+ N  n=1 k∥Zn∥∞+ N  n=1 k∥Yn∥∞.4 For each h∈H, we define5 xh=(x0,x1,x2,...,xN),6 xn=(xn 1,...,xn J+n−1)∈RJ+n−1,0≤n≤N,7 x0 j=xm+j h,1≤j≤J,8 xn j=x(tn;tn−1,xn−1 j−1), 1≤j≤J+n−1,1≤n≤N;(3.1)9 and we denote xn 0=xmand xn J+n=xM,n≥0. Recall that x(t;t∗,x∗)represents the theoretical solution to problem (1.5),10 t∗∈ [0,T],x∗∈ [xm,xM]. In addition, if urepresents the theoretical solution to (1.1)–(1.4) we define11 uh=(u0,u1,u2,...,uN),12 un=(un 0,un 1,...,un J+n−1)∈RJ+n,0≤n≤N,13 un j=u(xn j,tn), 0≤j≤J+n−1,0≤n≤N,(3.2)14 and we denote un J+n=0. Therefore, ˜ uh=(x0,u0,x1,u1,...,xN,uN)∈Ah.15 Next, we define the discretization operator. Let Rbe a positive constant and we denote by BAh(˜ uh,R hp)⊂Ahthe open16 ball with ∧ centre ˜ uhand radius R hp,1<p<2. Next, we introduce the operator17 8h:BAh(˜ uh,R hp)→Bh,18 8hy0,V0,...,yN,VN=Y0,P0,P0,Y1,P1,...,YN,PN,(3.3)19 defined by the following equations,20 Y0=y0−X0∈RJ,(3.4)21 P0=V0−U0∈RJ+1.(3.5)22 Vectors X0and U0represent approximations at t=0, respectively, to the initial grid nodes and to the theoretical solution23 at them. Also,24 Pn+1 0=Vn+1 0−Γ(tn+1)+Qn+1 0yn+1,α(yn+1,Vn+1)Vn+1 gxm,Qn+1 0yn+1,γ(yn+1)Vn+1,tn+1,(3.6)25 Yn+1 j+1=1 kyn+1 j+1−yn j−k 2gyn j,Qn jyn,γ(yn)Vn,tn +gyn+1,∗ j+1,Qn+1 j+1yn+1,∗,γ(yn+1,∗)Vn+1,∗,tn+1,(3.7)26 Pn+1 j+1=1 kVn+1 j+1−Vn jexp −k 2myn j,Qn jyn,γ(yn)Vn,Vn j,tn +myn+1,∗ j+1,Qn+1 j+1yn+1,∗,γ(yn+1,∗)Vn+1,∗,Vn+1,∗ j+1,tn+1,(3.8)27 0≤j≤J+n−1,0≤n≤N−1. Where, with the notation introduced in Section 2,28 yn+1,∗ 0=xm,yn+1,∗ J+n+1=xM,29 yn+1,∗ j+1=yn j+k g(yn j,Qn jyn,γ(yn)Vn,tn), 0≤j≤J+n−1,(3.9)30 Vn+1,∗ j+1=Vn jexp −k m yn j,Qn jyn,γ(yn)Vn,Vn j,tn,0≤j≤J+n−1,31 Vn+1,∗ J+n+1=0,(3.10)32 8L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 0≤n≤N−1, αn j(y,V)=α(yn j,Qn j(y,γ(y)V),tn), γ n j(y)=γ (yn j), 0≤j≤J+n,1≤n≤N;γn,∗ j(y)=γ (yn,∗ j), 0≤1 j≤J+n+1. On the other hand, we denote by 2 Qn 0(X,V)= J+n  l=1 qn,0 l(X)Vl,Qn j(X,V)= J+n  l=j qn,j l(X)Vl,3 the general quadrature rules employed at each time level. These quadrature rules use fixed values for nodes yn 0=xn 0=4 xm,yn J+n=xn J+n=xM, and for the solution Vn J+n=Vn J+n=0,0≤n≤N.5 Note that, 8htakes into account all the possible nodes and their corresponding solution values at each time level, and it 6 employs quadrature rules possibly based on a subgrid. If ˜ Uh=(X0,U0,X1,U1,...,XN,UN)∈BAh(˜ uh,R hp), satisfies 7 8h(˜ Uh)=0∈Bh,(3.11) 8 the nodes and the corresponding values to the solution at such nodes of ˜ Uhare a numerical solution to the scheme defined 9 by (2.14)–(2.19) when the quadrature rules are given by (2.12)–(2.13). On the other hand, the numerical solution to the 10 scheme defined by (2.14)–(2.19) satisfies (3.11).11 Henceforth, Cwill denote a positive constant, independent of h,k(k=r h), j(0≤j≤J+n)and n(0 ≤n≤N); C12 possibly has different values in different places. Now, we introduce the following properties that we suppose the quadrature 13 rules satisfy. These are the sufficient properties the quadrature rules have to satisfy to carry out our convergence analysis. 14 (P1) B(xn j,tn)−Qn j(xn,γ(xn)un)≤C h2, if h→0, where γj(xn)=γ (xn j), 0≤j≤J+n,0≤n≤N.15 (P2) xM xmα(x,B(x,tn), tn)u(x,tn)dx −Qn 0(xn,α(xn,Bn)un)≤C h2, if h→0,0≤n≤N. Where αn j(xn,Bn)=16 αxn j,B(xn j,tn), tn,0≤j≤J+n,0≤n≤N.17 (P3) |qn,j i(xn)| ≤ q h, where qis a positive constant independent of h,k(k=r h), j(0 ≤j≤J+n) and n(0 ≤n≤N), 18 j≤i≤J+n,0≤j≤J+n,0≤n≤N.19 (P4) Let Rand pbe positive constants with 1 <p<2. The functions qn,j iare Lipschitz continuous in B∞(xn,R hp), j≤i≤20 J+n,0≤j≤J+n,0≤n≤N.21 (P5) Let Rand pbe positive constants with 1 <p<2. If yn,zn∈B∞(xn,R hp)and Vn∈B∞(un,R hp),22 |Qn j(yn,γ(zn)Vn)−Qj(zn,γ(zn)Vn)| ≤ C∥yn−zn∥∞,23 if h→0,0≤j≤J+n,0≤n≤N.24 (P6) Let Rand pbe positive constants with 1 <p<2. If yn,zn∈B∞(xn,R hp)and Vn∈B∞(un,R hp),Q3 25 |Qn j(yn,α(zn,Vn)Vn)−Qn j(zn,α(zn,Vn)Vn)| ≤ C∥yn−zn∥∞,26 if h→0,0≤j≤J+n,0≤n≤N.27 The following result shows that operator (3.3) is well defined. 28 Proposition 1. Let us assume that hypotheses (H1)–(H6)about problem (1.1)–(1.4), and properties (P1)–(P6)of the quadrature 29 rules hold. If 30 X0,V0,...,XN,VN∈BAh(˜ uh,R hp), 31 where R is a fixed positive constant and 1<p<2, then, for h sufficiently small, 32 Qn j(Xn,γ(Xn)Vn)∈DB,(3.12) 33 0≤n≤N, and 34 Qn j(Xn,∗,γ(Xn,∗)Vn,∗)∈DB,(3.13) 35 0≤n≤N−1.36 Proof. The definition of Qn j, the hypotheses (H1)–(H6), the properties (P1)–(P6) and that Vnis bounded, allow us to obtain 37 Qn j(Xn,γ(Xn)Vn)−B(xn j,tn)≤C∥Xn−xn∥∞38 + J+n  l=j qn,j l(X)γXn l−γxn lVl+ J+n  l=j qn,j l(X)γxn lVl−un l+o(1)39 ≤C R hp+C q R hp+1(J+n−j+1)∥Vn∥∞+∥γ∥∞+o(1), (3.14) 40 0≤j≤J+n,0≤n≤N,h→0. Therefore, (3.12) holds for hsufficiently small. On the other hand, (3.13) is easily derived 41 following the same arguments. 42 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 9 On the other hand, the quadrature rules employed in the numerical methods presented in Section 2could be written into1 this notation as2 Qn 0(X,V)=Vλn 1Xλn 1−xm+ Jn−1  l=1 Xλn l+1−Xλn l 2Vλn l+Vλn l+1,(3.15)3 Qn j(X,V)=Vλn l(j)Xλn l(j)−Xj+ Jn−1  l=l(j) Xλn l+1−Xλn l 2Vλn l+Vλn l+1,(3.16)4 1≤j≤J+n−1, where λn l(j)corresponds to the first subindex in the subgrid with j≤λl(j). If we consider this choice of5 quadrature rules, it is shown that they satisfy the properties (P1)–(P6).6 Theorem 1. Let us assume that the hypotheses (H1)–(H6)about problem (1.1)–(1.4) hold. Let R and p be positive constants with7 1<p<2. If yn,zn∈B∞(xn,R hp)and Vn∈B∞(un,R hp)and xn λn lJn l=0,0≤n≤N, are subgrids with the property8 (SG) There exists a positive constant C such that, for h sufficiently small, xn λn l+1−xn λn l≤C h,0≤l≤Jn−1,xn λn 0=xm,xn λn Jn=xn J+n,9 with xn λn lJn l=0 contained in xn,0≤n≤N.10 Thus, the quadrature rules (3.15)–(3.16) satisfy properties (P1)–(P6).11 Proof. We note that properties (P1)–(P4) can be easily derived under our assumptions by means of property (SG) and the12 properties of the composite trapezoidal quadrature rule and the rectangular quadrature rule. On the other hand, proof of13 (P5) and (P6) is very similar. Therefore, we will establish (P5) and omit the proof of (P6).14 We will rewrite the formula in this way15 Qn j(yn,γ(zn)Vn)−Qn j(zn,γ(zn)Vn)=yn λl(j)−zn λl(j)−yn j−zn jγ (zn λl(j))Vn λl(j) 16 − Jn−1  l=l(j)+1 yn λl−zn λl 2γzn λl+1Vn λl+1−γzn λl−1Vn λl−1− yn λl(j)−zn λl(j) 2γzn λl(j)+1Vn λl(j)+1+γzn λl(j)Vn λl(j) 17 + yn λJn−zn λJn 2γzn λJ−1Vn λJ−1+γzn λJVn λJ.(3.17)18 Now, the hypotheses (H1)–(H6), ∥Vn∥∞<∞, and the property (SG), allow us to obtain19 γzn λl+1Vn λl+1−γzn λl−1Vn λl−1≤γzn λl+1−γxn λl+1Vn λl+1 20 +γzn λl−1−γxn λl−1Vn λl−1+γxn λl+1Vn λl+1−un λl+1 21 +γxn λl−1Vn λl−1−un λl−1+γxn λl+1un λl+1−γxn λl−1un λl−1 22 ≤Chp+h,(3.18)23 l(j)+1≤l≤Jn−1. Thus, by means of (3.17)–(3.18) and Jn≤(C1h)−1, we obtain24 |Qn j(yn,γ(zn)Vn)−Qn j(zn,γ(zn)Vn)| ≤ C∥yn−zn∥∞1+hp+h(Jn−1−l(j)) 25 ≤C∥yn−zn∥∞(3.19)26 as desired. 27 4. Consistency28 We define the local discretization error as29 lh=8h(˜ uh)∈Bh,30 and we say that the discretization (3.3) is consistent if, as h→0,31 lim ∥8h(˜ uh)∥Bh=lim ∥lh∥Bh=0.32 The following theorem establishes the consistency of the numerical scheme defined by Eqs. (3.4)–(3.8).33 Theorem 2. Let us assume that hypotheses (H1)–(H6)about problem (1.1)–(1.4), and properties (P1)–(P6)of the quadrature34 rule hold. Then, as h →0, the local discretization error satisfies,35 ∥8h(˜ uh)∥Bh= ∥u0−U0∥∞+∥x0−X0∥∞+O(h2+k2). (4.1)36 16 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx Theorem 5. Let us assume that hypotheses (H1)–(H6)about problem (1.1)–(1.4) hold, and that the considered quadrature rules 1 satisfy properties (P1)–(P6). Then, for h sufficiently small, the numerical method defined by Eqs. (3.4)–(3.8) has a unique solution 2 ˜ Uh∈B(˜ uh,Rh)and 3 ∥˜ Uh−˜ uh∥Ah≤C∥x0−X0∥∞+∥u0−U0∥∞+O(h2+k2).(6.1) 4 The proof of Theorem 5 is derived by means of consistency (Theorem 2), stability (Theorem 3) and Theorem 4.5 Next, we can establish an error bound for the numerical and the theoretical solution at the numerical values of the grid 6 nodes. 7 Theorem 6. Let us assume that hypotheses (H1)–(H6)about problem (1.1)–(1.4) hold, and that the considered quadrature rules 8 satisfy properties (P1)–(P6). For h sufficiently small, let u∗ h=(u0 ∗,u1 ∗,u2 ∗,...,uN ∗)∈N n=0RJ+n, defined by 9 un ∗=u(Xn 0,tn), u(Xn 1,tn), . . . , u(Xn J+n−1,tn)∈RJ+n,10 0≤n≤N, where Xn j,0≤j≤J+n−1,0≤n≤N, are the grid nodes given by scheme (3.4)–(3.8). Then, 11 ∥Un−un ∗∥∞≤C∥x0−X0∥∞+∥u0−U0∥∞+O(h2+k2).(6.2) 12 This theorem follows immediately from Theorem 5. In particular, if X0=x0and U0=u0, the proposed numerical scheme 13 is second-order accurate. 14 At this moment, we have obtained the convergence of the numerical method (2.14)–(2.19) which does not employ 15 selection at each time level. Also, we have proved the convergence of numerical methods which employ a selection ∧ criteria,16 whenever the positions, which are determined by the ∧ criteria we have chosen, lead us to subgrids which satisfy the property 17 (SG). For the criteria presented in this paper, this property could be shown in two stages. First, as proved in [25], for the 18 selection ∧ criteria given in (2.20), it leads us to subgrids with such aproperty, when we applied it over nodes which are in 19 a ∧ neighbourhood of the theoretical ones with radius R hp(for ∧ criteria (2.21) we could prove the same in a similar way). In 20 a second stage, it is proved that the nodes, which in fact the numerical method computes, are in such ∧ neighbourhoods. In 21 order to do this, it is enough to realize that such nodes could be seen, up to each level of time, as the solutions obtained by 22 a discrete operator which has the form of the one defined in (3.3).23 With respect to the convergence behaviour of these kind of numerical schemes when the theoretical solution has weak 24 singularities, we refer to the corresponding section of [25] because for both equations and numerical methods the behaviour 25 along characteristics curves are similar. 26 7. Numerical results 27 We have carried out numerical experiments with the scheme defined in Section 2. We have considered a theoretical test 28 problem that presents meaningful nonlinearities (both from a mathematical and biological point of view). The numerical 29 integration for the numerical experiment was carried out on the time interval [0,10]. The size interval was taken as [0,1].30 The size-specific growth, fertility and mortality moduli are chosen as g(x,z,t)=1−x 1+z, α(x,z,t)=x(1−x) 2(z+1),31 µ(x,z,t)=2(1+z) 1+1 2e−2t(1−x)2+(1−x)(1−x)2+4e2t.32 The weight function is taken as γ (x)=1 (we consider the total population) and the external inflow function as 33 Γ(t)=e−2t 1+1+2e2t √1+4e2t 1+e−2t 2(1+√1+4e2t)+2et−√1+4e2t 2−2e2tlog e−t 2(1+1+4e2t).34 Finally, we consider as the initial size-specific density the function 35 u0(x)=(1−x)+(1−x)2+2 (1−x)2+4.36 With the functions chosen, the problem (1.1)–(1.4) has the following solution 37 u(x,t)=e−2t(1−x)+(1−x)2+2e2t (1−x)2+4e2t.38 Since we know the exact solution to each problem, we can show numerically that our methods are second-order accurate 39 by means of error tables. In Tables 2–4, we present the results obtained for the test problem with the three different selection 40 procedures. 41 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 17 Table 2 Error and experimental order of convergence. Method without selection. k h 3.125e−2 1.563e−2 7.813e−3 3.906e−3 1.953e−3 3.125e−2 1.585e−4 8.297e−5 8.456e−5 8.496e−5 8.506e−5 1.563e−21.580e−4 4.223e−5 2.069e−5 2.109e−5 2.119e−5 1.908 2.004 2.004 2.004 7.813e−31.581e−4 4.033e−5 1.089e−5 5.166e−6 5.265e−6 1.970 1.955 2.002 2.002 3.906e−31.586e−4 3.987e−5 1.019e−5 2.766e−6 1.291e−6 1.987 1.985 1.977 2.001 1.953e−31.590e−4 3.982e−5 1.001e−5 2.560e−6 6.970e−7 1.994 1.994 1.993 1.989 Table 3 Error and experimental order of convergence. Method with selection (2.21) (constant number of grid nodes). k h 3.125e−2 1.563e−2 7.813e−3 3.906e−3 1.953e−3 3.125e−2 1.585e−4 8.085e−5 8.449e−5 8.495e−5 8.506e−5 1.563e−21.701e−4 4.223e−5 2.036e−5 2.107e−5 2.118e−5 1.908 1.989 2.004 2.004 7.813e−31.779e−4 4.229e−5 1.089e−5 5.111e−6 5.260e−6 2.008 1.955 1.994 2.002 3.906e−31.813e−4 4.287e−5 1.047e−5 2.766e−6 1.280e−6 2.053 2.015 1.977 1.997 1.953e−31.830e−4 4.320e−5 1.042e−5 2.599e−6 6.970e−7 2.069 2.040 2.010 1.989 Table 4 Error and experimental order of convergence. Method with selection (2.20) (asymptotic selection), β= 0.125. k h 3.125e−2 1.563e−2 7.813e−3 3.906e−3 1.953e−3 3.125e−2 1.585e−4 8.297e−5 8.456e−5 8.496e−5 8.506e−5 1.563e−21.580e−4 4.223e−5 2.069e−5 2.109e−5 2.119e−5 1.908 2.004 2.004 2.004 7.813e−31.581e−4 4.033e−5 1.089e−5 5.166e−6 5.265e−6 1.970 1.955 2.002 2.002 3.906e−31.586e−4 3.987e−5 1.019e−5 2.766e−6 1.291e−6 1.987 1.985 1.977 2.001 1.953e−31.590e−4 3.982e−5 1.001e−5 2.560e−6 6.970e−7 1.994 1.994 1.993 1.989 In each entry in columns two to seven of Tables 2–4 the upper value represents the global error1 eh,k=max max 0≤j≤J|u(X0 j,t0)−U0 j|,max 1≤n≤Nmax 0≤j≤J+1|u(Xn j,tn)−Un j| 2 and the lower number is the experimental order sof the method as computed from3 s=log(e2h,2k/eh,k) log(2).4 Each column and each row of the table correspond to different values of the spatial and time discretization parameter,5 respectively. The results in the ∧ table clearly confirm the expected second-order of convergence with every selection6 procedure.7 In Fig. 1, we present, ∧ on a logarithmic scale, an efficiency plot where we show the error as a function of the computational8 time (in seconds). We build, for each selection procedure in the numerical scheme, a table with the errors corresponding to9 different values of parameters kand h. Then, we discover the most efficient value of r(for each selection procedure) and,10 18 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx Fig. 1. Error vs. cpu-time. Selection strategy (2.21) (), selection strategy (2.20) (◦), no selection strategy (). finally, we compare all the procedures in Fig. 1. We can observe that the selection procedure in which a node is eliminated 1 at every time step shows the best behaviour. 2 On the other hand, we have considered a test problem which was employed in [14] in a framework in which a forest 3 population is structured by its d.b.h. The minimum and maximum size are xm=1 and xM=51, respectively. The size4 specific growth, fertility and mortality moduli are chosen as g(x,z,t)=7x1−x xM16, α(x,z,t)=0, µ(x,z,t)=0.1. 5 The weight function is taken as γ (x)=1 and the external inflow function as Γ(t)=R, suitable chosen in order to satisfy 6 the first compatibility condition. In this context, there is a constant inflow of newborns. Finally, we consider as the initial 7 size-specific density the function 8 u0(x)=                      5+5x−2 2+13 ,x≤2, 10 +15 x−2 21+x−2 2+30 x−2 23x−2 2−2,2<x≤4, 52−x−2 232+3x−2 2−11+2x−2 2−1,4<x≤6, 0x>6. 9 In this case, we do not know the exact solution to the problem (1.1)–(1.4) but we know that the problem has a stable singular 10 steady state given by 11 u∗(x)=Rexp(−µϕ(x)) g(x),x∈ [xm,xM],12 where the function ϕ(x)represents the time required for a tree to increase d.b.h. from xmto x,13 ϕ(x)=x x0 dσ g(σ).14 This problem presents a significant difficulty because the steady state has a singularity at the maximum size of the 15 population. Similar situations have been dealt with in other related problem [8]. 16 In Fig. 2, we present the results obtained when T=10, with the first selection strategy (2.20) and discretization 17 parameter values k=0.00390625 and h=0.125. On the left-hand side, we show the evolution of the total population 18 over time and, on the right-hand side, the density distribution at t=T. The experiment shows that the stable steady state 19 is reached ∧ in this time period. Clearly, the numerical method is able to obtain it. 20 This test shows that the best selection strategy depends on the problem we wish to solve. Its special feature means that 21 the selection strategy that avoids useless nodes close to the maximum size exhibits the best ∧ behaviour. This is experimentally 22 shown because when we applied the most efficient one (2.21), we also obtain this good approximation to the stable steady 23 state but we need the discretization parameters to take values lower than k=0.0009765625 and h=0.125, which 24 represents a more expensive experiment. Otherwise, oscillations or abnormal finishing of computation occur. These are 25 only due to the type of solution we are approaching that is beyond the required hypotheses in the convergence analysis. 26 We could also perform other tests in which the use of other higher order numerical methods might not be suitable. 27 However, in [5] and for an age-structured model, this lack of efficiency was shown. This is due to the fact that, in this 28 L.M. Abia et al. / Journal of Computational and Applied Mathematics xx (xxxx) xxx–xxx 19 Fig. 2. Left-hand side: Evolution over time of total population computed with the first selection strategy (2.20). Right-hand side: density distribution at t=T. biological framework, these kind of methods need several compatibility relationships between the initial and boundary1 conditions and, in reality, biological data barely satisfies the first compatibility relationship. In [5] we showed that they2 produce oscillations. However, characteristics methods follow perfectly the possible discontinuities in the function or in3 some of its derivatives.4 Acknowledgement5 The authors were supported in part by the Ministerio de Ciencia e Innovación (Spain), project MTM2011-25238.6 References7 [1] J.A.J. Metz, E.O. Dieckmann (Eds.), The Dynamics of Physiologically Structured Populations, in: Springer Lecture Notes in Biomathematics, vol. 68, Springer, Heildelberg, 1986. 8 [2] G.F. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985.9 [3] M. Iannelli, Mathematical Theory of Age-Structured Population Dynamics, in: Applied Mathematics Monographs, C.N.R., Giardini Editori e Stampatori, Pisa, 1994. 10 [4] J.M. Cushing, An Introduction to Structured Population Dynamics, in: CMB-NSF Regional Conference Series in Applied Mathematics, SIAM, 1998.11 [5] L.M. Abia, O. Angulo, J.C. López-Marcos, Age-structured population dynamics models and their numerical solutions, Ecol. Model. 188 (2005) 112–136.12 [6] L.M. Abia, O. Angulo, J.C. López-Marcos, Size-structured population dynamics models and their numerical solutions, Discrete Contin. Dyn. Syst. Ser. B 4 (2004) 1203–1222. 13 [7] M. Adimy, O. Angulo, F. Crauste, J.C. López-Marcos, Numerical integration of a mathematical model of hematopoietic stem cell dynamics, Comput. Math. Appl. 56 (2008) 594–606. 14 [8] 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. 15 [9] 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. 16 [10] 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. 17 [11] 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. 18 [12] T. Kohyama, Size-structured multi-species model of rain forest trees, Funct. Ecol. 6 (1992) 206–212.19 [13] M.A. Zavala, R. Bravo de la Parra, A mechanistic model of tree competition and facilitation for Mediterranean forests: scaling from leaf physiology to stand dynamics, Ecol. Mod. 188 (2005) 76–92. 20 [14] 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. 21 [15] J.M. Cushing, The dynamics of hierarchical age-structured populations, J. Math. Biol. 32 (1994) 705–729.22 [16] A. Calsina, J. Saldaña, Asymptotic behaviour of a model of hierarchically structured population dynamics, J. Math. Biol. 35 (1997) 967–987.23 [17] E.A. Kraev, Existence and uniqueness for height structured hierarchical population models, Nat. Resour. Model. 14 (2001) 45–70.24 [18] A.S. Ackleh, K. Ito, Measure-valued solutions for a hierarchically size-structured population, J. Differential Equations 217 (2005) 431–455.25 [19] A.S. Ackleh, K. Ito, A quasilinear hierarchical size-structured model: well-posedness and approximation, Appl. Math. Optim. 51 (2005) 35–59.26 [20] J. Shen, C.-W. Shu, M. Zhang, High resolution schemes for a hierarchical size-structured model, SIAM J. Numer. Anal. 45 (2007) 352–370.27 [21] A. Calsina, J. Saldaña, A model of physiologically structured population dynamics with a nonlinear individual growth rate, J. Math. Biol. 33 (1995) 335–364. 28 [22] O. Angulo, A. Durán, J.C. López-Marcos, Numerical study of size-structured population models: a case of Gambussia affinis, C. R. Biol. 328 (2005) 387–402. 29 [23] A.J. Fabens, Properties and fitting of the von Bertalanffy growth curve, Growth 29 (1965) 265–289.30 [24] D. Sulsky, Numerical solution of structured population models II. Mass structure, J. Math. Biol. 32 (1994) 491–514.31 [25] O. Angulo, J.C. López-Marcos, Numerical integration of fully nonlinear size-structured population models, Appl. Numer. Math. 50 (2004) 291–327.32 [26] 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. 33