Numerical Analysis of a Size-Structured Population Model with a Dynamical Resource

Original article Biomath 3 (2014), 1403241, 1–11 Bf Volume ░, Number ░, 20░░ BIOMATH ISSN 1314-684X Editor–in–Chief: Roumen Anguelov Bf BIOMATH http://www.biomathforum.org/biomath/index.php/biomath/ Biomath Forum Numerical Analysis of a Size-Structured Population Model with a Dynamical Resource Oscar Angulo∗, J.C. L´ opez-Marcos∗and M.A. L´ opez-Marcos∗ ∗Departamento de Matem´ atica Aplicada Universidad de Valladolid, Valladolid, Spain Emails: [email protected]a.es, [email protected]a.es, [email protected]a.es Received: 19 October 2013, accepted: 24 March 2014, published: 28 May 2014 Abstract—In this paper, we analyze the convergence of a second-order numerical method for the approximation of a size-structured population model whose dependency on the environment is managed by the evolution of a vital resource. Optimal convergence rate is derived. Numerical experiments are also reported to demonstrate the predicted accuracy of the scheme. Also, it is applied to solve a problem that describes the dynamics of a Daphnia magna population, paying attention to the unstable case. Keywords-structured population models; numerical methods; convergence; Daphnia magna I. Introduction Physiologically structured population models are based on the use of one or more attributes that structure the individuals in the population. Size is one of the most natural and important attributes structuring the population for many species: typical examples being fishes and trees. In such species the ability of an individual to obtain the necessary resources to survive and reproduce depends strongly on its size. Structured population models reflect the effect of physiological state of individuals on the population dynamics. In addition, the use of nonlinear structured population models allows us to take into account the effect of competition for natural resources in the structured-specific growth, mortality and fertility rates. We can find an extensive study of physiologically structured population models, with analytical studies of aspects such as existence and uniqueness, smoothness and the asymptotic behaviour of solutions in [1], [2], [3], [4], [5]. In this paper, we consider a size-structured population model nonlinearly coupled with an integroordinary differential equation accounting for substrate consumption and/or product formation. It was introduced first in [6] for modeling a Daphnia magna population. The model involves a nonlinear hyperbolic partial differential equation ut+(g(x,S(t),t)u)x=−µ(x,S(t),t)u,(1) 0<x<xM(t),t>0, a nonlinear and nonlocal boundary condition which reflects the reproduction process g(0,S(t),t)u(0,t)=ZxM(t) 0 α(x,S(t),t)u(x,t)dx,(2) t>0, and an initial size-distribution for the population u(x,0) =u0(x),0≤x≤xM(0).(3) Citation: Oscar Angulo, J.C. L´ opez-Marcos, M.A. L´ opez-Marcos, Numerical Analysis of a Size-Structured Population Model with a Dynamical Resource, Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 1 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... The influence of the environment on the life history of the individuals is given by a function S(t) which represents a physiological resource. Its dynamics is managed by the next initial value problem, S0(t)=f(S(t),I(t),t),t>0,(4) S(0) =S0, which is coupled with (1)-(3). The evolution of the resource also depends on the population, which is performed by means of the nonlocal term I(t) defined by I(t)=ZxM(t) 0 γ(x,S(t),t)u(x,t)dx,t≥0.(5) Also, we consider that the maximum size xM(t) that an individual could have at time t, changes with time and its dynamics is described by d dt xM(t)=g(xM(t),S(t),t),t>0, (6) xM(0) =xM. The independent variables xand trepresent size and time, respectively. The dependent variable u(x,t) is the size-specific density of individuals with size xat time t. In general, the size of any individual varies according to the following ordinary differential equation dx dt =g(x,S(t),t).(7) Functions g,αand µrepresent the growth, fertility and mortality rates, respectively. These are usually called the vital functions and define the life history of an individual. Functions αand µ are nonnegative. Note that all the vital functions (g,µand α) depend on size x(the structuring internal variable), on time tand on the value of the resource at time t, which can reflect the influence of the environmental changes on the vital functions. Function fon the right-hand side of (4) depends on the value of the resource at time t, on the total amount of individuals in the population by means of the weighted functional I(t) (which represents the way of weighting the size distribution density in order to model the different influence of individuals of different sizes on such dynamics) and on time t. The paper is structured as follows. In the next section we introduce some background on theoretical and numerical results. In section III, we describe the numerical method, which is completely analyzed in section IV. Finally, numerical results that confirm the expected order of convergence and show the biological example dynamics are included in section V. II. Preliminary results Theoretical analysis (1)-(5) which includes existence, uniqueness and long-time behaviour, is highly difficult. A theoretical study of this model appeared first time in a H.Thieme’s [7] presentation and authors in [5] also pointed that, for an analysis of such models, we had to perform as it was made in [8]. However, only a simpler model, in which a positive growth was employed, was analyzed in [9]. The numerical solution to the model, due to its obvious mathematical complexity, entails a serious challenge. In [10], a simple modification of the scheme succesfully employed in [11], for the solution of general size-structured population models, was considered. The Daphnia magna test in section V was also introduced. Furthermore, the theoretical steady state of the model for this test was provided. This scheme was shown not suitable for a long-time integration with this biological test. In this work, the numerical method described in section III was proposed. Some simulations was performed to show values of the parameters which led us to an asymptotically stable equilibrium and to an asymptotically stable periodic situation, in both cases solutions were bounded. The purpose of the current work is to validate such numerical method by means of an analysis of its convergence. On the other hand, a modification of this numerical procedure was introduced in [12] in order to approximate singular asymptotic states for these kinds of models. In this work, we showed stable and unstable steady state solutions. The study of these singular states is not the aim of present work. Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 2 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... In our convergence analysis, we shall use the general discretization framework introduced by L´ opez-Marcos et al. [13]. We present it rather tersely, the reader is referred to this paper and the references therein for a more critical and detailed treatment. Thus, once we introduced a fixed given problem concerning a differential equation, we denote ua theoretical solution to such a problem. We denote ˜ Uha numerical approximation to u. The subscript hreflects the dependency on a parameter h(mesh size) which takes values in a set Hof positive numbers with inf H=0. The approximation is reached by solving the discretized problem Φh(˜ Uh)=0(8) (this family of discrete problems with his referred as discretization), where, for each h∈H, the mapping Φhis fixed with domain Dh⊂ Ahand taking values in Bh. Here, Ahand Bhare vector spaces with the same finite dimension. We further assume that, for each h∈H, we have chosen a norm in both spaces and an element ˜ uh∈Dhwhich is a suitable discrete representation of uin Ah. Furthermore, we introduce the global and local discretization error, ˜ ehand lh, respectively, and the consistency, stability and convergence properties of a discretization. The following result is crucial and uses a deep topological lemma due to Stetter [14], Theorem 1. Assume that (8) is consistent and stable with thresholds Rh. If Φhis continuous in B(˜ uh,Rh)and klhkBh=o(Rh)as h →0, then: i) For h sufficiently small, the discrete equations (8) possess a unique solution in B(˜ uh,Rh). ii) As h →0, the solutions converge and k˜ ehkAh=O(klhkBh). III. The numerical method The numerical method we employ to approximate the solution to (1)-(5) is based on the discretization of a representation of the solution along the characteristic curves [10]. First of all we rewrite the partial differential equation (1) in a more suitable form for its numerical treatment. So we define µ∗(x,z,t)=µ(x,z,t)+gx(x,z,t). Thus equation (1) has the form ut+g(x,S(t),t)ux=−µ∗(x,S(t),t)u,(9) 0<x<xM(t),t>0. We denote by x(t;t∗,x∗) the characteristic curve of equation (9) that takes the value x∗at time t∗. Such a characteristic curve is the solution to the initial value problem d dt x(t;t∗,x∗)=g(x(t;t∗,x∗),S(t),t),t≥t∗, (10) x(t∗;t∗,x∗)=x∗. Now we consider the function that represents the solution to (9) along the characteristic curves u(x(t;t∗,x∗),t), t≥t∗, which satisfies the initial value problem d dtu(x(t;t∗,x∗),t)= −µ∗(x(t;t∗,x∗),S(t),t)u(x(t;t∗,x∗),t),t≥t∗, u(x(t∗;t∗,x∗),t∗)=u(x∗,t∗), and, therefore, it can be represented in the following integral form u(x(t;t∗,x∗),t)=(11) u(x∗,t∗)exp(− Zt t∗ µ∗ (x(τ;t∗,x∗),S(τ), τ)dτ),t≥t∗. 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, as a parameter related to the size variable which describes the number of points in the uniform initial grid. The diameter of 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 as an approximation to the density at initial time (t0), the grid restriction of the initial condition in (3), U0 j=u0(X0 j), 0 ≤j≤J. Also, we use S0 in (4) as the initial value of the resource. Then, the numerical method provides, at each discrete time level, a mesh grid on the size interval in Xn, the approximation to the density on such a Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 3 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... mesh grid in Unand the approximation of the value of the resource Sn, from the approximations we computed at the previous time level, by using discretizations of equations (10), (11), (2), (4), (5). Thus, for n=0,1,2, . . ., the numerical solution at time tn+1=tn+k, is obtained from the known values of the numerical solution at time tnas follows, Xn+1 0=0,(12) Xn+1 j+1=Xn j+k g(Xn+1/2 j+1,Sn+1/2,tn+1/2),(13) 0≤j≤J,(14) Sn+1=Sn+(15) k f (Sn+1/2 ,Q(Xn+1/2 ,γn+1/2.Un+1/2),tn+1/2),(16) Un+1 j+1=Un jexp n−kµ∗(Xn+1/2 j+1,Sn+1/2 ,tn+1/2)o,(17) 0≤j≤J,(18) Un+1 0=Q(Xn+1,αn+1.Un+1) g(Xn+1 0,Sn+1,tn+1),(19) where we have to compute approximations at time level tn+1/2=tn+k/2, Xn+1/2 0=0, Xn+1/2 j+1=Xn j+k 2g(Xn j,Sn,tn),0≤j≤J, Sn+1/2=Sn+k 2f(Sn,Q(Xn,γn.Un),tn), Un+1/2 j+1=Un jexp(−k 2µ∗(Xn j,Sn,tn)),0≤j≤J, Un+1/2 0=Q(Xn+1/2,αn+1/2.Un+1/2) g(Xn+1/2 0,Sn+1/2,tn+1/2). Note that the general step of the method increases the number of grid points and also the dimension of the vector with the numerical densities: at time tn, we have J+1 grid nodes in Xnand the (J+1)- dimensional vector Un, and at time tn+1we obtain J+2 grid nodes in Xn+1and the (J+2)-dimensional vector Un+1. In order to maintain the number of grid points suitable to perform the next step, we eliminate at time tn+1the first grid node Xn+1 l which satisfies |Xn+1 l+1−Xn+1 l−1|=min 1≤j≤J+1|Xn+1 j+1−Xn+1 j−1|.(20) We reproduce the same reduction in the corresponding vector Un+1. In the description of the method, we use the following notation; vectors αpand γpcontain the evaluations of the functions αand γin (2) and (5), respectively, at the grid points in Xp, at the resource value Spand at time tp. Products γp.Up and αp.Upmust be considered componentwise. In order to approximate integrals over the interval [0,xM(tp)], we use the composite trapezoidal quadrature rule based on the grid points Xp= [Xp 0,Xp 1,...,Xp J], that is Q(Xp,Vp)= J X j=1 Xp j−Xp j−1 2Vp j−1+Vp j.(21) Note that the method is implicit: all the expressions provide explicit equations for the numerical values at the highest time level, except those which involve the numerical density Up 0at the first grid point, but it is easy to implement the method in an explicit form. IV. Convergence Analysis Below, we will analyze numerical methods based on the integration along characteristics that use a general quadrature rule with suitable properties to approximate the integral terms. The proofs of every result are heavily laborious and they would be included in a more technical work. We assume the following regularity conditions on the data functions and the solution to the problem (1)-(5): (H1) u∈ C2([0,xM(t)] ×[0,T]), u(x,t)≥0, x∈ [0,xM(t)], t≥0. (H2) S∈ C2([0,T]), S(t)≥0, t≥0. (H3) γ∈ C2([0,xM(t)]×D×[0,T]), where Dis a compact neighbourhood of {S(t),0≤t≤T}. (H4) µ∈ C2([0,xM(t)] ×D×[0,T]), is nonnegative and Dis a compact neighbourhood of {S(t),0≤t≤T}. (H5) α∈ C2([0,xM(t)] ×D×[0,T]), is nonnegative and Dis a compact neighbourhood of {S(t),0≤t≤T}. (H6) f∈ C2(D×DI×[0,T]), is nonnegative, Dis a compact neighbourhood of Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 4 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... {S(t),0≤t≤T}, and DIis a compact neighbourhood of (ZxM(t) 0 γ(x,S(t)) u(x,t)dx,0≤t≤T). (H7) g∈ C3([0,xM(t)]×D×[0,T]), where Dis a compact neighbourhood of {S(t),0≤t≤T} and g(0,s,t)≥C>0, t≥0, s∈R. In addition, the characteristic curves x(t;t∗,x∗) defined in (7) are continuous and differentiable with respect to the initial values (t∗,x∗)∈ [0,T]×[0,xM(t)]. The above hypotheses may be based on three possible reasons. First, biological assumptions such as the nonnegativity of some of the vital functions or, in (H7), to reflect that any individual in the studied population could shrink [2]. Second, the mathematical requirements to obtain the existence and uniqueness of solutions for the problem (1)- (5) [2]. Finally, the regularity properties needed in the numerical analysis to derive optimal rates of convergence [11]. We also assume that the spatial discretization parameter, h, takes values in the set H={h>0 : h=xM/J,J∈N}. Now, we suppose that the time step, k, satisfies k=r h, where ris an arbitrary and positive constant, fixed throughout the analysis. In addition, we set N=[T/k]. For each h∈H, we define the spaces Ah= N Y n=0RJ+n×RJ+n+1×RN+1, Bh=RJ×RJ+1×R×RN× N Y n=1RJ+n×RJ+n×RN. Both spaces have the same dimension. In order to measure the size of the errors, we define kηk∞=max1≤j≤p|ηj|,η∈Rp,kVnk1= PJ+n j=0h|Vn j|,Vn∈RJ+n+1. Thus, we endow the spaces Ahand Bhwith the following norms. If y0,V0,...,yN,VN,a∈ Ah, then ky0,V0,...,yN,VN,akAh= max max 0≤n≤Nkynk∞,max 0≤n≤NkVnk∞,kak∞. On the other hand, if Y0,Z0,A0,Z0,Y1,Z1,...,YN,ZN,A∈ Bh, kY0,Z0,A0,Z0,Y1,Z1,...,YN,ZN,AkBh =kY0k∞+kZ0k∞+|A0|+kZ0k∞+ N X n=1 kkZnk∞ + N X n=1 kkYnk∞+ N X n=1 k|An|. Now, for each h∈H, we define xh= (x0,x1,x2,...,xN), xn=(xn 1,...,xn J+n)∈RJ+n, x0 j=j h, 1 ≤j≤J, xn j=x(tn;tn−1,xn−1 j−1),1≤j≤J+n,1≤n≤N. (22) Also, xn+1 2=(xn+1 2 1,...,xn+1 2 J+n+1)∈RJ+n+1, xn+1 2 j=x(tn+1 2;tn,xn j−1),1≤j≤J+n+1,(23) 0≤n≤N−1. We denote xn 0=xn+1 2 0=0, n≥0. In addition, if urepresents the theoretical solution to (1)-(5) we define uh=(u0,u1,u2,...,uN), un= (un 0,un 1,...,un J+n)∈RJ+n+1, un j=u(xn j,tn),0≤j≤J+n,0≤n≤N,(24) and un+1 2=(un+1 2 0,un+1 2 1,...,un+1 2 J+n+1)∈RJ+n+2, un+1 2 j=u(xn+1 2 j,tn),0≤j≤J+n+1,0≤n≤N−1. (25) Finally, if Sis the theoretical solution to (4) then we define sh=(s0,s1,s2,...,sN), sn=S(tn),0≤n≤N,(26) and sn+1 2=S(tn+1 2),0≤n≤N−1.(27) Therefore ˜ uh=(x0,u0,x1,u1,...,xN,uN,sh)∈ Ah. Next, we introduce the discretization operator. Let Rbe a positive constant and we denote by BAh(˜ uh,R hp)⊂ Ahthe open ball with center ˜ uh and radius R hp, 1 <p<2, Φh:BAh(˜ uh,R hp)→ Bh, Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 5 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... Φhy0,V0,...,yN,VN,a =Y0,P0,A0,P0,Y1,P1,...,YN,PN,A,(28) defined by the following equations: Y0=y0−X0∈RJ,(29) P0=V0−U0∈RJ+1,(30) A0=a0−S0∈R.(31) Vectors X0,U0and value S0represent approximations at t=0, respectively, to the initial grid nodes, to the theoretical solution at these points and to the initial resource. Also, Pn+1 0=Vn+1 0− Qyn+1,αn+1·Vn+1 g0,an+1,tn+1,(32) Yn+1 j+1=1 knyn+1 j+1 −yn j−k g(yn+1 2,∗ j+1,an+1 2,∗,tn+1 2),(33) Pn+1 j+1=1 knVn+1 j+1 (34) −Vn jexp −kµ∗yn+1 2,∗ j+1,an+1 2,∗,tn+1 2,(35) 0≤j≤J+n−1, An+1=1 knan+1−an −k f an+1 2,∗,Q(yn+1 2,∗,γn+1 2,∗·Vn+1 2,∗),tn+1 2o,(36) 0≤n≤N−1. Where, with the notation introduced in Section III, yn+1 2,∗ j+1=yn j+k 2g(yn j,an,tn),(37) Vn+1 2,∗ j+1=Vn jexp −k 2µ∗yn j,an,tn!,(38) 0≤j≤J+n−1, Vn+1 2,∗ 0=Q(yn+1 2,∗,αn+1 2,∗·Vn+1 2,∗) g(0,an+1 2,∗,tn+1 2) ,(39) an+1 2,∗=an+k 2f(an,Q(yn,γn·Vn),tn),(40) 0≤n≤N−1. We denote by Q(X,V)= M X l=0 ql(X)Vlthe general quadrature rule employed in (32)-(40). Note that, Φhtakes into account all the possible nodes and their corresponding solution values at each time level, and it employs quadrature rules possibly based on a subgrid. If ˜ Uh=(X0,U0,X1,U1,...,XN,UN,S)∈ BAh(˜ uh,R hp), satisfies Φh(˜ Uh)=0∈ Bh,(41) the nodes and the corresponding values of the solution at such nodes of ˜ Uhare a numerical solution to the scheme defined by (13)-(19) when the composite trapezoidal quadrature rule is given. On the other hand, the numerical solution of the scheme defined by (13)-(19) satisfies (41). Henceforth, Cwill denote a positive constant, independent of h,k(k=r h), j(0 ≤j≤J+n) and n(0 ≤n≤N); Cmay have different values in different places. Now, we assume that the quadrature rules satisfy the following properties: (P1) |I(tn)− Q (xn,γn·un)|≤C h2, when h→ 0, 0 ≤n≤N. (P2) I(tn−1 2)− Q xn−1 2,γn−1 2·un−1 2≤C h2, when h→0, 1 ≤n≤N. (P3) ZxM(tn) 0 α(x,S(tn),tn)u(x,tn)dx − Q (xn,αn·un) ≤C h2, when h→0, 0 ≤n≤N. (P4) ZxM(tn−1 2) 0 α(x,S(tn−1 2),tn−1 2)u(x,tn−1 2)dx − Q xn−1 2,αn−1 2·un−1 2≤C h2, when h→0, 1≤n≤N. (P5) |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, 0≤n≤N. (P6) |qj(xn−1 2)| ≤ 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, 1≤n≤N. (P7) Let Rand pbe positive constants with 1<p<2. The quadrature weights qjare Lipschitz continuous functions on B∞(xn,R hp), Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 6 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... 0≤j≤J+n, 1 ≤n≤Nand on B∞(xn−1 2,R hp), 0 ≤j≤J+n, 1 ≤n≤N. (P8) 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), then  J+n X i=0 (qi(yn)−qi(zn))γ(zn i,an,tn)Vn i ≤ Ckyn−znk∞, when h→0, 0 ≤n≤N. (P9) 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), then  J+n X i=0 (qi(yn)−qi(zn))αzn i,an,tnVn i ≤ Ckyn−znk∞, when h→0, 0 ≤n≤N. (P10) Let Rand pbe positive constants with 1 < p<2. If yn−1 2,zn−1 2∈B∞(xn−1 2,R hp), Vn−1 2∈ B∞(un−1 2,R hp), and an−1 2∈B∞(sn−1 2,R hp), then  J+n+1 X i=0 qi(yn−1 2)−qi(zn−1 2) γzn−1 2 i,an−1 2,tn−1 2Vn−1 2 i ≤Ckyn−1 2−zn−1 2k∞, when h→0, 1 ≤n≤N. (P11) Let Rand pbe positive constants with 1 < p<2. If yn−1 2,zn−1 2∈B∞(xn−1 2,R hp), Vn−1 2∈ B∞(un−1 2,R hp), and an−1 2∈B∞(sn−1 2,R hp), then  J+n+1 X i=0 qi(yn−1 2)−qi(zn−1 2) αzn−1 2 i,an−1 2,tn−1 2Vn−1 2 i ≤Ckyn−1 2−zn−1 2k∞, when h→0, 1 ≤n≤N. These are the enough properties the quadrature rules have to satisfy to carry out our convergence analysis. The following result establishes that the composite trapezoidal rule used in our experiments satisfies them. Theorem 2. Assume that the hypotheses (H1)- (H7) hold. If the quadrature rules are the composite trapezoidal quadrature on subgrids xn jn lM(n) l=0 , 0≤n≤N with the property (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, with xn jn lM(n)−1 l=1 contained in xn,0≤n≤N. Then, properties (P1)-(P11) hold. Now we introduce the following result over the numerical values at the half-level time. Proposition 1. Assume that the hypotheses (H1)- (H7) hold and that the considered quadrature rules satisfy properties (P1)-(P11). Let be yn∈B∞(xn,R hp),Vn∈B∞(un,R hp)and an∈B∞(sn,R hp). Then, as h →0,yn+1 2,∗∈ B∞(xn+1 2,R0hp), an+1 2,∗∈B∞(sn+1 2,R0hp), and Vn+1 2,∗∈B∞(un+1 2,R0hp)where xn+1 2, sn+1 2, and un+1 2are defined by (23), (27) and (25), respectively. Now, we define the local discretization error as lh=Φh(˜ uh)∈ Bh, and we say that the discretization (28) is consistent if, as h→0, lim kΦh(˜ uh)kBh=lim klhkBh=0. The following theorem establishes the consistency of the numerical scheme defined by equations (29)-(36). Theorem 3. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P11). Then, as h →0, the local discretization error satisfies, kΦh(˜ uh)kBh=ku0−U0k∞+kx0−X0k∞ +|s0−S0|+O(h2+k2).(42) Another notion that plays an important role in the analysis of the numerical method is the stability with h-dependent thresholds. For h∈H, let Rhbe a real number ( the stability threshold) with Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 7 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... 0<Rh<∞: we say that the discretization (28) is stable for ˜ uhrestricted to the thresholds Rh, if there exist two positive constants h0and S ( the stability constant) such that, for any h∈Hwith h≤h0, the open ball BAh(˜ uh,Rh) is contained in the domain of Φh, and, for all ˜ Vh,˜ Whin that ball, k˜ Vh−˜ Whk ≤ SkΦh(˜ Vh)−Φh(˜ Wh)k. Below, we introduce the theorem that establishes the stability of the discretization defined by equations (29)-(36). Theorem 4. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P11). Then, the discretization is stable for ˜ uhwith Rh=R hp,1<p<2. Finally, we define the global discretization error as ˜ eh=˜ uh−˜ Uh∈ Ah, We say that the discretization (28) is convergent if there exists h0>0 such that, for each h∈H with h≤h0, (41) has a solution ˜ Uhfor which, as h→0, lim k˜ uh−˜ UhkAh=lim k˜ ehkAh=0. We propose the following theorem which establishes the convergence of the numerical method defined by equations (29)-(36). Theorem 5. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P11). Then, for h sufficiently small, the numerical method defined by equations (29)-(36) has a unique solution Uh∈ B(uh,Rh)and kUh−uhkAh≤Ckx0−X0k∞+ku0−U0k∞ +|s0−S0|+O(h2+k2).(43) The proof of Theorem 5 is immediately derived by means of the consistency (Theorem 3), the stability (Theorem 4) and Theorem 1. Next, we can establish an error bound for the the numerical and the theoretical solution at the numerical values of the grid nodes. Theorem 6. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P11). For h sufficiently small, let be u∗ h=(u0 ∗,u1 ∗,u2 ∗,...,uN ∗)∈ N Y n=0 RJ+n+1, defined by un ∗=u(Xn 0,tn),u(Xn 1,tn),...,u(Xn J+n,tn)∈RJ+n+1, 0≤n≤N, where Xn j,0≤j≤J+n, 0≤n≤N, are the grid nodes given by the scheme (29)-(36). Then, kUn−un ∗k∞≤Ckx0−X0k∞+ku0−U0k∞ +|s0−S0|+O(h2+k2).(44) This theorem follows immediately from Theorem 5. Specifically, if X0=x0,U0=u0and S0=s0, the proposed numerical scheme is secondorder accurate. At this moment, we have obtained convergence of the numerical method (29)-(36) which does not employ selection at each time level. Also, we have proven the convergence of numerical methods which employ a selection criterion, whenever the positions, which are determined by the criterion we have chosen, lead us to subgrids which satisfy property (SR). For the criterion presented in this paper, this property may be shown in two stages. First, as proved in [11], it leads us to subgrids with such a property when we applied it over nodes which are in a neighbourhood of the theoretical ones with radius R hp. In a second stage, it is proven that the nodes, which in fact the numerical method computes, are in such neighbourhoods. In order to do this, it is enough to realize that such nodes could be seen, up to each time level, as the solutions obtained by a discrete operator which has the form of that defined in (28). V. Numerical results We have carried out different numerical experiments with the scheme defined in Section III. We have considered a theoretical test problem with meaningful nonlinearities (both from a mathematical and biological point of view). The numerical integration for the numerical experiment was Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 8 of 11 O Angulo et al., Numerical Analysis of a Size-Structured Population Model... carried out on the time interval [0,10]. The size interval was taken as [0,1]. The size-specific growth, fertility and mortality moduli are chosen as g(x,z,t)=λ 2 1+z zz 1+z2−x2+xr 1+z29 30 −z k, α(x,z,t)=3 2λ 1+ z C(29 30 −z k)!−29λ 30r 1+2 z C(29 30 −z k)!−29λ 30r ,µ(x,z,t)= λ1+z zz 1+z+2x−3r 1+z29 30 −z k. The weight function is taken as γ(x,z,t)=x2and, finally, f(z,i,t)=rz 1−z k−rzi (1+z)5 z5(1+4e−λt). With this choice of data functions, the problem (1)-(5) has the following solution u(x,t)= S(t) 1+S(t)−x!2 −e−λt       S(t) 1+S(t)!2 −x2      , S(t)=29 30 Ce29rt/30 1+Ce29rt/30/k, r=0.1, C=24, k=5, λ=0.3. Since we know the exact solution to the problem, we can show numerically that our method is second-order accurate by means of an error Table. In Table I, each entry in columns two to five represents, at the upper value, the global error 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|) and, at the lower number, the experimental order s of the method as computed from s=log (e2h,2k/eh,k) log 2 . Each column and each row of the Table correspond to different values of the spatial and time discretization parameter, respectively. The results in the Table clearly confirm the expected secondorder convergence. On the other hand, the numerical integration of the model with an efficient method allows us to consider a more realistic test problem. The property of convergence in finite time interval is important to carry out experiments in which the long-time behaviour of the population is investigated. Therefore, the numerical method has been employed to describe the dynamics of a population k\h1.25e-2 6.25e-2 3.13e-2 1.56e-2 1.25e-2 1.67e-4 1.27e-4 1.22e-4 1.21e-4 6.25e-2 2.08e-4 4.18e-5 3.15e-5 3.03e-5 2.00 2.01 2.01 3.13e-2 2.02e-4 5.26e-5 1.05e-5 7.85e-6 1.98 2.00 2.00 1.56e-2 2.08e-4 5.05e-5 1.32e-5 2.62e-6 2.00 1.99 2.00 TABLE I Error and experimental order of convergence. of ectothermic invertebrates. This is the case of the water flea, Daphnia magna. In this particular case, the functions data are given by g(x,S,t)= g(S 1+S−x), µ(x,S,t)=µ,α(x,S,t)=αS 1+Sx2, f(S,i,t)=rS 1−S K−iS 1+S,γ(x,S,t)=x2, and the values of the parameters are given by g=1, µ=0.1, α=0.75, r=3, K=8.3 and xM=1 [6]. This set led us to unbounded solutions [12]. Nevertheless, for the parameter value g=0.0075, the solution is bounded. As it was pointed in [10], as the value of the parameter Kincreases, the equilibrium state becomes unstable. We have performed a numerical experiment with the discretization parameters k=0.0625, J=4000 and the interesting value K=9.64. In this case, we observe the unstability of the equilibrium and the solution evolving towards a cycled situation (Figure 1). Taking into account that the numerical solution is attracted to a limit cycle, considering a sufficiently large time, the numerical solution obtained after this long time integration lies practically on such a cycle. In this way, the numerical method provides an approximation to the limit cycle. In Figure 2, the representation of such a cycle in the tridimensional space defined by the total population, the maximum individual size and the dynamical resource, is drawn. From the numerical results obtained for the total population, the maximum size and the dynamical resource, we can obtain an in depth analysis of these quantities throughout a period of the limit cycle. For example, we have estimated the period of the solution by interpolation and it is about 64.6824. Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 Page 9 of 11