Numerical integration of an age-structured population model with infinite life span

Applied Mathematics and Computation 434 (2022) 127401 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc Numerical integration of an age-structured population model with infinite life span L.M. Abia a , O. Angulo b , ∗, J.C. López-Marcos a , M.A. López-Marcos a a Departamento de Matemática Aplicada & IMUVA, Facultad de Ciencias, Universidad de Valladolid, Valladolid, Spain b Departamento de Matemática Aplicada & IMUVA, ETS de Ingenieros de Telecomunicación, Universidad de Valladolid, Paseo de Belén, 15 - Campus Miguel Delibes, 47011 Valladolid, Spain a r t i c l e i n f o Article history: Received 19 May 2021 Revised 14 January 2022 Accepted 10 July 2022 Available online 5 August 2022 2010 MSC: 92D25 65M25 65M12 Keywords: Age-structured population Convergence analysis Unbounded life-span Numerical methods Nicholson’s blowflies a b s t r a c t The choice of age as a physiological parameter to structure a population and to describe its dynamics involves the election of the life-span. The analysis of an unbounded lifespan age-structured population model is motivated because, not only new models continue to appear in this framework, but also it is required by the study of the asymptotic behaviour of its dynamics. The numerical integration of the corresponding model is usually performed in bounded domains through the truncation of the age life-span. Here, we propose a new numerical method that avoids the truncation of the unbounded age domain. It is completely analyzed and second order of convergence is established. We report some experiments to exhibit numerically the theoretical results and the behaviour of the problem in the simulation of the evolution of the Nicholson’s blowflies model. ©2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction Population dynamics modelization entails making different choices. Once it is decided that the individuals in the population are non homogeneous, the variables in which these are organized in the population have to be chosen. This is crucial because the variables that determine the structure of the population are largely responsible for its dynamics. Age is one of the most natural and important parameters structuring a population. It plays a main role in demography but, nowadays, it is not restricted to it and also appears in other fields such as ecology, epidemiology and cell growth. Other structuring variables (usually known as size), in general, are intimately dependent on age. A second concern relates to the description of the phenomena that influence the dynamics of the population (mortality, fertility, migration, etc.). It goes back and forth between a discrete age-time setting and a continuous one. We consider the model within the continuous framework. ∗Corresponding author. E-mail addresses: [email protected] (L.M. Abia), oscar[email protected] (O. Angulo), [email protected] (J.C. López-Marcos), [email protected]a.es (M.A. López- Marcos) . https://doi.org/10.1016/j.amc.2022.127401 0 096-30 03/© 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 In a continuous age-structured population model with finite age-span [0 , a † ) , the natural interpretation for the mortality rate function ( μ= μ(a ) ) is through the survival probability to age a , (a ) = exp − a 0 μ(σ) dσ. Then, it is assumed that (a † ) = 0 , in consistency with the model. As a consequence, the mortality rate ( μ) must be unbounded at a † [1,2] . Although the maximum finite age condition has become a benchmark of the study on age-structured problems [2–4] , the analysis of an unbounded age-span under the conditions in Iannelli and Milner [2] is useful. The model we study is described by the following first order hyperbolic partial differential equation, u t + u a = −μ(a, I μ(t ) , t ) u, a > 0 , t > 0 , (1.1) which represents how the dynamics of the individuals in the population evolve. Here independent variables a and trepresent age and time, respectively. The dependent function u (a, t) denotes the population density with respect to the age a at time t, and the function μ(a, I μ(t) , t) ≥0 is the nonnegative age-specific mortality rate function. The partial differential Eq. (1.1) is supplemented with a nonlinear and nonlocal boundary condition which represents the individuals birth law, given by u (0 , t) =  ∞ 0 β(a, I β(t) , t) u (a, t) da , t > 0 , (1.2) where β(a, I β(t ) , t ) ≥0 is the nonegative age-specific fertility rate function. An initial condition u (a, 0) = u 0 (a ) , a ≥0 , (1.3) where u 0 (a ) is a nonnegative function, is also given. In this general nonlinear model, the nonnegative vital functions μand βalso depend on the functionals I μ(t) =  ∞ 0 γμ(a ) u (a, t) da, I β(t) =  ∞ 0 γβ(a ) u (a, t) da , t > 0 , (1.4) to take into account the influence of the age distribution of the population on the life history of the individuals. In these functionals (1.4) , the nonnegative function parameters γμ(a ) and γβ(a ) give an account of the different influence of the individuals in the population, depending on their age, on the mortality and birth rates. The model described by (1.1) –(1.3) was introduced first in Gurtin and MacCamy by Gurtin and McCamy, with the use of vital functions ( μand β) only depending on the age and the total population P (t) =  ∞ 0 u (a, t) da . In that work, they proved the existence and uniqueness of a nonnegative solution u ∈ C 1 (R + ×[0 , T ]) under appropriate compatibility properties between initial and boundary conditions, continuity of the derivatives of mortality and fertility rates, continuity of the initial condition and u  0 ∈ L 1 (R + ) . They also demonstrated the global existence of the solution u for all time when the fertility was uniformly bounded and the existence of a unique equilibrium age-distribution. The study of stability of such equilibrium was only made by considering particular expresions of the vital functions due to the high difficulty in dealing with the general problem. In [6] , the analysis of the problem (1.1) –(1.3) , assuming that the vital functions depend on time and on a weighted population, was addressed in a similar approach as in Gurtin and MacCamy [5] . An extensive summary about the theoretical results for problem (1.1) –(1.3) can be found in Iannelli and Milner [2] . Moreover, we will require additional regularity hypotheses on the vital rates functions and the initial condition to prove convergence of our proposal. We should mention that, for age-structured population models with infinite life span, Webb [7] and Iannelli and Milner [2] considered situations with a birth function of the form β(a ) = e −pa  m j=1 βj a j , that corresponds to a reproductive process in which individuals of the species continue to reproduce as long as they survive. The model (1.1) –(1.3) includes enough information to describe the evolution of a population structured by age. However, the explicit analytical solution of such a general model is unaffordable and can only be obtained in very special cases. Moreover, its asymptotic behaviour requires to be careful with the shape of the age-specific birth and death rates at infinity (see [2,5,6] , where specific form for the rates βand μwere considered). This means that numerical integration is a valuable and feasible alternative. From a numerical point of view, the model (1.1) –(1.3) has been studied in different situations (linear and nonlinear models and applications, for example, in epidemiology) for the last 40 years. We can find an outline of the most significant numerical methods developed so far in Abia et al. [8] . These numerical studies were based on the hypothesis of an initial condition with compact support [0 , A ] , and a finite integration time T , which allowed us to integrate in a finite age interval [0 , A + T ] . Under these assumptions, the numerical schemes provided are accurate. Nevertheless, close attention needs to be paid to the unbounded life span model. On the one hand, new models still arise in this framework as, for example, [9] . On the other hand, the asymptotic study of these models makes the use of an unbonded life span essential: the finite setting does not allowed a far enough numerical integration, this is why numerical methods that approximate the solution in an unbounded life-span should be addressed. The paper is organized as follows: in Section 2 , we reformulate the age-structured population problem (1.1) –(1.3) by introducing a change in both the structured dependent and independent variables that transforms the infinitum life span agestructured problem into a mathematical equivalent structured population model with an artificial structure. In Section 3 , 2 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 we propose a second order numerical method which is described in detail, and we carry out its convergence analysis in Section 4 : it is the highest order method analyzed for this setting. Finally, Section 5 is devoted to the numerical experimentation where we introduce some numerical tests to reveal the experimental second order of convergence, and some features of the long-time integration of the Nicholson’s blowflies model. 2. An artificial size-structured population model The framework we propose to analyze the infinite life-span model involves an artificially structured population model, defined over a bounded state interval, which is equivalent to the original unbounded age-specific one. Then, we integrate numerically the transformed model to obtain the solution to the initial problem. First of all, we will consider a change in both the structured dependent and independent variables to transform the age-structured problem into an equivalent population model structured by an artificial size variable. Thus, the problem is reformulated over a bounded domain of the new artificial structuring variable. For this kind of problem, different numerical methods have been considered in the recent literature [10] . Then, we can apply a numerical method for this wellknown structured population model. Finally, from this approximation we will recover a numerical solution to the original unbounded life-span problem. The new independent variable can be seen as a computational variable in terms of which the numerical procedure, that we will describe in the following section, makes a truncation of the unbounded age interval that changes with the discretization parameter. Therefore, we propose the following change of variables: we relate the age with a new artificial structuring variable x by means of a = α(x ) , x ∈ [0 , 1) , with αan increasing function such that α(0) = 0 and lim x → 1 −α(x ) = + ∞ . Next, the new unknown density function v (x, t) is defined by the identity g(x ) v (x, t) = u (α(x ) , t) , x ∈ [0 , 1) , t ∈ R , where g(x ) = 1 α (x ) , x ∈ [0 , 1) . Thus, the new population density, v (x, t) , satisfies the following structured model v t + (g(x ) v ) x = −μ∗(x, I ∗ μ(t ) , t ) v , 0 < x < 1 , t > 0 , (2.1) g(0) v (0 , t) =  1 0 β∗(x, I ∗ β(t ) , t ) v (x, t ) dx , t > 0 , (2.2) v (x, 0) = 1 g(x ) u 0 (α(x )) , 0 ≤x < 1 , v (1 , 0) = 0 , (2.3) where, we could represent I ∗ μ(t) and I ∗ β(t) as I ∗ μ(t) =  1 0 γ∗ μ(x ) v (x, t) dx, I ∗ β(t) =  1 0 γ∗ β(x ) v (x, t) dx , t > 0 , (2.4) with μ∗(x, I ∗ μ(t ) , t ) = μ(α(x ) , I μ(t ) , t ) ; β∗(x, I ∗ β(t ) , t ) = β(α(x ) , I β(t ) , t ) ; γ∗ μ(x ) = γμ(α(x )) and γ∗ β(x ) = γβ(α(x )) . In this size-structured population model, the variable that structures the population, x , is artificial without biological meaning. Thus, function αis linked with the (size-specific) growth rate function in the associated size-structured population model. In this setting, the selection of function αhas to follow some properties that allow us to define a nonnegative continuous size-specific growth rate function on [0,1], with g(1) = 0 . Then, α∈ C 1 ([0 , 1)) , α (x ) > 0 , x ∈ [0 , 1) (that is, strictly increasing) and lim x → 1 −α (x ) = + ∞ . Note also that lim x → 1 −g(x ) v (x, t) = 0 = lim x → 1 −u (α(x ) , t) . In order to consider the new problem for the new unknown function, we assume lim x → 1 − u 0 (α(x )) g(x ) = 0 . Note that this condition is not very restrictive: an initial data with compact support on [0 , + ∞ ) verifies this condition In this paper, our main objective is to introduce a second order numerical method to solve (1.1) –(1.3) through the discretization of the associated size-structured population problem (2.1) –(2.3) . In the following, we specify the minimum requirements upon the data functions and the solution of the problem for the later convergence analysis. Thus, let be T , σ∈ R + , we define A σ= { f ∈ C 2 ([0 , ∞ )) ;f (n ) (a ) = O(e −σa ) , as a → ∞ , n = 0 , 1 , 2 } and A b σ= { f ∈ C 2 ([0 , ∞ )) ;f bounded and f (n ) (a ) = O(e −σa ) , as a → ∞ , n = 1 , 2 } . Throughout the paper, we assume the following hypotheses, (H1) u ∈ C 2 ([0 , ∞ ) ×[0 , T ]) , is nonnegative and, given t ∈ [0 , T ] , the function u (t) (a ) := u (a, t) , a ≥0 , is such that u (t) (a ) ∈ A σ. 3 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 (H2) β∈ C 2 ([0 , ∞ ) ×D β×[0 , T ]) , is nonnegative and, given (z, t) ∈ D β×[0 , T ] , the function β(z,t) (a ) := β(a, z, t) , a ≥0 , is such that β(z,t) (a ) ∈ A b σ, where D βis a compact neighbourghood of   ∞ 0 γβ(a ) u (a, t) da, 0 ≤t ≤T  . (H3) μ∈ C 2 ([0 , ∞ ) ×D μ×[0 , T ]) , is nonnegative and, given (z, t) ∈ D μ×[0 , T ] , the function μ(z,t) (a ) := μ(a, z, t) , a ≥0 , is such that μ(z,t) (a ) ∈ A b σ, where D μis a compact neighbourghood of   ∞ 0 γμ(a ) u (a, t) da, 0 ≤t ≤T  . (H4) γμ, γβ∈ A b σ, are nonnegative. These assumptions should be enough to get second order convergence with the numerical method proposed. Also, they are reasonable requests to make. On the one hand, under the assumptions (H1) at t = 0 , and (H2)–(H4), the problem (1.1) – (1.3) has a unique nonnegative solution, which is global in time (note that we have included some additional restrictions on the hypotheses given in Gurtin and MacCamy [5] ). Furthermore, it could be proved that u ∈ C 2 (R + ×[0 , T ]) if an additional second order compatibility condition is satisfied. On the other hand, hypotheses (H1) at t = 0 , and (H2)–(H4), ensure u (·, t) ∈ L 1 (R + ) . Once the change of variables is made, the next step is to consider a numerical method to approximate the solution to (2.1) –(2.3) in a fixed time domain [0 , T ] and the bounded size interval [0,1]. Therefore, if we employ a second order numerical scheme, to reach the optimal rate of convergence, we require some regularity conditions on the functions data ( g, μ∗, β∗, γ∗ μand γ∗ β) and the solution v on their corresponding domains similar as (H1)–(H4), that are obtained from them. That is, (H5) v ∈ C 2 ([0 , 1] ×[0 , T ]) , is nonnegative. (H6) g ∈ C 3 ([0 , 1]) , is nonnegative and  1 0 g −1 (x ) dx = ∞ . (H7) β∗∈ C 2 ([0 , 1] ×D β∗×[0 , T ]) , is nonnegative and D β∗is a compact neighbourghood of  1 0 γ∗ β(x ) v (x, t) dx, 0 ≤t ≤T . (H8) μ∗∈ C 2 ([0 , 1] ×D μ∗×[0 , T ]) , is nonnegative and D μ∗is a compact neighbourghood of  1 0 γ∗ μ(x ) v (x, t) dx, 0 ≤t ≤T . (H9) γ∗ μ, γ∗ β∈ C 2 ([0 , 1]) , are nonnegative. All in all, the proposed change of variable has to satisfy (H10) α∈ C 4 ([0 , 1)) , α(0) = 0 , α (x ) > 0 , x ∈ [0 , 1) , lim x → 1 −α(x ) = ∞ , lim x → 1 −α (x ) = ∞ and lim x → 1 α (x ) u 0 (α(x )) = 0 . 3. Numerical approximation In this section we introduce a numerical method for the original age-structured model (1.1) –(1.3) through the transformation that converts it into the size-structured model (2.1) –(2.3) . The new model could be solved with the usual numerical techniques employed to obtain the solution of size-structured models in the case of a bounded size interval (see Abia et al. [10] ). We will approach the solution of the size-structured problem and, finally, we return to the solution of the original problem. In this work, we will integrate (2.1) –(2.3) by means of a scheme based on the numerical approximation to the solution along the characteristic curves. It is a natural and efficient way to discretize (2.1) . To this end, we rewrite it as v t + g(x ) v x = −μ∗(x, I ∗ μ(t ) , t ) + g  (x ) v , 0 < x < 1 , t > 0 . (3.1) Now, we denote by x (t ;t ∗, x ∗) the characteristic curve associated with Eq. (3.1) which takes the value x ∗at time t ∗. It is the solution to the following initial value problem x  (t ;t ∗, x ∗) = g(x (t ;t ∗, x ∗)) , t ≥t ∗, x (t ∗;t ∗, x ∗) = x ∗. (3.2) Then, we define w (t ;t ∗, x ∗) = v (x (t ;t ∗, x ∗) , t ) , t ≥t ∗, that is the solution of w  (t ;t ∗, x ∗) = −μ∗(x (t ;t ∗, x ∗) , I ∗ μ(t ) , t ) + g  (x (t ;t ∗, x ∗)) w (t ;t ∗, x ∗) , t ≥t ∗, w (t ∗;t ∗, x ∗) = v (x ∗, t ∗) . (3.3) 4 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 The solution along the characteristic curves can be written as w (t ;t ∗, x ∗) = v (x ∗, t ∗) exp − t t ∗μ∗(x (τ;t ∗,x ∗) ,I ∗ μ(τ) ,τ) + g  (x (τ;t ∗,x ∗)) dτ,t ≥t ∗(3.4) Note that the viability of the change of variables proposed in the previous section, is independent of the particular choice of the function αwhich relates the age variable with the new structuring variable. However, the efficiency of the numerical method used for the approximation of the solution of (2.1) –(2.3) , and the corresponding analysis of convergence, would be affected by this election. Then, when considering a particular choice of α, we will take into account that the problem obtained after the change of variables (2.1) –(2.3) is a typical size-structured population model. Here we propose to employ α(x ) = −1 K α log (1 −x ) , x ∈ [0 , 1) , (where K αis a positive parameter that will be fixed), which provides the well-known Bertalanffy’s growth rate [11] g(x ) = K α(1 −x ) , x ∈ [0 , 1] . (3.5) For the numerical integration on the bounded time interval [0 , T ] , we introduce N ∈ N as the number of time steps. The time discretization parameter is given by k = T /Nand the time levels are represented by t n = n k , 0 ≤n ≤N. Then, we introduce the grid over the size interval. We use the natural grid in which the grid points are computed by means of the solution to (3.2) . Actually, it can be made analytically due to the simple expression of the growth-rate function (3.5) . Therefore, we define x j , j ≥0 as x j = 1 −e −K αkj , j = 0 , 1 , 2 , . . . (3.6) The numerical method requires a finite number of grid nodes { x j } J+1 j=0 . So, we denote as J, the first positive integer such that 1 −x J ≤K 1 k , with K 1 a suitable positive number independent on k . Taking into account (3.6) , Jis the first integer that satisfies J > −log (K 1 k ) K αk . Finally, x J+1 = 1 . This choice of α(and g) is enough to ensure the existence of x J [12] . On the other hand, we should note that this grid satisfies x j −x j−1 ≤C k , j = 1 , . . . , J + 1 , for a fixed positive constant C, independent on k , which is a crucial property for the use of a quadrature, based on these nodes, in order to approximate the nonlocal terms. We obtain the approximations to the solution at the grid nodes by means of an explicit second order numerical method. We refer to the grid point x j by a subscript jand to the time level t n by a superscript n . Let V n j be a numerical approximation to v n j = v (x j , t n ) , 0 ≤j ≤J + 1 , 0 ≤n ≤N. For simplicity, we use the vectorial representation V n = (V n 0 , V n 1 , . . . , V n J+1 ) , 0 ≤n ≤ N. Therefore, starting from an approximation to the initial data (2.3) , for example, the grid restriction of the function v (x, 0) , the numerical solution at a new time level t n +1 , 0 ≤n ≤N −1 , is described in terms of the previous one. Such a general time step is given by the following formulae, that we obtain from the discretization of (3.4) V n +1 , ∗ j+1 = V n j exp −k μ∗x j , Q ∗ k γ∗ μ·V n , t n −K α, 0 ≤j ≤J −1 , (3.7) V n +1 j+1 = V n j exp −k 2 μ∗x j , Q ∗ k ( γ∗ μ·V n ) , t n  + μ∗x j+1 , Q ∗ k ( γ∗ μ·V n +1 , ∗) , t n +1 −2 K α, 0 ≤j ≤J −1 , (3.8) and the boundary condition (2.2) , V n +1 0 = Q ∗ k ( β∗,n +1 (V ) ·V n +1 ) /K α, (3.9) where ( γ∗ μ) j = γ∗ μ(x j ) , ( γ∗ β) j = γ∗ β(x j ) , β∗,n (V ) j = β∗(x j , Q ∗ k ( γ∗ β·V n ) , t n ) , 0 ≤j ≤J, 1 ≤n ≤N. Finally, we fix V n, ∗ J+1 = V n J+1 = 0 , 0 ≤n ≤N. In (3.7) –(3.9) , Q ∗ k represents a quadrature rule to approximate the integrals over the size interval. In order to save computational effort when computing the quadrature approximations, initially we fix a subgrid of the natural grid { x j l } M+1 l=0 , satisfying x j 0 = 0 , x j M+1 = 1 , and C 0 k ≤x j l+1 −x j l ≤C 1 k , 0 ≤l ≤M + 1 , where C 0 and C 1 are positive constants independent of k , for further details we refer to Angulo and López-Marcos [12] . Then, to obtain an explicit numerical method, we choose a quadrature formula based the composite trapezoidal rule, on the subintervals defined by the subgrid, combined with the use of a right-hand rectangular rule on the first interval and a left-hand rectangular rule on the last one. Therefore, the resulting formula is given by Q ∗ k (V ) = x j 1 V j 1 + M−1  l=1 x j l+1 −x j l 2 V j l + V j l+1 + (1 −x j M ) V j M , (3.10) 5 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 where V = V 0 , V 1 , . . . , V J+1 . It keeps the second order of accuracy. Lastly, γ∗ μ·V n , γ∗ μ·V n +1 , ∗and β∗,n +1 (V ) ·V n +1 represent the componentwise product of the corresponding vectors. By this election, V n +1 0 is always well-defined by formula (3.9) as the boundary nodal value V n +1 0 does not contribute to the quadrature approximation. We emphasize that, even the quadrature rule were chosen closed at the end x j 0 = 0 , again formula (3.9) will provide implicitly an approximation V n +1 0 , for enough small value of the discretization parameter k , because then, the coefficient of the contribution of V n +1 0 to the quadrature approximation is decreasing to zero with k . The formulae previously presented could be described in terms of age by returning to the original variables, and would represent an approximation to the age-structured problem. The points in the age-grid are given by a j = α(x j ) , 0 ≤j ≤J. So, it is easy to show that a j = j k , 0 ≤j ≤J, which are equally spaced (and correlated with time as usual). The condition introduced to finish the computations at the natural grid induces that the last point in the age-grid satisfies a J > −1 K α log (K 1 k ) . Let U n j be a numerical approximation to u (a j , t n ) , 0 ≤j ≤J, 0 ≤n ≤Nand U n = (U n 0 , U n 1 , . . . , U n J ) , 0 ≤n ≤N. Thus the numerical method is completed with the following formula U n j = K αe −K αkj V n j , 0 ≤j ≤J, 0 ≤n ≤N. (3.11) 4. Convergence In this section we analyze the numerical method (3.6) –(3.10) that approaches the solution to the associated sizestructured model (2.1) –(2.3) presented in Section 3 . We finish with the convergence of the transformed values (3.11) to the solution to the age-structured model (1.1) –(1.3) . We begin with a Lemma that establishes the smoothness properties of the associated size-structured problem. Lemma 1. Assume that hypotheses (H1)–(H4) hold and α(x ) = −1 K α log (1 −x ) , x ∈ [0 , 1) . Then, the functions data and the solution to (2.1) –(2.3) satisfy (H5)–(H9) when 3 K α< σ. Proof. Due to the hypotheses (H1)–(H4) and the smoothness of α, we ensure such a smoothness properties when we restrict them to the interval [0,1). We extend the smoothness properties on the interval [0,1] by means of the definition of each function and derivatives at x = 1 as the corresponding limit if there exists. In the case of the solution to (2.1) –(2.3) lim x → 1 − v (x, t) = lim a →∞ u (a, t) K αe −K αa = 0 , t > 0 , then v (1 , t) = 0 , t > 0 . The partial derivatives of v with respect to time at x = 1 are computed with, lim x → 1 − ∂ n v ∂t n (x, t) = lim a →∞ ∂ n u ∂t n (a, t) K αe −K αa = 0 , t > 0 , n = 1 , 2 , then ∂ n v ∂t n (1 , t) = 0 , t > 0 , n = 1 , 2 . The first order partial derivative with respect to size at x = 1 satisfies lim x → 1 − v x (x, t) = lim x → 1 − u a (α(x ) , t) −u (α(x ) , t) g  (x ) g 2 (x ) = lim a →∞ u a (a, t) + K αu (a, t) K 2 αe −2K αa = 0 , t > 0 , then v x (1 , t) = 0 , t > 0 . And finally, the second orden derivative with respect to size at x = 1 is given by lim x → 1 − v xx (x, t) = lim x → 1 − u aa (α(x ) , t) + 3 K αu a (α(x ) , t) + 2 K 2 αu (α(x ) , t) g 3 (x ) = lim a →∞ u aa (a, t) + 3 K αu (a, t) + 2 K 2 αu (a, t) K 3 αe −3K αa = 0 , t > 0 , then v xx (1 , t) = 0 , t > 0 . These ensure that v ∈ C 2 ([0 , 1] ×[0 , T ]) . With similar arguments we derive the smoothness properties of μ∗, β∗, γ∗ μand γ∗ β.  The convergence result will be the consequence of a study of the consistency and the nonlinear stability. We shall employ the discretization framework introduced by López-Marcos et al. [13] to analyze the numerical method described by (3.6) – (3.10) . We consider that the discretization parameter k takes values in the set H = { k > 0 : k = T /N, N ∈ N } , and for each k ∈ Hwe define the vector spaces X k = (R J+1 ) N+1 and Y k = (R J+1 ) ×R N ×(R J ) N . If V = V 0 , V 1 , . . . , V J ∈ R J+1 , and  V  ∞ ,J+1 = max 0 ≤j≤J | V j | , we define the following norm on X k , if (V 0 , V 1 , . . . , V N ) ∈ X k  (V 0 , V 1 , . . . , V N )  X k = max 0 ≤n ≤N  V n  ∞ ,J+1 . 6 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 In addition, if Z ∈ R N , W ∈ R J , and  Z  ∞ ,N = max 1 ≤n ≤N | Z n | ,  W  ∞ ,J = max 1 ≤j≤J | W j | , then for (Z 0 , Z 0 , Z 1 ∗, . . . , Z N ∗) ∈ Y k , we define  (Z 0 , Z 0 , Z 1 ∗, . . . , Z N ∗)  Y k =  Z 0  ∞ ,J+1 +  Z 0  ∞ ,N + N  n =1 k  Z n ∗ ∞ ,J , where Z n ∗= (Z n 1 , Z n 2 , . . . , Z n J ) , 1 ≤n ≤N. Finally, if V = (V 0 , V 1 , . . . , V J ) , we also consider  V  1 ,M+1 = M  l=0 k | V j l | , where { x j l } M+1 l=0 are the nodes that define the subgrid. For each k ∈ H, we consider the element v k = (v 0 , v 1 , . . . , v N ) ∈ X k , v n = (v n 0 , v n 1 , . . . , v n J ) ∈ R J+1 , v n j = v (x j , t n ) , 0 ≤j ≤J, 0 ≤n ≤N, where v is the solution of (2.1) –(2.3) . Let R k be a fixed positive constant and denote by B(v k , R k ) ⊂X k , the open ball with center v k and radious R k . Next, we introduce the mapping k : B(v k , R k ) → Y k defined by the equations k (W 0 , W 1 , . . . , W N ) = (Z 0 , Z 0 , Z 1 ∗, . . . , Z N ∗) , (4.1) Z 0 = W 0 −V 0 , (4.2) Z n +1 j+1 = 1 k W n +1 j+1 −W n j exp −k 2 μ∗x j , Q ∗ k ( γ∗ μ·W n ) , t n + μ∗x j+1 , Q ∗ k ( γ∗ μ·W n +1 , ∗) , t n +1 −2 K α, 0 ≤j ≤J −1 , (4.3) Z n +1 0 = W n +1 0 −Q ∗ k ( β∗,n +1 (W ) ·W n +1 ) /K α, (4.4) 0 ≤n ≤N −1 , where W n +1 , ∗ j+1 = W n j exp −k μ∗x j , Q ∗ k ( γ∗ μ·W n ) , t n −K α, (4.5) 0 ≤j ≤J −1 , 0 ≤n ≤N −1 . Then V k = (V 0 , V 1 , . . . , V N ) ∈ X k , is a solution of the scheme (3.6) –(3.10) if and only if k (V k ) = 0 . (4.6) In the following we study the consistency, stability, existence of solutions and convergence of our scheme. In our analysis we shall assume the hypotheses (H5)–(H9). Hence we choose the radius R k of the open ball B(v k , R k ) ⊂X k such that, for k sufficiently small, if (V 0 , V 1 , . . . , V N ) ∈ B(v k , R k ) then Q ∗ k ( γ∗ μ·V n ) , Q ∗ k ( γ∗ μ·V n, ∗) ∈ D μ∗, and Q ∗ k ( γ∗ β·V n ) , Q ∗ k ( γ∗ β·V n, ∗) ∈ D β∗, (4.7) 0 ≤n ≤N. At this point, we should prove the second order of accuracy of the cuadrature rule given in (3.10) , which is straightforward. In order to demonstrate (4.7) , we obtain, as k → 0 , |Q ∗ k γ∗ μ·V n −I ∗ μ( t n ) | ≤|Q ∗ k γ∗ μ·V n −Q ∗ k γ∗ μ·v n | + |Q ∗ k γ∗ μ·v n −I ∗ μ( t n ) | ≤ γ∗ μ ∞ R k + o ( 1 ) , 0 ≤n ≤N, and a similar inequality holds for Q ∗ k ( γ∗ β·V n ) . In the case of Q ∗ k ( γ∗ μ·V n, ∗) ∈ D μ∗, we prove, with the use of the formula in (4.5) , that | V n, ∗ j −v n j | ≤exp ( K αk )  V n −1 j−1 exp −k μ∗x j−1 , Q ∗ k γ∗ μ·V n −1 , t n −1 −v n −1 j−1 exp − t n t n −1 μ∗x τ;t n −1 , x j−1 , I ∗ μ( τ) , τdτ ≤exp ( K αk ) | V n −1 j−1 −v n −1 j−1 | exp −k μ∗x j−1 , Q ∗ k γ∗ μ·V n −1 , t n −1  + | v n −1 j−1 | | exp −k μ∗x j−1 , Q ∗ k γ∗ μ·V n −1 , t n −1 −exp −k μ∗x j−1 , I ∗ μt n −1 , t n −1 | + | v n −1 j−1 | | exp −k μ∗x j−1 , I ∗ μt n −1 , t n −1 −exp − t n t n −1 μ∗x τ;t n −1 , x j−1 , I ∗ μ( τ) , τdτ|  ≤exp ( K αk ) ( R k + o ( 1 ) ) , and then |Q ∗ k γ∗ μ·V n, ∗−I ∗ μ( t n ) | ≤|Q ∗ k γ∗ μ·V n, ∗−Q ∗ k γ∗ μ·v n | + |Q ∗ k γ∗ μ·v n −I ∗ μ( t n ) | ≤ γ∗ μ ∞ exp ( K αk ) ( R k + o ( 1 ) ) + o ( 1 ) , 0 ≤n ≤N, and a similar inequality holds for Q ∗ k ( γ∗ β·V n, ∗) ∈ D β∗. Therefore, Eq. (4.7) is derived if R k = o ( 1 ) , as k → 0 . From now on, C will denote a positive constant which is independent of k , n (0 ≤n ≤N) and j(0 ≤j ≤J) ; Chas possibly different values in different places. 7 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 4.1. Consistency We define the local discretization error as l k = k (v k ) ∈ Y k , (4.8) and we say that the discretization (4.1) is consistent if, as k → 0 , lim k → 0  k (v k )  Y k = lim k → 0  l k  Y k = 0 . (4.9) The next theorem establishes the consistency of the numerical scheme (3.6) –(3.10) . Theorem 2. Assume hypotheses (H5)–(H9) on the functions data and the solution to (2.1) –(2.3) . Then, as k → 0 , the local discretization error satisfies  k (v k )  Y k =  v 0 −V 0  ∞ + O (k 2 ) . (4.10) Proof. We begin with the auxiliary values. They correspond with a first order approximation at the next time level. We employ the first order accuracy of the rectangular quadrature rule, the second order accuracy of Q ∗ k , hypotheses (H5), (H8) and (H9), to arrive at | v n +1 j+1 −v n +1 , ∗ j+1 | = | v n j | e K αk  exp − t n +1 t n μ∗x τ;t n , x j , I ∗ μ( τ) , τdτ−exp −k μ∗x j , Q ∗ k γ∗ μ·v n , t n  ≤C  t n +1 t n μ∗x τ;t n , x j , I ∗ μ( τ) , τdτ−k μ∗x j , I ∗ μ( t n ) , t n  + C k | μ∗x j , I ∗ μ( t n ) , t n −μ∗x j , Q ∗ k γ∗ μ·v n , t n | ≤C k 2 + C k | I ∗ μ( t n ) −Q ∗ k γ∗ μ·v n | ≤C k 2 . (4.11) Next, we need to bound the goodness of the approach when we combine Q ∗ k with the auxiliary values. Then, the order of convergence of Q ∗ k , hypotheses (H5) and (H9), allow us to derive | I ∗ μt n +1 −Q ∗ k γ∗ μ·v n +1 , ∗| ≤| I ∗ μt n +1 −Q ∗ k γ∗ μ·v n +1 | + |Q ∗ k γ∗ μ·v n +1 −v n +1 , ∗| ≤C k 2 . (4.12) Again, the second order of accuracy of the trapezoidal and Q ∗ k rules, hypotheses (H5), (H8) and (H9), and (4.11) and (4.12) imply that  t n +1 t n μ∗(x (τ;t n ,x j ) ,I ∗ μ(τ) ,τ) dτ−k 2 μ∗(x j , Q ∗ k ( γ∗ μ·v n ) , t n )+μ∗(x j+1 , Q ∗ k ( γ∗ μ·v n +1 , ∗) , t n +1 )  ≤ t n +1 t n μ∗(x (τ;t n ,x j ) ,I ∗ μ(τ) ,τ) dτ−k 2 μ∗(x j , I ∗ μ(t n ) , t n ) + μ∗(x j+1 , I ∗ μ(t n +1 ) , t n +1 )  + k 2 μ∗(x j , I ∗ μ(t n ) , t n ) −μ∗(x j , Q ∗ k ( γ∗ μ·v n ) , t n )  + k 2 μ∗(x j+1 , I ∗ μ(t n +1 ) , t n +1 ) −μ∗(x j+1 , Q ∗ k ( γ∗ μ·v n +1 , ∗) , t n +1 )  ≤C k 3 + C k I ∗ μ(t n ) −Q ∗ k ( γ∗ μ·v n ) + Ck I ∗ μ(t n +1 ) −Q ∗ k ( γ∗ μ·v n +1 , ∗)  ≤C k 3 . (4.13) Now, we denote k (v k ) = (L 0 , L 0 , L 1 , . . . , L N ) , thus (4.3) , (H5) and (4.13) allow us to obtain, | L n +1 j+1 | = 1 k  v n +1 j+1 −v n j exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·v n ) , t n ) + μ∗(x j+1 , Q ∗ k ( γ∗ μ·v n +1 , ∗) , t n +1 ) −2 K α = 1 k | v n j | e K αk  exp − t n +1 t n μ∗(x (τ;t n ,x j ) ,I ∗ μ(τ) ,τ) dτ −exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·v n ) , t n )+μ∗(x j+1 , Q ∗ k ( γ∗ μ·v n +1 , ∗) , t n +1 )  = O (k 2 ) , (4.14) 8 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 0 ≤j ≤J −1 , 0 ≤n ≤N −1 . Then, we take into account (4.4) , the properties of Q ∗ k , hypotheses (H5), (H7) and (H9), to arrive at | L n 0 | = | v n 0 −Q ∗ k ( β∗,n (v n ) ·v n ) /K α| = 1 K α 1 0 β∗(x, I ∗ β(t n ) , t n ) v (x, t n ) dx −Q ∗ k ( β∗,n (v n ) ·v n )  ≤1 K α 1 0 β∗(x, I ∗ β(t n ) , t n ) v (x, t n ) dx −Q ∗ k ( β∗,n I ·v n )  + 1 K αQ ∗ k β∗,n I −β∗,n (v n ) ·v n  ≤C k 2 + C I ∗ β(t n ) −Q ∗ k ( γ∗ β·v n )  ≤C k 2 , (4.15) 1 ≤n ≤N, where we have introduced the following auxiliar notation ( β∗,n I ) j = β∗(x j , I ∗ β(t n ) , t n ) , 0 ≤j ≤J + 1 . Thus, the combination of (4.14) and (4.15) proves (4.10) .  4.2. Stability Another notion that plays an important role in the analysis of our numerical method is stability . For each k ∈ H, let R k be a positive real number or + ∞ (the stability threshold) . We say that the discretization (4.1) is stable for v k restricted to the thresholds R k , if there exist two positive constants k 0 and S(the stability constant) such that, for any k in Hwith k ≤k 0 the open ball B(v k , R k ) is contained in the domain of k and for all V k , W k ∈ B(v k , R k )  V k −W k  X k ≤S  k (V k ) −k (W k )  Y k . Theorem 3. Assume the hypotheses of Theorem 2 and let R k be a positive constant such that (4.7) holds. Then the discretization (4.2) –(4.6) is stable for v k with threshold R k = R k . Proof. Let (V 0 , V 1 , . . . , V N ) , (W 0 , W 1 , . . . , W N ) be in the ball B(v k , R k ) of the space X k . We set E n = V n −W n ∈ R J+1 , 0 ≤n ≤N; k (V 0 , V 1 , . . . , V N ) = (Z 0 , Z 0 , Z 1 , . . . , Z N ) , k (W 0 , W 1 , . . . , W N ) = (S 0 , S 0 , S 1 , . . . , S N ) . There exists a positive constant Csuch that, for k sufficiently small, |Q ∗ k ( γ∗ β·(V n −W n )) | ≤C V n −W n  1 ,M+1 , (4.16) |Q ∗ k ( γ∗ μ·(V n −W n )) | ≤C  V n −W n  1 ,M+1 , (4.17) |Q ∗ k (V n ·W n ) | ≤C V n  ∞ ,J+1  W n  1 ,M+1 , (4.18) 0 ≤n ≤N. By (4.3) , we can write E n +1 j+1 = k (Z n +1 j+1 −S n +1 j+1 ) + e K αk V n j exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·V n ) , t n ) + μ∗(x j+1 , Q ∗ k ( γ∗ μ·V n +1 , ∗) , t n +1 )  −W n j exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·W n ) , t n ) + μ∗(x j+1 , Q ∗ k ( γ∗ μ·W n +1 , ∗) , t n +1 )  = k (Z n +1 j+1 −S n +1 j+1 ) + e K αk E n j exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·V n ) , t n ) + μ∗(x j+1 , Q ∗ k ( γ∗ μ·V n +1 , ∗) , t n +1 )  + W n j exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·V n ) , t n ) + μ∗(x j+1 , Q ∗ k ( γ∗ μ·V n +1 , ∗) , t n +1 )  −exp −k 2 μ∗(x j , Q ∗ k ( γ∗ μ·W n ) , t n ) + μ∗(x j+1 , Q ∗ k ( γ∗ μ·W n +1 , ∗) , t n +1 ) , (4.19) 9 L.M. Abia, O. Angulo, J.C. López-Marcos et al. Applied Mathematics and Computation 434 (2022) 127401 equilibria and, when applicable, it allows us to obtain a representation of stable limit cycles. We have applied this approach to the Nicholson’s blowflies model and we have obtained similar results as predicted by other authors, however, with a significant difference as the stable limit cycle approached is in accordance with original experimental values. This problem could have more interesting dynamics than shown in previous works which remain unexplored. Finally, the numerical method proposed to study this model invites us to rethink and question other problems that originally we proposed with unbounded life-span but they were simplified with a suitable truncated age-interval [9,17] . Acknowledgements Authors are very grateful to the referees for their careful reading of the original manuscript and their aid to improve it. 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, grant PID2020-113554GB-I00/AEI/10.13039/50110 0 011033 grant RED2018-102650-T funded by MCIN/AEI/10.13039/50110 0 011033, and research grant VA193P20 of the Junta de Castilla y León, Spain, and European FEDER funds (EU) and VA138G18 of the Junta de Castilla y León. References [1] M. Iannelli, F. Milner, On the approximation of the Lotka–McKendrick equation with finite life-span, J. Comput. Appl. Math. 136 (2001) 245–254 . [2] M. Iannelli, F. Milner, The Basic Approach to Age-Structured Population Dynamics, Springer, Dordrecht, The Netherlands, 2017, doi: 10.1007/ 978- 94- 024- 1146- 1 . [3] L. Abia, O. Angulo, J. López-Marcos, M. López-Marcos, Approximating the survival probability in finite life-span population models, J. Comput. Appl. Math. 330 (2018) 783–793, doi: 10.1016/j.cam.2017.05.004 . [4] D. Breda, M. Iannelli, S. Maset, R. Vermiglio, Stability analysis of the Gurtin–MacCamy model, SIAM J. Numer. Anal. 46 (2) (2008) 980–995, doi: 10. 1137/070685658 . [5] M. Gurtin, R. MacCamy, Nonlinear age-dependent population dynamics, Arch. Ration. Mech. Anal. 54 (1974) 281–300 . [6] C. Rorres, Local stability of a population with density-dependent fertility, Theor. Popul. Biol. 16 (3) (1979) 283–300, doi: 10.1016/0040-5809(79)90018-2 . [7] G. Webb, Theory of Nonlinear Age-Dependent Population Dynamics, Marcel Dekker, New York, 1985 . [8] L. Abia, O. Angulo, J. López-Marcos, Age-structured population models and their numerical solution, Ecol. Model. 188 (1) (2005) 112–136, doi: 10.1016/ j.ecolmodel.20 05.05.0 07 . [9] M. Adimy, O. Angulo, F. Crauste, J. López-Marcos, Numerical integration of a mathematical model of hematopoietic stem cell dynamics, Comput. Math. Appl. 56 (3) (2008) 594–606, doi: 10.1016/j.camwa.2008.01.003 . [10] L. Abia, O. Angulo, J. López-Marcos, Size-structured population dynamics models and their numerical solutions, Discrete Contin. Dyn. Syst. - Ser. B 4 (4) (2004) 1203–1222, doi: 10.3934/dcdsb.2004.4.1203 . [11] L.V. Bertalanffy, A quantitative theory of organic growth (inquiries on growth laws. II), Hum. Biol. 10 (1938) 181–213 . [12] O. Angulo, J. López-Marcos, Numerical schemes for size-structured population equations, Math. Biosci. 157 (1–2) (1999) 169–188, doi: 10.1016/ s0 025-5564(98)10 081-0 . [13] 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 . [14] A. Nicholson, The self-adjustement of population to change, Cold Spring Harbor Symp. Quant. Biol. 22 (1957) 153–173 . [15] D. Sulsky, Numerical solution of structured population models I. Age structure, J. Math. Biol. 31 (8) (1993) 817–848 . [16] W.S.C. Gurney, R.M. Nisbet, J.H. Lawton, The systematic formulation of tractable single-species population models incorporating age structure, J. Anim. Ecol. 52 (2) (1983) 479–495, doi: 10.2307/4567 . [17] O. Angulo, J. López-Marcos, M. López-Marcos, A numerical integrator for a model with a discontinuous sink term: the dynamics of the sexual phase of monogonont rotifera, Nonlinear Anal. RWA 6 (5) (2005) 935–954, doi: 10.1016/j.nonrwa.20 04.11.0 07 . 16