The Convergence Analysis of a Numerical Method for a Structured Consumer-Resource Model with Delay in the Resource Evolution Rate

In this paper, we go through the development of a new numerical method to obtain the solution to a size-structured population model that describes the evolution of a consumer feeding on a dynamical resource that reacts to the environment with a lag-time response. The problem involves the coupling of the partial differential equation that represents the population evolution and an ordinary differential equation with a constant delay that describes the evolution of the resource. The numerical treatment of this problem has not been considered before when a delay is included in the resource evolution rate. We analyzed the numerical scheme and proved a second-order rate of convergence by assuming enough regularity of the solution. We numerically confirmed the theoretical results with an academic test problem.


Introduction
We consider a model that describes the evolution of a size-structured consumer in an environment inhabited by a single unstructured resource. We assume that the resource responds to the environment with a constant time delay. The model is composed of two nonlinear coupled problems. On the one hand, a size-structured population model governed by a first-order hyperbolic equation with a nonlocal and nonlinear boundary condition, which represents the evolution over time of the consumer: Variables x and t represent the size and the time, respectively. Size is the variable which structures the individuals in the consumer population, x 0 denotes the newborn individual's size (the size at birth) that is assumed to be constant and positive and x M (t) represents the maximum size of individuals at time t. Then, x M (0) = x 0 M is the initial maximum size. The dependent variables u and s are the density of individuals of size x and the amount of the resource available at time t, respectively. The so-called vital rates, the mortality, the fertility and the growth rates, are given by µ, α and g, respectively. The mortality and the fertility are nonnegative functions and the growth rate has no-sign restriction, although we assume that g(x 0 , ·, ·) > 0, which means that each individual increases in size at birth. Vital rates depend on the structuring variable, the time and the amount of the resource available, to take into account the influences of these factors on the dynamics of the population. The size of the individuals changes according to the differential equation x (t) = g(x(t), s(t), t), t > 0.
As we allow a negative growth rate g, we are considering the case in which an individual in the population could shrink in size under a food shortage environmental condition. In particular, the maximum size, x M (t), is not fixed in the model and evolves following the corresponding characteristic curve of (1) x M (t) = g(x M (t), s(t), t), t > 0, On the other hand, the unstructured resource s(t), t ≥ 0, which provides feeding for the individuals of the population, evolves with time according to the following initial value problem for a delay ordinary differential equation: The rate of change in the evolution of the resource is given by a function f that includes dependence on time and on the consumer population through the nonlocal term where γ is a function that represents the individual rate of consumption for individuals of a determined size. It also depends on the amount of the resource available at times t and t − τ, s(t) and s(t − τ) respectively. The memory effect, s(t − τ), can be seen as a deferred influence on the environment of the resource affordability. A common biological situation for this phenomenon occurs when the resource population is close to the carrying capacity of the environment or near extinction, which can make the population react with a certain delay [1]. Another situation in which this deferred influence occurs is when the resource is formed only by the adult individuals of the population, and the delay is due to the maturation time [2]. The solution of the model (1), (3) is determined once we know the initial conditions u(x, 0) = u 0 (x), x 0 ≤ x ≤ x 0 M , and s(t) = s 0 (t), t ∈ [−τ, 0]. Although a large number of papers on consumer models have been considered in past years, not so many address the case in which one of the populations is structured. Consumer-resource models have been studied from the very initial work of Kooijman and Metz [2]. These authors presented a mathematical model for the development of an ectothermic population (Daphnia magna, water flea) in which the amount of food they are supplied with represents a regulatory mechanism for the population density. The model was length-structured and also included a stage of maturity: a differentiation among juveniles and adults. This work is considered as the origin of the modern dynamic energy budget theory. The well-posedness of the model equations and the continuously dependence on model data was studied by Thieme [3]. The dynamical properties of the model were explored numerically in [4], and later, analytically in [5][6][7]. However the description of the evolution of the resource as a delay differential equation within a physiological structured consumer-resource model has not been developed theoretically yet.
The theoretical treatment of this model is not easy; therefore, its numerical analysis is a valid tool and sometimes the only one affordable. The integration of the model without delay has been developed by means of the Excalator Boxcar Train (EBT) [8] and a characteristics method [9]. This last work included the convergence proof of the method while the convergence of the popular EBT method was recently considered [10,11]. However, the numerical treatment of the coupled model with a resource evolving according to a delay differential equation remained unexplored until our study, so our goal was to provide a numerical method to perform the integration of such a model and its convergence analysis.
The remainder of the paper is organized as follows. Section 2 is devoted to the introduction of the numerical method to integrate problem (1), (3). We employ the technique of integration along the characteristic curves in order to obtain the new numerical scheme. In Section 3, we carry out the convergence analysis of the scheme. It is based on a theoretical framework that involves the properties of consistency and stability of the numerical method. We also pay attention to the properties required by the numerical quadrature rule used in the integration. Finally, in Section 4, we present a test done in order to numerically confirm the theoretical order of convergence.

Numerical Method
We begin with the description of the numerical scheme. We must employ the integration along the characteristic curves; therefore, we should remind the reader of the definition of x(t; t * , x * ), that represents the characteristic curve that at time t * starts at x * . Now, we consider w(t; t * , where µ * (x, s(t), t) = µ(x, s(t), t) + g x (x, s(t), t), whose solution is given by the following formula For the numerical method, we discretize this expression, together with the boundary condition and the initial data in (1), and the initial value problem in (3).
The integration is carried out on a finite time interval [0, T]. Then, given J, L ∈ N, we define the discretization parameters h = (x 0 M − x 0 )/J, k = τ/L, and the number of discrete time levels N = [T/k], which are given by t n = n k, −L ≤ n ≤ N. We begin the description of the numerical method with the initial values, which, in this case, are the initial size discretization X 0 , X 0 j = x 0 + j h, 0 ≤ j ≤ J, the initial condition on the initial grid, U 0 , U 0 j = u 0 (X 0 j ), 0 ≤ j ≤ J, and the resource, that has to be initialized on the interval [−τ, 0], with the values S n = s 0 (t n ), −L ≤ n ≤ 0.
We compute, for 0 ≤ n ≤ N − 1, the approximation at the general level t n+1 : {X n+1 , S n+1 , U n+1 }, , 0 ≤ m ≤ n. We should point out that we design a second-order method based on the modified Euler scheme and on the trapezoidal quadrature rule. This means that we need an intermediate stage in which we compute a first-order approximation: {X n+1, * , S n+1, * , U n+1, * }, X n+1, * = (X n+1, * 0 , X n+1, * 1 , . . . , X n+1, * J+n+1 ), U n+1, * = (U n+1, * 0 , U n+1, * 1 , . . . , U n+1, * J+n+1 ), given by the equations = X n j + k g(X n j , S n , t n ), 0 ≤ j ≤ J + n, S n+1, * = S n + k f (S n , S n−L , Q(X n , γ n (X, S) · U n ), t n ), As soon as we have computed the intermedium quantities {X n+1, * , S n+1, * , U n+1, * } we can obtain the approximations at the advanced time level with In (8) and (9), at each time level, γ n (X, S) and α n (X, S) represent vectors with components γ n j (X, S) = γ(X n j , S n , t n ) and α n j (X, S) = α(X n j , S n , t n ), 0 ≤ j ≤ J + n, respectively, and γ n, * (X, S) and α n, * (X, S) vectors with components γ n, * j (X, S) = γ(X n, * j , S n, * , t n ) and α n, * j (X, S) = α(X n, * j , S n, * , t n ), 0 ≤ j ≤ J + n + 1, respectively. The notation γ n (X, S) · U n , α n (X, S) · U n , γ n, * (X, S) · U n, * and α n, * (X, S) · U n, * represent the componentwise products of the corresponding vectors, and Q is an appropriate second-order quadrature rule, as given below. We observed that, in this procedure, the total number of nodes in the grid increases in one at each time step due to the new node that fluxes from the left boundary into the domain. This feature led us to name this kind of scheme the aggregation grid nodes method (AGN) [12]. The algorithm is explicit and totally straightforward and can be resumed as in the pseudocode given in Box 1.
S n = s 0 (t n ) end % Advancing the solution on time, step by step do n=0, N-1 Formulae in (8) % First stage to X n+1, * , S n+1, * , U n+1, * Formulae in (9) % Second stage to X n+1 , S n+1 , U n+1 % Optional (difference between AGN and SGN methods) Formulae in (10) % Selection procedure in case of SGN method end In the numerical experimentation, we use a more efficient procedure in order to keep the number of nodes constant (and consequently, to not increase the computational cost) at each time step. Toward that goal, we eliminate a chosen node from the grid at each time step. The selection procedure we use in the method depends on the dynamics of the grid, as it was introduced in [13]. We refer to this kind of scheme as the selection grid nodes method (SGN) [12]. The node picked out is the first X n+1 Once it is elected, we rename the nodes to keep the index labels from j = 0 to j = J. Finally, we pay attention to the nonlocal terms; they are discretized by means of a composite quadrature rule, based on the trapezoidal formula, on the grid points X = (X 0 , X 1 , . . . , X ℵ ), ℵ ∈ N,

Convergence Analysis
The convergence analysis here follows the discretization framework developed in [14]. In this framework, suitable discrete normed spaces and operators are introduced to formulate the equations of the numerical method. Then, appropriate properties of consistency and nonlinear stability are established. The numerical method is an adaptation of the one given in [9], and so is the analysis. However, for the sake of completeness, we provide the details of the numerical method formulation, and for the sake of simplicity, we present only those parts of the convergence analysis proofs associated with estimates derived from the discretization of the delay term in the differential equation. That means that the step by step recurrences for the errors at each time level increase the order of the difference equation as the time discretization parameter tends to zero. For the theoretical analysis, we consider the AGN method.
With this aim, we fix T > 0 and assume that the problem (1), (3) has a unique solution u(x, t), x ∈ [x 0 , x M (t)], and s(t), t ∈ [−τ, T], with the following regularity assumptions: and is nonnegative.
Additionally, we assume that there exists a compact neighborhood D of {s (t) , 0 ≤ t ≤ T}, such that and is nonnegative.
and is nonnegative.
and is nonnegative.
and there exists a positive constant C such that g(x 0 , s, t) ≥ C, s, t ≥ 0. In addition, the characteristic curves x(t; t * , x * ) defined by (5), are continuous and differentiable with respect to the initial values (t * , Finally, we suppose that there are compact neighborhoods, D f of {s (t) , −τ ≤ t ≤ T}, and D I of and is nonnegative. Now, we choose L ∈ N and assume that the time discretization parameter, k, takes values in the set Then for k ∈ K, we set N = [T/k] and choose J ∈ N such that h = (x 0 M − x 0 )/J, and the ratio r := k/h is a positive constant fixed throughout the analysis. Then, x 0 M is always the last node in the size grid. For each k ∈ K, we define the space where vectors y 0 , V 0 , . . . , y N , V N , a ∈ A k are used to describe the following approximations: to the inner and right-hand boundary grid nodes, given by y n ∈ R J+n , and to the theoretical solution at the complete grid, V n ∈ R J+n+1 , for each time level t n , 0 ≤ n ≤ N; and to the theoretical solution to the delay differential equation at time levels t n , −L ≤ n ≤ N, provided by a ∈ R L+N+1 . We also consider the space describes the residuals of the discrete solution when it is substituted on the discrete equations that define the numerical method: we discriminate Y 0 , P 0 , A 0 = (A −L , A −L+1 , . . . , A 0 ), given by the approximations to the initial conditions; P 0 , provided by the solution at the left boundary node at each time level; Y n , P n , 0 ≤ n ≤ N, given by the grid nodes and solution to them; and A = (A 1 , A 2 , . . . , A N ), due to the solution to the delay differential equation, as we show below. Both spaces, A k and B k , have the same dimension.
In order to measure the size of the errors, we define Thus, we endow the spaces A k and B k with the following norms. If y 0 , V 0 , . . . , y N , V N , a ∈ A k , Now, for each k ∈ K, we define and we denote x n 0 = x 0 , n ≥ 0. Recall that x(t; t * , x * ) represents the theoretical solution to problem (5), ]. In addition, if u represents the theoretical solution to (1) we define Finally, if s is the theoretical solution to (3) then we define In the following, we introduce the discretization operator. Let R be a positive constant and let we denote by B A k (ũ k , R k p ) ⊂ A k the open ball with centerũ k and radius R k p , 1 < p < 2. Then, we define given by the following equations, first where a 0 = (a −L , a −L+1 , . . . , a 0 ). Vectors X 0 and U 0 represent approximations at t = 0, respectively, to the initial grid nodes and to the theoretical solution at these nodes. Vector S 0 represents an approximation to the initial resource in the interval [−τ, 0]. Second, g(y n j , a n , t n ) + g(y n+1, * j+1 , a n+1, * , t n+1 ) , µ * y n j , a n , t n + µ * y n+1, * j+1 , a n+1, * , t n+1 , f a n , a n−L , Q(y n , γ n (y, a) · V n ), t n + f a n+1, * , a n−L+1 , Q(y n+1, * , γ n+1, * (y, a) · V n+1, * ), t n+1 , 0 ≤ n ≤ N − 1. Where, with the notation introduced in Section 2, y n+1, * j+1 = y n j + k g(y n j , a n , t n ), V n+1, * j+1 = V n j exp −k µ * y n j , a n , t n , a n+1, * = a n + k f a n , a n−L , Q(y n , γ n (y, a) · V n ), t n , We rewrite the quadrature rule employed in (18)-(25) in the next general form Q(X, V) = J+n ∑ l=0 q l (X) V l . Written in this way, we highlight that the number of nodes considered at each time level t n is J + n + 1, counting both the boundary node X n 0 = x n 0 = x 0 , 0 ≤ n ≤ N and the interior grid nodes X n j , 1 ≤ j ≤ J + n, 0 ≤ n ≤ N. This notation is also valid even when we consider quadrature rules whose nodes are, at each time level, a subgrid of X n , 0 ≤ n ≤ N, by asumming that the corresponding weights q l (X) of some of the nodes are zero .
Note that Φ k takes into account all the grid nodes and their corresponding solution values at each time level, and it employs quadrature rules possibly based on a subgrid. IfŨ k = (X 0 , U 0 , the nodes, X n and 0 ≤ n ≤ N and the corresponding values of the solution, U n , 0 ≤ n ≤ N, at such nodes are numerical solutions to the scheme defined by (9) when the composite trapezoidal quadrature rule is used on a subgrid at each time step. On the other hand, the numerical solution to the scheme defined by (9) satisfies (26). Henceforth, C will denote a positive constant, independent of k, h (k = r h), j (0 ≤ j ≤ J + n) and n (−L ≤ n ≤ N), and C may have different values at different places.
As in [12], it is possible to establish the following property that we called (SG) for the selection procedure given by (10) The property (SG) condenses the essential information about the adaptive selection procedure to allow the proof, under the hypotheses (H1)-(H7), of the following general properties of the composite trapezoidal quadrature Q (X, V) based on the subgrids, where q is a positive constant independent of h, k, j (0 ≤ j ≤ J + n) and n (−L ≤ n ≤ N), for 0 ≤ j ≤ J + n, 0 ≤ n ≤ N. (P4) Let R and p be positive constants with 1 < p < 2. The quadrature weights q j are Lipschitz continuous functions on B ∞ (x n , R k p ), 0 ≤ j ≤ J + n, 0 ≤ n ≤ N. (P5) Let R and p be positive constants with 1 < p < 2. If y n , z n ∈ B ∞ (x n , R k p ), V n ∈ B ∞ (u n , R k p ) and a n ∈ B ∞ (s n , R k p ); then when k → 0, 0 ≤ n ≤ N. (P6) Let R and p be positive constants with 1 < p < 2. If y n , z n ∈ B ∞ (x n , R k p ), V n ∈ B ∞ (u n , R k p ) and a n ∈ B ∞ (s n , R k p ); then J+n ∑ i=0 (q i (y n ) − q i (z n )) α (z n i , a n , t n ) V n i ≤ C y n − z n ∞ , when h → 0, 0 ≤ n ≤ N.
Then, the properties (P1)-(P6) allow us to establish that the nonlinear discrete operators describing Φ k are well defined for k ∈ K small enough, as we formulate in the following theorem (we refer to [9] about details of the proof). Theorem 1. Assume that hypotheses (H1)-(H7) hold and that the quadrature rules used in (18)-(25) satisfy the properties (P1)-(P6). If where R is a fixed positive constant and 1 < p < 2, then, for k sufficiently small, 0 ≤ n ≤ N. Furthermore, there exists a positive constant R , independent of k, such that, as k → 0, X n, * ∈ B ∞ (x n , R k p ), S n, * ∈ B ∞ (s n , R k p ) and V n, * ∈ B ∞ (u n , R k p ), and Q(X n, * , γ n, * (X, S) · V n, * ) ∈ D I , Now, we define the local discretization error as and we say that the discretization (14) is consistent if, as k → 0, The following theorem establishes the consistency of the numerical scheme defined by Equations (26).
Another piece of notion that plays an important role in the analysis of the numerical method is the stability with k-dependent thresholds. For k ∈ K, let R k be a real number ( the stability threshold) with 0 < R k < ∞. We say that the discretization (26) is stable forũ k restricted to the thresholds R k , if there exist two positive constants k 0 ∈ K and S ( the stability constant) such that, for any k ∈ K with k ≤ k 0 , the open ball B A k (ũ k , R k ) is contained in the domain of Φ k " and for allṼ k ,W k in that ball, We introduce the following auxiliary result, proved in [9] where the same quadrature rule and the same procedure of selection of the subgrid at each time step were used. Proposition 1. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P6). Let be y n , z n ∈ B ∞ (x n , R k p ), V n , W n ∈ B ∞ (u n , R k p ) and a n , b n ∈ B ∞ (s n , R k p ), where R is a positive constant and 1 < p < 2. Then, as k → 0, |Q(y n, * , γ n, * (y, a) · V n, * ) − Q(z n, * , γ n, * (z, b) · W n, * )| ≤ C ( V n, * − W n, * 1 + y n, * − z n, * ∞ + |a n, * − b n, * |) , Now, we introduce the theorem that establishes the stability of the discretization defined by Equations (26).

Theorem 4.
Assume that (26) is consistent and stable with thresholds R k . If Φ k is continuous in B(ũ k , R k ) and l k B k = o(R k ) as k → 0, then: For k sufficiently small, the discrete Equations (26) possess a unique solution in B(ũ k , R k ). ii) As Finally, we propose the next theorem which establishes the convergence of the numerical method defined by Equations (26).
Theorem 5. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P6). Then, for k sufficiently small, the numerical method defined by Equations (26) has a unique solutionŨ k ∈ B(ũ k , R k ) and The proof of Theorem 5 is immediately derived by means of consistency (Theorem 2), stability (Theorem 3) and Theorem 4. Specifically, if X 0 = x 0 , U 0 = u 0 and S 0 = s 0 , the proposed numerical scheme is second-order accurate.
Next, we also establish an error bound for the differences between the theoretical solution u evaluated at the numerical values of the grid nodes, and the approximation obtained with the numerical method. Theorem 6. Assume that hypotheses (H1)-(H7) hold and that the considered quadrature rules satisfy properties (P1)-(P6). Then, for k sufficiently small, let be defined by u n, † = u(X n 0 , t n ), u(X n 1 , t n ), . . . , u(X n J+n , t n ) ∈ R J+n , 0 ≤ n ≤ N, where X n j , 0 ≤ j ≤ J + n, 0 ≤ n ≤ N, are the grid nodes given in (26). Then, the solutionŨ k satisfies Proof. We only have to consider the triangle inequality |u(X n j , t n ) − U n j | ≤ |u(X n j , t n ) − u(x n j , t n )| + |u(x n j , t n ) − U n j |, 0 ≤ j ≤ J + n, 0 ≤ n ≤ N. The smoothness hypothesis (H1) on u is enough to derive that the first term on the right hand side of the inequality is O(|X n j − x n j |) = O(h 2 + k 2 ), as proved in Theorem 5. Additionally, the second-order rate of convergence proved in such a Theorem shows that the second term is O(h 2 + k 2 ).
The last convergence theorem remains true even if the method employs a cuadrature rule at each time-step in which only a selection of nodes are involved, whenever the subgrid of these cuadrature nodes satisfies the (SG) property. In the case of the selection procedure given by (10), the convergence can be demonstrated in two stages. First, as proven in [12], we can establish that the selection procedure (10) creates subgrids with the property (SG) if the procedure is applied to a grid whose nodes are in the neighborhood of the theoretical ones within a radius R k p . This is just the case, as the convergence Theorem 6 establishes, step by step, if the solutions at each fixed time level are obtained by the discrete operator (14).

Numerical Results
We carried out different simulations with the aim of exploring the behavior of the numerical method for the structured consumer population with a delayed resource model (1), (3). We tried to corroborate the convergence theoretical results with an academical test. It consisted of a test problem with meaningful nonlinearities both from mathematical and biological points of view in which the delay differential equation for the resource employed τ = 1. We used the following size-specific growth, fertility and mortality rates The weight function in (4) was γ(x, z, t) = x 2 , and finally, the function that drove the dynamics of the resource (3) was chosen to be With this choice of data functions, the problem (1), (3) has the following solution: Λ e 29 ρ t/30 1 + Λ e 29 ρ t/30 /κ , where Λ = 24; κ = 5; and ρ and λ are parameters that must be fixed. In the experiment, we employed ρ = 0.1 and λ = 0.3. The knowledge of the exact solution to the problem allowed us to compare it with the numerical solution and to compute the error caused by the numerical approximation. Then, given h and k the discretization parameters in size and time, and the corresponding numerical solution with (X 0 , U 0 , X 1 , U 1 , . . . , X N , U N , S), we computed Note that the exact positions of the grid nodes, x n j , 0 ≤ j ≤ J, 0 ≤ n ≤ N, are unknown, so we compared them with the solution on the computed grid X n j , 0 ≤ j ≤ J, 0 ≤ n ≤ N. In Figure 1, we show the evolution of the computed error in both consumer and resource populations. We observe that the error is produced mainly at the birth (minimum size) and increases with time.
We can also obtain the experimental order of convergence by means of the following well-known formula: We can show the convergence of our method and that it is second-order accurate by means of Table 1. In the test, the initial size interval was taken as [0, x 0 M ], with x 0 M = 0.8, and the numerical integration was carried out on the time interval [0, T], with T = 10. Each column and each row of the table represent a computation with different values of the time and size discretization parameters, respectively. The upper number of each entry in columns two to six of said table represents the error computed from (64) and the lower quantity is the experimental order of convergence. We also remark that the ratio r := k/h was kept fixed as soon as we considered approximations along each diagonal of the table. We realize that the results in Table 1 clearly confirm the expected (theoretically) second-order convergence for these approximations. Different values of parameters ρ and λ confirm the same behavior of the error.

Conclusions
We proposed a numerical method specially adapted to integrate a model that describes the evolution of a size-structured consumer population feeding on a dynamical resource. The model couples a first-order hyperbolic partial differential equation with nonlocal and nonlinear boundary condition and a differential equation driving the dynamics of the resource. The novelty of the model is the appearance of a delay in the resource evolution rate. Numerical methods are the only feasible approach, except in exceptional cases, for the approximation of the solution of the problem and eventually studying the dynamical behavior of the model. The authors have in preparation [15], using this model, a numerical study on the effect of the delay in the resource evolution rate on the dynamics of a population of the Daphnia magna (water flea).
The analyzed numerical method integrated the problem along its characteristic curves. We proved optimal second-order of convergence to the solution; this was confirmed numerically by means of an academic test problem. This is the first convergence analysis to our knowledge for the discretization of this kind of problem.
Author Contributions: All authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded in part by project MTM2017-85476-C2-1-P of the Spanish Ministerio de Economía y Competitividad and European FEDER Funds. VA138G18 of the Consejería de Educación, Junta de Castilla y León, Spain.