Finite difference schemes for a structured population model in the space of measures

Abstract: We present two finite-difference methods for approximating solutions to a structured population model in the space of non-negative Radon Measures. The first method is a first-order upwind-based scheme and the second is high-resolution method of second-order. We prove that the two schemes converge to the solution in the Bounded-Lipschitz norm. Several numerical examples demonstrating the order of convergence and behavior of the schemes around singularities are provided. In particular, these numerical results show that for smooth solutions the upwind and high-resolution methods provide a first-order and a second-order approximation, respectively. Furthermore, for singular solutions the second-order high-resolution method is superior to the first-order method.


Introduction
Transport equations have been rigorously studied for many years [1][2][3]. These equations are particularly useful in modeling the dynamics of a population with some physiological structure such as age or size [2,[4][5][6]. Classically, these models take place in a smooth or integrable setting [2,[7][8][9]. These frameworks help model populations using densities with respect to the structure (e.g., size, age) variable. However, it is often useful to consider populations with concentrated masses. As proposed in [5], these cases can be handled by considering the framework of measures, particularly, of Radon measures. The space of Radon measures is of particular interest as it contains both distributions with densities (i.e., absolutely continuous with respect to the Lebesgue measure) and those with concentrated mass (i.e., Dirac measures). This allows for joint analysis of both discrete and continuous distributions in the structure (e.g., size, age) variable .
In this article, we are interested in the following structured population model: Here, given a Borel subset A ⊂ R + , µ(t)(A) represents the number of individuals at time t of structure (e.g., size, age) x in A, and the functions g and d represent the growth and death rate of individuals at time t of structure x, respectively. Likewise, the function β represents the reproduction rate of these individuals. More precisely, at time t and distribution µ(t) of the structure, an individual with structure x produces offspring at rate β(t, µ(t))(x). Thus, D dx µ(0) denotes the Radon-Nikodym derivative of µ(t) with respect to the Lebesgue measure, dx, at the point x = 0. The first equation in the model (1.1) describes how the number of individuals with structure x, µ(t)(x) informally, change in time t due to the combination of two effects: the transport term ∂ x (g(t, µ)µ) which moves the distribution µ at velocity g and the death rate which removes individuals from the system at rate d. The second equation models the inflow of individuals at the boundary due to birth. The third equation simply states the regularity of the initial condition.
The model functions, g, d, β, are assumed to be influenced by both time, t, and solution to the population model, µ(t). In applications (e.g., see [6,10,11]), it is common to choose β, g and µ to depend on a weighted mean of the population in the following form: β(t, µ)(x) = B t, x, R + K B (y)dµ(y) , g(t, µ)(x) = G t, x, R + K G (y)dµ(y) , d(t, µ)(x) = D t, x, R + K D (y)dµ(y) , for given maps B, D, G : [0, T ] × R + × R + → R + and K B , K G , K D : R + → R + . Common physically motivated model functions utilize Beverton-Holt type [12] or Ricker type [13] nonlinearities with respect to the weighted mean of the population and of a Von Bertalanffy type [14] model with respect to structure x. For example, utilizing Ricker type nonlinearity with a Von Bertalanffy model for size [10] with individual maximum size x max and r being a parameter related to the inherent growth rate at size 0 and low population levels, the function G may be chosen to be G = r(x max − x) exp(− R + K G (y)dµ(y)) x ≤ x max 0.
x > x max . (1.4) Well-posedness of such a model was studied in [15] and later expanded to a more general model in [6]. Numerical schemes based on Split-Up (SU) and the Escalator Boxcar Train (EBT) algorithms were compared in [16], originally proposed in [17] and [18]. While the schemes in [17,18] provide an accurate approximation of the solution, they are only shown numerically to have a convergence rate of first order (i.e. the error is proportional to the mesh size) even when the solution is of smooth density. Higher accuracy is desired in the study of inverse problems and optimal control problems where solving the equation multiple times is required.
In this paper, we present a different approach than the above mentioned numerical schemes which is based on finite difference methods. Such methods have been well studied both in a smooth setting and in an integrable setting [8,9,[19][20][21]. We follow the frame work laid out in the space of integrable functions presented in [9,21] and provide convergence results in the space of Radon measures. While the schemes we present are in the spirit of those presented in [21], we would like to point out some of the key differences of the current work with that presented in [21]. First, the model studied in this paper is different than the model studied in [21] as here we do not consider the type of hierarchical dependency (which is a more general form of dependency) on the structure considered in that work. In particular, in [21] they consider a size-structured model on finite domain [0, 1] and they assume the vital rates g, β, d depend on where u represents the density of individuals of size x and time t. They formulate a model for the density u(x, t) in the spirit of (1.1) on the space of integrable functions. Through finite-difference approximation method they establish the well-posedness of a PDE that is governed by Q (not u) in the space of integrable functions by utilizing the compact embedding of the space of bounded variation functions in the space of integrable function to achieve existence of solutions and prove that the limit satisfies an entropy condition to establish uniqueness. They also show that their first-order finite difference approximation converges in the weak* topology to a (unique) limit u that satisfies (1.5) but did not prove that u satisfies the size-structure model governed by u. In that work, the initial condition regularity assumed for u is that of bounded variation and cohorts of individuals with a single point structure x in the form of Dirac measures are not considered (clearly in the present work such initial conditions are permitted). To our knowledge, the formulation of the model considered in [21] on the space of finite signed measures is still an open question. Secondly, the first order scheme presented in [21] is implicit while the scheme presented in this paper is explicit. Here, we extend the convergence results to explicit schemes on the space of measures and to demonstrate the difference of the schemes (explicit and implicit), we compare the computation times and errors in section 5.1 below. Lastly, while the authors in [21] provided numerical examples of a high-resolution second-order method, no stability estimates or convergence results were given for this method. In this paper, we establish convergence of the second order method to the solution of (1.1). This paper is organized as follows: in Section 2, we present notation associated with the space of Radon measures used in this paper and we recall a useful result. In Section 3, we state assumptions needed for the model (1.1) to be well-posed which follow from the results in [15]. In Section 4, we present two numerical schemes. In particular, in Section 4.1 we provide an explicit upwind scheme and prove the scheme converges to the solution of (1.1) (Theorem 4.2). In Section 4.2, we provide a high-resolution second-order scheme in the spirit of the method presented in [21], but formulated in the space of Radon measures. We then prove convergence results in this space for the high-resolution method (Theorem 4.3). In Section 5, we provide numerical examples comparing the explicit first order method developed in this paper, the high-resolution second order method presented here, and the implicit first order scheme studied in [21]. We test our methods against examples covered in [16,21] (see, example 5.1). We also provide different examples of the flexibility of our scheme in handling discontinuous and singular densities (see, examples 5.2 and 5.3). In Section 6, we provide an application where competition for resources between individuals may influence the maximum size of individuals. This example also has the benefit of demonstrating the regularity of the stationary-state solution discussed in [22]. Finally, in Section 7 we provide some concluding remarks.

Notation and background
We employ the following notation for commonly used spaces: R + = [0, ∞), C k (X; Y) denotes the space of functions from X to Y that are differentiable up to order k, C c (X; Y) denotes the space of continuous functions from X to Y that have compact support, C b (X; Y) denotes the space of bounded continuous functions from X to Y, M + (X) denotes the space of all finite non-negative Radon measures on X, and W 1,∞ (X; Y) denotes the usual Sobolev space on X with codomain Y. When talking about spaces of functions, if the codomain Y = R, we will omit Y from the notation above.
Unless otherwise specified, we equipped M + (R + ) with the Bounded-Lipschitz Norm (BL-norm) defined as The BL-norm has also been referred to as the flat norm [23,24], Dudley norm [25,26], and Fortet-Mourier norm [27,28]. Another norm commonly associated with M + (R + ) is the total variation norm (TV-norm) given by It should be noted that while over nonnegative measures they are equivalent as norms, the induced metrics from the BL-norm and TV-norm do not agree. In general, for ν, µ ∈ M + (R + ), In other words, the BL-norm and TV-norm are different on the space of signed measures. For more information on these norms, we refer the reader to [29] and the references therein.
We say a sequence of non-negative Radon measures, (µ n ), converges weakly if there is a µ ∈ M + (R + ) such that for every f ∈ C b (R + ), In M + (R + ), we additionally have that the BL-norm metricizes weak convergence. For more detail, see [15,30]. We make use of the following compactness result which follows from an application of the wellknown Riesz representation and Alaoglu Theorems: Theorem 2.1. For any compact set K ⊂ R + and any C > 0, the set {µ ∈ M + (K) : µ T V ≤ C} is compact for weak convergence.

Model assumptions
We also assume that the model functions satisfy the following: (A1) For any R > 0, there exists L R > 0 such that for all µ i T V ≤ R and t i ∈ [0, T ] (i = 1, 2) the following hold Remark 3.1. As noted in [15], the Lipschitz in time assumptions can be weakened to the functions having a modulus of continuity with respect to time, ω R (|t 1 − t 2 |). We define a solution to (1.1) as follows: Definition 3.1. Given T ≥ 0, we say a function µ ∈ C([0, T ], M + (R + )) is a weak solution to (1.1) if for all φ ∈ (C 1 ∩ W 1,∞ )([0, T ] × R + ), the following holds:

Discretization methods
In this section we present two methods for approximating the solutions of model (1.1). We first present a first-order upwind-based scheme and then we present a second-order high-resolution scheme.
We assume that the initial measure µ(0) ∈ M + (R + ) has a compact support in [0, L]. By standard range of influence arguments, the solution at time T will have support contained in [0, X] where X = L + ζT .
Remark 4.1. The assumption of a compactly supported initial condition is made for numerical computation purposes. In general, this assumption is not needed since the initial condition belongs to M + (R + ), and thus is inherently tight. In this case, one would use the tight version of Theorem 2.1 which can be found in [15].

Upwind scheme
We approximate the initial measure by a linear combination of Dirac Delta masses. More precisely, . Letting j L = L ∆x , we assume additionally that L is chosen large enough such that m 0 j = 0 for all j ≥ j L (it is possible since µ(0) is assumed to be compactly supported).
The main idea of the proposed method is to use a finite difference scheme based on Eq (1.1) to generate coefficients for the following time steps, approximating µ(k∆t) with . We remark that due to the model functions dependency of µ, their discretization at time k depend on every m k j , j = 1, 2, 3, . . . . We present two different explicit schemes. From here on, we additionally assume for all (t, µ) ∈ [0, T ] × M + (R + ) and x ∈ (0, X), The nonnegativity assumption on g is enforced because we are using an upwind-type scheme below. This assumption can be relaxed but in such a case one needs to apply a combination of an upwind and downwind scheme depending on the sign of g (see for e.g., [31]).
In this section, we study the following explicit scheme: which can be rewritten in the form for j = 1, 2, . . . and k = 1, 2, . . . , K − 1. Notice that m k 0 , which approximate the boundary value D dx µ(k∆t), depends on all the m k j , j ≥ 1, and that m k+1 j , the approximation of µ(k∆t + ∆t)((x j−1 , x j ]), is computed explicitely and directly from the knowledge of the result of the previous time step m k 0 , m k 1 , .... Throughout this section, we assume the following Courant-Friedrichs-Lewy (CFL) condition on our mesh: for ζ defined by (A2), The method defined by (4.3) is a simple Upwind scheme which serves two purposes. First, it provides an easily implemented method where numerical computations presented here suggest it is of first order when the initial condition is absolutely continuous (see Section 5). Secondly, many of the proofs for the higher order method are similar to the easier proofs presented in this section. We begin with a simple, but useful result.
Lemma 4.1. Under the assumptions above we have, for every k, µ k ∆x ∈ M + (R + ).
Proof. Since µ k ∆x is a linear combination of Dirac measures with summable coefficients m k j , it suffices to show that for every k, m k j ≥ 0 for all j. Notice that since µ(0) ∈ M + (R + ), we have m 0 j ≥ 0 for all j = 1, 2, . . . and by the positivity assumptions on g k 0 and β, we have m 0 0 ≥ 0. The result then follows from mathematical induction using formulation (4.3) along with the CFL condition (4.4) and the positivity assumptions (4.1).
An immediate advantage of Lemma 4.1 is the following: In the next lemma, we show that the approximations µ k ∆x have compact support for every k. The same argument can be used to show the approximations are tight for non compactly supported initial conditions.
More precisely m k j = 0 for all j ≥ j L + k.
Proof. We will show that for every k there is an index j k such that ∞ j= j k m k j ≤ 0. Since the m k j are non-negative, m k j = 0 for all j ≥ j k . We accomplish this through induction. We claim for each k, By our assumption on L, the claim holds for k = 0. Assuming the claim is true for k, then from (4.3), the CFL condition (4.4), "telescoping" sums, and the fact that m k In view of this result, we restrict the approximation to a finite sum. That is Then there is a constant C * independent of ∆x and ∆t such that for every k = 0, 1, . . . , K, Notice the "telescoping" term in (4.5), the second equation in scheme (4.2), and the CFL condition We next show the regularity of the scheme.
Lemma 4.4. There exists an L > 0 independent of ∆x and ∆t such that for l > p, Using "summation-by-parts" and the second equation in (4.3), Using that |φ j | ≤ 1 and where we used Lemma 4.3.
Our main goal is to use Ascoli-Arzela's Theorem to guarantee the existence of a convergent subsequence of our scheme. However, the theorem takes place on the set of continuous functions from a compact Hausdorff space into a metric space (for more detail, see [32]) and we have only provided discrete measures. We address this by defining a family of curves µ ∆t ∆x : where χ E (t) denotes the characteristic function over the set E and where ∆t and ∆x satisfy (4.4). One can check that µ ∆t ∆x ∈ C([0, T ], M + ([0, X])) and, with the assistance of Lemmas 4.3 and 4.4, the following two lemmas hold: Here the constant C * is the same from Lemma 4.3.
Lemma 4.6. There exist an L > 0 independent of ∆x and ∆t such that for all t 1 , is precompact for convergence in the BL-norm as well. Owing to Ascoli-Arzela, we obtain the following Theorem.
We essentially follow the proof presented in section 12.5 of [20], while enjoying the nice properties of Dirac measures. Since each curve is determined by its values on the points (t k , x j ), we restrict ourselves to the numerical scheme (4.2). For the moment, take φ ∈ ( We may assume φ is compactly supported as Lemma 4.2 shows the measures µ ∆t ∆x are also of compact support. We multiply (4.2) by φ k j ∆t = φ(t k , x j )∆t and sum over j = 1, 2, . . . , J and k = 0, 1, . . . , K − 1 to arrive at Starting with the right-hand side of (4.6), notice For convenience when working with the left-hand side, we group the related terms and work with them separately. Making use of "summation by parts" and that m k J = 0 for all k we have Thus, by combining (4.7), (4.8), (4.9), and rearranging terms, Eq (4.6) becomes We will now simplify this expression. To this end we first note that for some θ(x) ∈ (0, 1). Assuming that φ has compact support in [0, T ] × R + , we deduce where o(1) → 0 as ∆t → 0 uniformly in k and x. In the same way where o(1) → 0 as ∆x → 0 uniformly in k and x. Moreover using Lemma 4.3 and recalling that K∆t = T , we have ∆t K−1 k=0 R + o(1)dµ k ∆x (x) = o(1). By the above results and with adding and subtracting some terms, we obtain where o(1) −→ 0 as ∆t, ∆x −→ 0. We now verify that for any k and any t ∈ [t k , t k+1 ] there holds Moreover in view of assumption (A1) and Lemma 4.4, we see Finally, using assumption (A2) and Lemma 4.4, we have Since φ ∈ C 2 ([0, T ] × R + ), we have that ∂ x φ(t k , ·) W 1,∞ ≤ C, and so we obtain From the above calculations, we deduce (4.11). Integrating (4.11) for t ∈ [t k , t k+1 ] and summing over k yields The other terms in the right-hand side of (4.10) can be handled in an analogous way. We then deduce from (4.10) that Passing the limit as ∆t, ∆x −→ 0 along a converging subsequence, we then obtain that Eq (3.1) holds for any φ ∈ (C 2 ∩ W 1,∞ )([0, T ] × R + ) with compact support. A standard density argument shows that Eq (3.1) holds for any φ ∈ ( Borrowing uniqueness of the solution of Eq (3.1) from [15], we have in fact that the whole sequence µ ∆t ∆x converges to the solution of (1.1).

Second order high-resolution scheme
In this section, we aim to provide a high-order method. We will provide this approximation over a finite domain as special care is needed around boundary points. We approximate the initial condition as follows: where, the coefficients m 0 j are generated as in Section 4.1. We generate the coefficients with a second-order high-resolution scheme similar to those studied in [9,21]  where the flux term is given by 13) and the sum on the boundary is given by Here we denote by mm(a, b) the minmod function given by mm(a, b) = sgn(a) + sgn(b) 2 min(|a|, |b|) and we denote by ∆ + m k j and ∆ − m k j the forwards and backwards difference in the variable x (index j), respectively. We impose the following CFL condition on the scheme: for ζ defined in (A2), In a similar fashion to the work cited before, it is useful to define the following coefficients:

Scheme (4.12) can then be rewritten as
.

(4.15)
From this formulation and the CFL condition (4.14), the following Lemmas are immediate. The proofs follow the arguments presented in the related Lemmas presented in Section 4.1.
Following the ideas from the Section 4.1, we next prove that our approximations are uniformly bounded and Lipschitz in time.
Proof. Here it is more useful to use formulation (4.12) of the high-resolution scheme. We have As in Lemma 4.3, the "telescoping" term in (4.16), the second equation in scheme (4.12), and the CFL condition (4.14) imply Lemma 4.10. There exists an L * > 0 independent of ∆x and ∆t such that for l > p, Proof. Let φ ∈ W 1,∞ (R + ) be with φ W 1,∞ ≤ 1 . It follows from the same arguments presented in Lemma 4.4 that Notice from definition (4.13), f k Thus, we have As with the upwind scheme, after interpolation, Theorem 2.1 and Ascoli-Arzela's Theorem guarantee the existence of a convergent subsequence of the scheme. Which leads us to the following , we arrive at the following Notice, the only terms different in (4.17) from (4.6) in Section 4.1 are the ones involving the flux. Therefore, we only need to address these terms and follow the argument presented in Theorem 4.2. Using "summation by parts" we have At the moment, we have only provided a scheme that is second order in space only (except at the boundary points). To make the scheme second order in time, we use the following second order in time total variation diminishing Runge-Kutta time discretization studied in [33] which preserves all convergence results. Let

Numerical examples
In this section, we test our schemes against a variety of problems. We compare the performance of the explicit first order method developed in this paper, the high-resolution second order method presented here, and the implicit first order scheme studied in [21]. We do this considering three specific problems with known solution varying the regularity of the inital condition. More specifically, we first consider a problems whose solution is a smooth function, then a problem whose solution is an irregular L 1 function, and eventually a problem whose solution is highly singular being the sum of two Dirac masses. For the examples whose solutions are absolutely continuous, we present the error and the order of the three methods (first order implicit, first order explicit and second order) and also record the execution times.
However, since computing the flat metric exactly is not a simple task, we use the same metric presented in [34] to approximate the BL distance. The metric is given by where W 1 is the 1-Wasserstein distance calculated with the algorithm presented in [34]. We approximate the order of accuracy, q, with the standard calculation: where µ represents the exact solution of the examples considered. As proved in [34], when the support of µ and ν are contained in a compact set, there exists a constant C such that Thus, the order of accuracy is maintained between the two metrics. For the second order scheme, local truncation error analysis of the flux, f k j+ 1 2 , suggests it is more accurate to adjust the mesh size for the interval (x 1 , x 2 ] from ∆x to 3 2 ∆x and the mesh size for the interval (x J−2 , x J−1 ] from ∆x to 1 2 ∆x without changing the values of x j (see, e.g., [8]). We implement these changes in our numerical computations below for the second order method. These changes do not affect the convergence results shown in Section 4.2.

Smooth density
Here we consider a problem where the initial measure is generated by a smooth density function. The purpose of this example is to show that the scheme can achieve second order accuracy under sufficiently nice conditions. We set where x ∈ [0, 1]. One can check that the solution to this problem is given by µ(t) = exp(t − x)dx. In Table 1, we present the error of our method from approximating the solution at time T = 1. We show in Figure 1 the running time of our algorithms.   Table 1 for the explicit upwind (red '− • −'), implicit upwind (purple ' * '), and high-resolution (blue '− −') schemes in a log-log plot.
As we observe in Table 1, the explicit upwind scheme is considerably faster than the implicit upwind scheme developed in [21] when the mesh size satisfies the CFL condition 4.4. For example, to achieve an error of 4.1 × 10 −4 , the explicit upwind scheme takes about 0.7 seconds while the implicit upwind scheme takes about 270 seconds. Meanwhile, the higher-resolution scheme out preforms both methods, achieving an error of about 3.2 × 10 −4 in around 0.11 seconds.
Following [6], we test the influence of nonlocal terms in the model function with the example where µ(0) = (1 + 0.5 cos(x))dx and x ∈ [0, 1]. The exact solution is given by µ(t) = exp(−t)(1 + 0.5 cos(x))dx. We present the BL-error and computation time of both methods in Table 2. Once again we observe that the Upwind scheme has order 1 and that the high-resolution scheme has order 2. As before, in Figure 2 we plot the error of the methods against computation time.   Table 2 for the explicit upwind (red '− • −'), implicit upwind (purple ' * '), and high-resolution (blue '− −') schemes in a log-log plot.

L 1 densities
In this section, we demonstrate the accuracy of both methods for discontinuous densities. It has been shown that the upwind scheme is of order 1 2 in the L 1 norm [35]. In the setting of probability measures, the upwind scheme was shown to be of order 1 in the Wasserstein metric [36]. Below, we present evidence that this holds true for the BL-norm in the setting of Radon measures. We use the model functions g(x) = 1, d(x) = 0, β(x) = 0, and µ(0) = χ [0.5,1] (x)dx (5.4) whose solution is µ(t) = χ [0.5+t,1+t] (x)dx, and test the accuracy of the methods in Table 3. Since d = 0 and β = 0 we are dealing with a conservative equation in the sense that the total mass is preserved in time. So in that case the distance ρ defined in (5.1) is equivalent to W 1 . In view of the results in [36], we thus expect the upwind schemes to have order 1, which is consistent with the results displayed in Table 3. This is in contrast to the 1 2 order achieved by the upwind scheme in the space of integrable functions [35]. As expected the high resolution scheme still out performs both upwind schemes. However, the numerical results suggests that the high resolution scheme does not achieve a second order in this case. Instead, these results suggest that the order of convergence is around 4 3 .

Singular measures
A common phenomena with second order methods when discontinuities and singularities are expected is dispersion [20]. Since our method is based on flux limiter techniques, it behaves well when the density functions are discontinuous or singular. To demonstrate this, we present the following simple example: g(x) = 1, d(x) = 0, β(x) = 0, and µ(0) = δ x=0.5 + δ x=1. 5 (5.5) which has exact solution µ(t) = δ x=t+0.5 + δ x=t+1.5 . We present the approximate solutions below by measuring intervals of length 0.04 (the aggregation intervals). Since in the previous examples the implicit upwind method had similar errors as the explicit upwind (albeit being much slower computationally), we omit showing the implicit upwind scheme results for this example for clarity of the figures. In Figure 3, we show how the two schemes behave in the presence of singular measures. The error and order of the method for this problem is presented below. From the data in Table 4, we observe a drop in order in both methods. This seems to correlate with the drop of order discussed in [35] where it was proven that the optimal L 1 -order for the explicit upwind scheme is 1 2 for nonsmooth initial conditions in L 1 . Analogously, the numerical results in Table 4 also suggest order 1 2 for the upwind scheme when the initial condition is not an absolutely continuous measure. Finally, from these results we observe that the higher-resolution scheme still outpreforms the upwind scheme with an order of approximately 2 3 .

Application
In this section, we present examples of an interesting phenomenon that can be handled by our numerical scheme. In particular, it is well known that availability of resources in the environment may impact the maximum size of the individual due to competition [37,38]. We utilize the model (1.1) and our high-resolution second order numerical scheme (4.12) to understand the dynamics of the population and the maximum size of the individual under such a scenario.  For this example, the birth and death functions depend solely on the total size of the population. The purpose of this is to simulate an environment which dictates a maximum size for individuals due to limiting resources. To this effect, we choose an environment that grows harsher as the population size increases. In this particular example, we observe the population building up at a critical size over long time. This critical size plays the role of a new maximum size different from the initial data. Letting P(t) = µ(t) T V , we take , with µ(0) = e −x dx for x ∈ [0, 1]. In Figure 4, we mesh the measure over the time interval [0,10] aggregated at each step over intervals of length 0.01. Notice that β(P) = d(P) precisely when P = 5.4. These choices of model parameters allows for a logistic-type dynamics for the population and hence population growth is inhibited for large populations due to limiting resources. In fact, since the birth and death rates depend only on the total population, we can observe the change of the total population according to the following differential equation: d dt P(t) = (β(P(t)) − d(P(t)))P(t).  This phenomenon is demonstrated by the mesh in Figure 4. Over time, the population builds up at a particular size. This critical size can be calculated using the model functions and the equilibrium P * = 5.4 to be x * = 5 17 . This build up causes the formation of a shock and thus the classical smooth framework is insufficient. However, this case is easily handled in the setting of L 1 densities and Radon measures.
From results proven in [22], taking the mortality function, d, strictly positive guarantees the stationary state, µ ∞ , is absolutely continuous with respect to the Lebesgue measure regardless of the regularity of µ(0). As t → +∞ we expect this stationary state to have the equilibrium mass of Eq (6.1), that is 1 0 d µ ∞ = P * = 27/5. Then g(P * )(x) = 0 for x ≥ x * where x * is such that x * f (P * ) = 1 i.e., x * = 5/17 0.29. Thus the mass in [x * , 1] as t → +∞ is flowing at a diminishing rate and at the same time is disappearing due to death. This implies µ ∞ ((x * , 1]) = 0 . We can then calculate the stationary state by looking for a solution in the form µ ∞ = u ∞ (x)dx for x ∈ [0, x * ] and 0 otherwise.

Singular initial condition
In this section, we demonstrate the phenomenon recently proven in [22] and discussed in the previous example. We take the same model functions as before, but with singular initial condition given by µ(0) = δ 0.1 + δ 0.5 + δ 0.75 , and we present the solution of the second-order finite-difference scheme in Figure 6. We observe the stationary solution is identical to the previous case; likewise, the formula of u ∞ (x) is independent of the initial condition. This demonstrates that the non-local nature of the model (1.1) has a "smoothing" effect on the regularity of the solution. Furthermore, this example shows that cohorts beyond the (long-term) maximum size, determined by the environment to be 5/17, eventually vanish and all remaining individuals are those with sizes below the maximum size.

Conclusion
In summary, the setting of Radon measures has many benefits to the study of structured population models. One major benefit is the mathematical study of the evolution of cohorts in the form of Dirac measures. Cohorts are a common tool used in biology to model populations of concentrated size. This type of distribution is not readily handled in the setting of smooth or L 1 functions, but is naturally integrated in the setting of Radon measures. When some additional structure such as selection-mutation or hierarchical competition it has been observed that measure solutions form even when the initial condition is smooth. In the case of selection-mutation, measure solutions have been shown to form in infinite time [4,39]. In the case of hierarchical competition, it has been shown that measure solutions will form in finite time [21,40]. These examples show a need for the study of population models in measure spaces.
In this paper, we have provided two schemes which approximate the solution to model (1.1) in the BL-norm. We have provided several numerical examples of these methods demonstrating their accuracy with different regularity on the initial conditions. The high-resolution scheme is shown to be of higher order than the upwind scheme throughout these regularity conditions. High-order and fast methods are critical for parameter estimation problems which involve minimizing a cost functional that requires solving the model (1.1) many times until a minimizer is realized. We remark that our proof of convergence also provides an alternative proof of the existence of solutions to the model (1.1) to those presented in [6,15].
We have also provided an example of an interesting biological phenomenon where the environment determines a maximum size for the population due to competition for limited resources. These examples are also interesting from a mathematical perspective as they show a change in regularity in both a "roughening" and "smoothing" manner to the initial distribution. These examples demonstrate the asymptotic behavior and regularity of the stationary solution studied in [22].
Future work includes rigorously proving the order of accuracy of these methods as well as extending these methods for other biologically relevant models discussed in [4]. Specifically, this scheme can be extended and used to study sensitivity analysis as in [41][42][43].