The Escalator Boxcar Train Method for a System of Aged-structured Equations in the Space of Measures

The Escalator Boxcar Train (EBT) method is a well known and widely used numerical method for one-dimensional structured population models of McKendrick-von Foerster type. Recently the method, in its full generality, has been applied to aged-structured two-sex population model (Fredrickson-Hoppensteadt model), which consists of three coupled hyperbolic partial differential equations with nonlocal boundary conditions. We derive the simplified EBT method and prove its convergence to the solution of Fredrickson-Hoppensteadt model. The convergence can be proven, however only if we analyse the whole problem in the space of nonnegative Radon measures equipped with bounded Lipschitz distance (flat metric). Numerical simulations are presented to illustrate the results.


Introduction.
1.1. Problem formulation. The Escalator Boxcar Train (EBT) algorithm is a numerical method for solving structured population models. It is based on representing the solution as a sum of masses localised in discrete points and tracing their evolution due to transport and growth. The algorithm has been used in applications for a long time [4], however convergence of the scheme was shown only recently in ref. [1] (without the rate of the convergence), and later in ref. [6] (the rate of the convergence) using the metric space approach proposed in ref. [7,8]. So far the method has been established only for single equation models. The aim of this paper is to derive the EBT algorithm for a system of structured population equations.
We focus on the Fredrickson-Hoppensteadt model, which is a two-sex population model with age as a structure variable. The model was originally formulated in ref. [5] and later developed in ref. [11]. The model consists of three population equations with structure coupled through nonlocal boundary terms and a nonlocal and nonlinear source term. Functions u m (t, x) and u f (t, x) describe the distribution of males and females, while u c (t, x, y) is the number of couples at time t. Structural variables x and y denote the age of males and females, respectively. The following system of equations describes the dynamics of the population of males, females and couples: ∂ t u f (t, y) + ∂ y u f (t, y) + c f (t, y)u f (t, y) = 0, x, y)u c (t, x, y)dxdy, ∂ t u c (t, x, y) + ∂ x u c (t, x, y) + ∂ y u c (t, x, y) + c c (t, x, y)u c (t, x, y) = T (t, x, y), u c (t, x, 0) = u c (t, 0, y) = 0, u c (0, x, y) = u c 0 (x, y).
Functions c m , c f and c c describe the rate of disappearance of individuals. In the case of males and females disappearance equals death, while in the case of couples it equals divorce or death of one of spouses. Functions β m and β f are birth rates of males and females. In the general case the functions may depend nonlinearly on ecological pressure, which is usually modelled by a nonlocal operator depending on the distribution of males, females and couples. However, in this paper we restrict our considerations to the case where disappearance and birth rates depend only on time and on the structure variable.
The marriage function T models the number of new marriages of males and females of age x and y, respectively, at time t. It depends nonlinearly on the distribution of individuals. The choice of this function is a subject of ongoing discussions, see [9], [10], [16], [13], due to the properties like heterosexuality, homogeneity, consistency or competition. In this paper we follow the formulation proposed in ref. [12], namely: Function Θ(x, y) describes the marriage rate of males of age x and females of age y. Notice that u m (t, x) − ∞ 0 u c (t, x, y)dy is the amount of unmarried males EBT FOR A SYSTEM OF AGE-STRUCTURED EQUATIONS 3 and u f (t, y) − ∞ 0 u c (t, x, y)dx is the number of unmarried females. Function h ∈ L 1 (R + ) ∩ L ∞ (R + ) describes the distribution of eligible males on the marriage market. Function g is of the same regularity and describes the distribution of eligible females on the marriage market.
1.2. EBT method for a scalar equation. Before moving to the derivation of the numerical scheme for system (1), we briefly discuss the analytical and numerical results for a simpler model of the McKendrick -von Foerster type [14]: which is a one-sex population model with structure variable x. The function b(t, x) is the rate at which the individuals change their state. In particular case, when x represents age, function b(t, x) equals one. One of the many numerical methods applied to this problem is the Escalator Boxcar Train introduced in ref. [4]. EBT method is based on representing the solution as a sum of masses (cohorts) localised in discrete points and tracing their spatio-temporal evolution using the following algorithm: The procedure consists in solving the system of ordinary differential equations (ODEs) (4)-(5) on a sufficiently short time interval [t k , t k+1 ], iteratively with respect to k = 0, 1, . . .. The details of the EBT algorithm are discussed in Section 2 for system (1) and therefore we omit them at this stage.
As mentioned before, the rate of convergence of the method for a nonlinear version of model (3) has been recently shown in ref. [6]. The results are based on Lipschitz semiflows in metric spaces developed in ref. [3] and applied to structured population models in ref. [8].
which means that µ n (t) is a linear combination of Dirac deltas with masses localised in points x i (t).
The proof of convergence of the ETB scheme is based on the following proposition, the proof of which is a modification of the proof presented in ref. [2].
where ρ is a corresponding metric.
Application of Proposition 1 reduces the proof of the convergence of the measure given by formula (7) to the solution of the model (3) to the estimate of the right hand side of inequality (8). In ref. [6] it was proven analytically and confirmed numerically that the method is of the first order.
1.3. Organisation of the paper. This paper is devoted to derivation of a corresponding EBT method for the two-sex age-structured population model (1). The convergence of the method follows from stability results obtained in ref. [17] and its proof is deferred to a forthcoming paper. Section 2 is devoted to derivation of the EBT algorithm for the two-sex population model. In Section 3, we provide wellposedness results for the structured population model and for the EBT scheme, and discuss a sketch of the convergence proof.
2. Derivation of the EBT algorithm for the two-sex population model. This section is devoted to the formal derivation of the EBT scheme for system (1). We assume that the solution u f , u m i u c to system (1) are continuously differentiable with respect to time and twice time continuously differentiable with respect to the structure variables. Additionally, the solution itself and its first derivative with respect to time and first and second derivatives with respect to the structure variables are assumed to be bounded. The same regularity is assumed on all model coefficients appearing in the underlying system (1).
The first step of the EBT method is grouping the initial distribution of males, females and couples into the so-called cohorts, with the grouping performed with respect to the age of the individuals. While in the case of males and females the groups form one-dimensional intervals including individuals of certain age, in the case of couples we deal with two-dimensional intervals to take into account the age of both partners. Each cohort is characterised by two quantities: the mass and its location. The mass (m m (t), m f (t), m c (t)) is the number of individuals or couples within the cohort at time t, and its location (x m (t), y f (t) and (x c (t), y c (t)), respectively) is the average value of the structure variable over the underlying cohort. Cohorts are evolving in time along the characteristic lines of the underlying problem. We wish to track how masses and locations of cohorts change in time.
Once individuals or couples are assigned to a certain cohort they remain there till disappearance (death or split-up in case of couples). The influx of new couples is described with marriage function, while the influx of new males and females is given in the boundary conditions. Each time step can be treated as an internalisation moment, where a new boundary cohort appears. In case of couples the boundary cohorts are empty and in case of individuals the masses in boundary cohorts result from the boundary conditions.
Following this argument, we start with introducing the initial cohorts. At time t = 0 we impose space mesh points l m i (0) and l f j (0), where i, j = B 0 , . . . , J in such a way that l m B0 (0) = l f B0 (0) = 0. We define J − B 0 initial intervals for males and females The above partitioning requires the condition that the supports of initial functions u m 0 , u f 0 and u c 0 , which are the initial distributions of individuals and couples, are contained in the sums of all corresponding cohorts that is Bearing in mind that the cohorts evolve in time, we write and understand that functions l m i (t) and l f i (t) are characteristics defined by the transport operators, thus they satisfy differential equations d dt l m i (t) = 1 and d dt l f j (t) = 1 and are straight lines: We impose a mesh on the time variable t ∈ [0, T ) in such a way that t 0 = 0 and Mesh points t n , where t n = 1, . . . , N , are called internalisation moments.
Altogether, we have J − B 0 + 1 initial cohorts for males and females and (J − B 0 + 1) 2 initial cohorts of couples.
At each point t n , we define B n = B 0 − n and l Bn (t) = t − t n for t n ≤ t < t n+1 . Now, for t ∈ [t n , t n+1 ), we define internal cohorts: . . , J − 1 and boundary cohorts: This means that at each internalisation moment t n we are adding new boundary cohorts to account for the influx of new individuals. At the same time boundary cohorts from the previous time interval [t n−1 , t n ) became internal ones on [t n , t n+1 ). It is essential to observe that at time step t n we deal with J − B 0 + n + 1 male and female cohorts and (J − B 0 + n + 1) 2 couple cohorts.
For t n ≤ t < t n+1 masses and their locations are defined as follows: , otherwise y f j (t) = where Π m Bn (t) = Ω m Bn (t) xu m (t, x)dx and Π f Bn (t) = Ω f Bn (t) yu f (t, y)dy. Let us notice that we deal with 3 unknowns, thus 3 equations, for each male and female cohort and with 2 unknowns for each couple cohort. Knowing the amount of cohorts explicitly we deduce that at each time step t n we have 3(J − B 0 + n + 1) equations for males and females and 2(J − B 0 + n + 1) 2 equations for couples. Consequently, at t n ≤ t < t n+1 we have 8 + 2(J − B 0 + n) equations more then we had at t n−1 ≤ t < t n .
As our aim is to track the dynamics of the above functions describing masses and their locations, we wish to derive the appropriate system of ordinary differential equations. To shorten the notation, from now on, we assume that t ∈ [t n , t n+1 ) and we write B instead of B n . Let us start with formulating ODEs for masses and locations of male individuals and symmetric results for female individuals and then move to the more demanding case of the couples. We start with differentiating m m i (t), x m i (t), i = B + 1, . . . , J, and boundary functions m m B (t) and Π m B (t): To differentiate the location function x m i (t), for i = B + 1, . . . , J, we check how the first momentx m i (t) = Ω m i (t) xu m (t, x)dx evolves.
Having d dtx m i (t) calculated, it is easy to differentiate the location functions for the internal cohorts: We also need to differentiate the first moment, Π m B (t), in the boundary cohort: Because we need to obtain the closed form of the schemes, we use the following facts that hold for i = B + 1, . . . , J: for f ∈ C 2 (R + ).
While the first one can be easily justified based on the definition of x m i (t) the proof of the second one requires equation (12) and the first order Taylor approximation Furthermore, we observe that +O max Equality (15) holds because Now we are ready to derive the system of equations for the male population: As we mentioned before, one can easily conclude that the system of ODEs for the female population is of similar form, namely |y − y f j (t)| 2 , [T (t, x, y) − c c (t, x, y)u c (t, x, y)] dxdy.
As previously, to differentiate the locations function (x c ij , y c ij )(t) = 1 m c ij (t) Ωij (t) (x, y)u c (t, x, y)dxdy we start with its first moment that is (   (1, 1)u c (t, x, y)dxdy.
As the derivation of the closed form for couples goes in line with reasoning for the male and female populations, we apply equations (14) and (15). We omit details concerning approximation of particular functions, because similar calculations were performed in the case of male population.
For a generic marriage function T (t, x, y) we are not able to obtain the closed form, and most we can conclude is: − (x c ij (t), y c ij (y))c c (t, x c ij (t), y c ij (t))m c ij (t) + (1, 1)m c ij (t) Hence, in further investigations we consider a specific marriage function given in equation (2). To obtain the desired form of ODEs describing dynamics of couples, we need to approximate the integral of the underlying function over cohort Ω c ij (t):     To obtain the approximation of the first order with respect to the structure variables, we apply formulas (12)- (15). Observe that the denominator of function F does not depend on structure variables and hence, there is no need of integrating it over Ω c ij (t). Let us start with the integral of the numerator. It consists of four parts which can be approximated separately: