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

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 al-lows 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 u t + (g(x, S (t), t) u) x = −µ(x, S (t), t) u, (1) 0 < x < x M (t), t > 0, a nonlinear and nonlocal boundary condition which reflects the reproduction process g(0, S (t), t) u(0, t) = x M (t) 0 α(x, S (t), t) u(x, t) dx, (2) and an initial size-distribution for the population u(x, 0) = u 0 (x), 0 ≤ x ≤ x M (0).
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, S (t) = f (S (t), I(t), t), t > 0, (4) S (0) = S 0 , 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 x M (t) 0 γ(x, S (t), t) u(x, t) dx, t ≥ 0. ( 5) Also, we consider that the maximum size x M (t) that an individual could have at time t, changes with time and its dynamics is described by x M (t) = g(x M (t), S (t), t), t > 0, (6) x M (0) = x M .
The independent variables x and t represent size and time, respectively.The dependent variable u(x, t) is the size-specific density of individuals with size x at time t.In general, the size of any individual varies according to the following ordinary differential equation 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 t and on the value of the resource at time t, which can reflect the influence of the environmental changes on the vital functions.Function f on 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.
In our convergence analysis, we shall use the general discretization framework introduced by López-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 u a theoretical solution to such a problem.We denote Ũh a numerical approximation to u.The subscript h reflects the dependency on a parameter h (mesh size) which takes values in a set H of positive numbers with inf H = 0.The approximation is reached by solving the discretized problem (this family of discrete problems with h is referred as discretization), where, for each h ∈ H, the mapping Φ h is fixed with domain D h ⊂ A h and taking values in B h .Here, A h and B h are 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 ũh ∈ D h which is a suitable discrete representation of u in A h .Furthermore, we introduce the global and local discretization error, ẽh and l h , 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 R h .If Φ h is continuous in B( ũh , R h ) and l h B h = o(R h ) as h → 0, then: i) For h sufficiently small, the discrete equations (8) possess a unique solution in B( ũh , R h ).ii) As h → 0, the solutions converge and

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) + g x (x, z, t).

Thus equation (1) has the form
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 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 and, therefore, it can be represented in the following integral form Given a constant step k > 0, we introduce the discrete time levels t n = n k, n = 0, 1, 2, . ... We also take J a 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 = x M /J, and the initial grid nodes are X 0 j = j h, 0 ≤ j ≤ J.In order to start the integration, we consider as an approximation to the density at initial time (t 0 ), the grid restriction of the initial condition in (3), U 0 j = u 0 (X 0 j ), 0 ≤ j ≤ J. Also, we use S 0 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 X n , the approximation to the density on such a Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 mesh grid in U n and the approximation of the value of the resource S n , 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 t n+1 = t n + k, is obtained from the known values of the numerical solution at time t n as follows, ), ( 16) where we have to compute approximations at time level ) .
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 t n , we have J + 1 grid nodes in X n and the (J + 1)dimensional vector U n , and at time t n+1 we obtain J+2 grid nodes in X n+1 and the (J+2)-dimensional vector U n+1 .In order to maintain the number of grid points suitable to perform the next step, we eliminate at time t n+1 the first grid node X n+1 l which satisfies We reproduce the same reduction in the corresponding vector U n+1 .
In the description of the method, we use the following notation; vectors α p and γ p contain the evaluations of the functions α and γ in (2) and (5), respectively, at the grid points in X p , at the resource value S p and at time t p .Products γ p .U p and α p .U p must be considered componentwise.In order to approximate integrals over the interval [0, x M (t p )], we use the composite trapezoidal quadrature rule based on the grid points 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 U p 0 at 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): In addition, the characteristic curves x(t; t * , x * ) defined in (7) are continuous and differentiable with respect to the initial values (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 = x M /J, J ∈ N}.Now, we suppose that the time step, k, satisfies k = r h, where r is an arbitrary and positive constant, fixed throughout the analysis.In addition, we set N = [T/k].For each h ∈ H, we define the spaces Both spaces have the same dimension.
In order to measure the size of the errors, we define . Thus, we endow the spaces A h and B h with the following norms.If y 0 , V 0 , . . ., y N , V N , a ∈ A h , then = 0, n ≥ 0. In addition, if u represents the theoretical solution to (1)-( 5) we define and u n+ 1 2 = (u (25) Finally, if S is the theoretical solution to (4) then we define s h = (s 0 , s 1 , s 2 , . . ., s N ), and Next, we introduce the discretization operator.Let R be a positive constant and we denote by B A h ( ũh , R h p ) ⊂ A h the open ball with center ũh and radius R h p , 1 < p < 2, Φ h y 0 , V 0 , . . ., y N , V N , a = Y 0 , P 0 , A 0 , P 0 , Y 1 , P 1 , . . ., Y N , P N , A , (28) defined by the following equations: Vectors X 0 , U 0 and value S 0 represent approximations at t = 0, respectively, to the initial grid nodes, to the theoretical solution at these points and to the initial resource.Also, Where, with the notation introduced in Section III, q l (X) V l the general quadrature rule employed in (32)-(40).Note that, Φ h takes 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.
the nodes and the corresponding values of the solution at such nodes of Ũh are 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, C will denote a positive constant, independent of h, k (k = r h), j (0 ≤ j ≤ J + n) and n (0 ≤ n ≤ N); C may have different values in different places.Now, we assume that the quadrature rules satisfy the following properties: where q is a positive constant independent of h, k, j (0 ≤ j ≤ J +n) and n (0 where q is 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 R and p be positive constants with 1 < p < 2. The quadrature weights q j are Lipschitz continuous functions on (P8) Let R and p be positive constants with 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.
Then, properties (P1)-(P11) hold.Now we introduce the following result over the numerical values at the half-level time.

Now, we define the local discretization error as
and we say that the discretization (28 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, 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 R h be a real number ( the stability threshold) with Biomath 3 (2014), 1403241, http://dx.doi.org/10.11145/j.biomath.2014.03.241 0 < R h < ∞: we say that the discretization (28) is stable for ũh restricted to the thresholds R h , if there exist two positive constants h 0 and S ( the stability constant) such that, for any h ∈ H with h ≤ h 0 , the open ball B A h ( ũh , R h ) is contained in the domain of Φ h , and, for all Ṽh , Wh in that ball, 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 Finally, we define the global discretization error as ẽh We say that the discretization (28) is convergent if there exists h 0 > 0 such that, for each h ∈ H with h ≤ h 0 , (41) has a solution Ũh for which, as h → 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 U h ∈ B(u h , R h ) and 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 defined by where X n j , 0 ≤ j ≤ J + n, 0 ≤ n ≤ N, are the grid nodes given by the scheme (29)-(36).Then, This theorem follows immediately from Theorem 5. Specifically, if X 0 = x 0 , U 0 = u 0 and S 0 = s 0 , 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 h p .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 The weight function is taken as γ(x, z, t) = x 2 and, finally, . With this choice of data functions, the problem (1)-( 5) has the following solution 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 and, at the lower number, the experimental order s of the method as computed from s = log (e 2h,2k /e h,k ) log 2 . 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 of ectothermic invertebrates.This is the case of the water flea, Daphnia magna.In this particular case, the functions data are given by g

Each column and each row of the
, and the values of the parameters are given by g = 1, µ = 0.1, α = 0.75, r = 3, K = 8.3 and x M = 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 K increases, 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.Limit cycle appearing in the case of a bounded unstable steady state.

VI. Conclusions
We have analyzed a second-order numerical method for a problem that describes a population with a possible shrinking size and with a dependency on the environment managed by the evolution of a vital resource.The second-order convergence has been theoretically proven by means of an argument of consistency and stability of the scheme.We have reported numerical experiments which demonstrate the predicted accuracy of the scheme.
This knowledge leads us to employ this second-order numerical method that, experimentally, shows good stability properties for the study of the behaviour of a real population.We have consider the long time numerical approximation of a particular model which describes the dynamics of a Daphnia magna population.Moreover, when the steady state is unstable we have used it to analyzed the dynamics, appearing a limit cycle, and a good approximation to the bifurcation process has been obtained.

Theorem 2 .
Assume that the hypotheses (H1)-(H7) hold.If the quadrature rules are the composite trapezoidal quadrature on subgrids x n ≤ N with the property (SR) There exists a positive constant C such that, for h sufficiently small, x n j

Fig. 1 .
Fig. 1.Evolution of the numerical solution.Case of a bounded unstable steady state.
Table correspond to different values of the spatial and time discretization parameter, respectively.The results in the Table clearly confirm the expected secondorder convergence.