FINITE DIFFERENCE APPROXIMATIONS FOR MEASURE-VALUED SOLUTIONS OF A HIERARCHICALLY SIZE-STRUCTURED POPULATION MODEL

We study a quasilinear hierarchically size-structured population model presented in [4]. In this model the growth, mortality and reproduction rates are assumed to depend on a function of the population density. In [4] we showed that solutions to this model can become singular (measure-valued) in finite time even if all the individual parameters are smooth. Therefore, in this paper we develop a first order finite difference scheme to compute these measure-valued solutions. Convergence analysis for this method is provided. We also develop a high resolution second order scheme to compute the measurevalued solution of the model and perform a comparative study between the two schemes.


Introduction.
In this paper we study the following model that describes the dynamics of a size-structured population with hierarchy between individuals: u t + (g(x, Q(t, x)) u) x + m(x, Q(t, x))u = 0, (t, x) ∈ (0, T ) × (0, 1), g(0, Q(t, 0)) u(t, 0) = C(t) + which is referred to thereafter as the environment.The value of the constant α represents the relative competition effect of lower ranks on the individual (in particular, α = 0 implies that individuals with lower ranks do not have any competition effect on individuals with higher ranks and as α increases the effect of such competition increases).For example, in a population of trees competing for sunlight, taller trees have an advantage over shorter ones.An age-structured (i.e., g = 1) version of (1.1) was proposed in [9].Existenceuniqueness and asymptotic behavior of solutions to this age-structured model was studied in [9,13].In [5,14], the existence and uniqueness results were extended to the size-structured populations, while the asymptotic behavior for the size structured model with rates that do not depend on size x was studied in [6].To investigate the asymptotic behavior of solutions, in these articles the nonlocal partial differential equation is transformed into a local one by means of variable change.Then an ordinary differential equation is obtained for the total population P (t) = b a u(t, y)dy.This treatment provides a way to study the asymptotic behavior of the total population P (t) and, in turn, of the solution u(t, x) (see [6,9,13]).The above mentioned analytical methods are difficult to apply to the case where the rates depend explicitly on size x.Thus, numerical methodology remains a valuable tool for understanding the dynamics of such complex structured models.
In [2] it was shown that if ∂ Q g(x, Q) ≤ 0 then a bounded variation weak solution, u, to model (1.1)-(1.2) exists globally in time.However, such a condition is restrictive from the biological point of view as it does not include growth rates exhibiting an Allee effect [1].But, without this monotonicity condition on the growth rate g, one can argue (e.g., [2,14]) that the environment Q become discontinuous in finite time even if all the individual vital rates and the initial data are very smooth, and each discontinuity in Q corresponds to a Dirac measure in u.Therefore, to extend the solution beyond this time, measure-valued solutions have to be considered.Measure-valued solutions for structured models with Q = Q(t) (i.e., no hierarchy α = 1) have been studied [7,10,12] and convergence of numerical approximations have been established [8,11].We point out that there is a key difference between the scrambled competition case α = 1 and the contest competition (hierarchy) case 0 ≤ α < 1.In particular, in the scrambled case the solutions tend to be smooth and have bounded total variation [3] when the model ingredients are smooth, while in the contest case the solutions will develop singularities (point measures) [14,2] even if the model parameters are smooth.This difference presents a mathematical difficulty that usually requires one to rely on totally different techniques to answer the desired question be it related to existence-uniqueness, asymptotic behavior or convergence of approximations.In [4], a vanishing viscosity approach was used to establish the existence of measure-valued solutions for the general model (1.1)- (1.2).This is the only result we know of that treats measure-valued solutions for the general environment Q(t, x) given in (1.2).
The main focus of this paper is the development of finite-difference approximations to compute solutions of the model (1.1)-(1.2).In Section 2 we propose a first order implicit approximation which extends that developed in [3] for the case α = 1.In Section 3 we establish convergence of this scheme to the unique measure-valued solution of the model (1.1)-(1.2) in the weak-star topology.In Section 4 we develop a slightly modified high-order finite-difference approximation similar to that presented in [16] to approximate the measure valued solutions.In Section 5 we present extensive numerical results comparing the two schemes.Concluding remarks are presented in Section 6.

2.
A first-order finite-difference scheme.The following regularity conditions will be imposed on our model parameters throughout the paper: (H1) u 0 (x) is a non-negative function with bounded total variation.(H2) m(x, Q) is non-negative twice continuously differentiable with respect to x and continuously differentiable with respect to Q on [0, 1) × [0, ∞).(H3) β(x, Q) is non-negative twice continuously differentiable with respect to x and continuously differentiable with respect to Q on [0, 1) × [0, ∞).(H4) g(x, Q) is twice continuously differentiable function with respect to x and continuously differentiable with respect to Q, g(x, Q) > 0, x ∈ [0, 1) and The following notation will be used throughout this paper: ∆x = 1 N and ∆t = T M denote the spatial and time mesh size respectively.For x j = j∆x, j = 0, 1, 2, • • • , N and t k = k∆t, k = 0, 1, 2, • • • , M , we denote by u k j and Q k j the difference approximations of u(t k , x j ) and Q(t k , x j ) and we define , and C k = C(t k ).We define the difference operator and the 1 norm by The implicit scheme we propose is given by where 2.1.Constructing a difference scheme for Q: A simple case.For the sake of illustration and convenience to the reader we first assume that the parameters satisfy the following: and Multiplying the first equation of (2.1) by ∆x and summing over i = j + 1, . . ., N and using g k N = 0, we obtain Observing from (2.3) that u (2.4) In this specific case, (2.4) is a first-order implicit upwind scheme for where Clearly solutions to (2.5) could develop a discontinuity in a finite time and thus the corresponding density u = −Q x becomes singular.2.2.Constructing a difference scheme for Q: The general case.For the general case, where we assume that our model parameters satisfy our assumptions (H1)-(H5).Consider with Multiplying the first equation of (2.1) by α ∆x and summing over i = 1, 2..., j and multiplying the same equation by ∆x and summing over i = j + 1, ..., N we readily see that (2.8) In this general case, (2.8) is a first-order implicit upwind scheme for where (2.10) 3. Convergence analysis.If we define ) is equivalently written as the following system of linear equation for Thus, for the remainder of this section we assume that (3.1) holds.Note that (3.1) holds, for example, if ∆x and ∆t are chosen sufficiently small such that ∆x for 0 ≤ j ≤ N and a fixed L 1 ∈ (0, 1).We assume that Lemma 3.1.Assume that L 2 ∆t < 1.We have the following estimate Without assuming the positivity of solutions Lemma 3.1 still holds under the assumption that In fact, multiplying the first equation of (2.1) by ∆x sgn(u k+1 j ) and summing over the indices, j = 1, 2, • • • , N , we obtain where we used the fact that and the index set Jump is defined by for some positive constant L 3 .
To this end we have for some positive constant L 4 .Thus {Q k j } has bounded total variation.Observe that by the smoothness assumptions on g and m there exists a fixed constant ζ such that sup Lemma 3.2.There exists an L 5 > 0, independent of ∆x, ∆t, such that for any m > p ≥ 0, we have Proof.Using the second equation of (2.8) and Lemma 3.1, we have for some positive constant ω.Thus, Define a family of functions {Q ∆t,∆x } by Then, the set of functions {Q ∆t,∆x } is compact in the topology of L 1 ((0, T ) × (0, 1)) and we have the following theorem.
Theorem 3.3.There exists a subsequence {Q ∆ti,∆xi } ⊂ {Q ∆t,∆x } which converges to a BV ([0, T ] × [0, 1]) function Q(t, x) in the sense that for all t > 0 and Proof.Using the results of Lemmas 3.1 -3.2 together with the proof of Lemma 16.7 in [18] page 276 establishes the result.Let for every test function Remark.Following a similar argument used in the proof of Lemma 16.9 on page 280 of [18], it can be shown that the limit of the difference approximations in Theorem 3.3 is a weak solution to problem (2.9).
Next we define entropy solution in the sense of [15].To this end, we define and if there exists a measure zero subset E of [0, T ] such that for all t ∈ [0, T ] \ E the function x → Q(t, x) is defined a.e. and lim In Lemma 3.6 and Theorem 3.7 below we restrict the arguments to the subsequence in Theorem 3.3, i.e., for simplicity of notation we let Q ∆t,∆x = Q ∆ti,∆xi in these arguments.In order to establish that the limit obtained in Theorem 3.3 is an entropy solution we first prove the following auxiliary lemma.
Proof.Using the fact that multiplying the first equation in (2.8) on both sides by Q k+1 j ∆x ∆t and summing (3.5)By Lemma 3.1 and Theorem 3.3 one can easily verify that the sequence Q ∆t,∆x (t, x) constructed via the difference scheme (2.1) converges (subsequentially) to the weak solution Q(t, x) established in Theorem 3.3 in L 2 norm.From this fact and assumption (H4) one can argue that 1 2 Therefore, since δ 3 → 0 as ∆x, ∆t → 0.
The next theorem shows that the limit function Q(t, x) defined in Theorem 3.3 is an entropy solution of (2.9).
Theorem 3.7.The limit function Q(t, x) defined in Theorem 3.3 is an entropy solution of (2.9).
Multiplying the first equation of (2.8) by sign 0 (Q k+1 j − c)φ k j ∆x ∆t, and summing over j = 1, . . ., N and k = 0, . . ., M − 1, we obtain (3.10)We consider each term of (3.10) separately.The first term on the left-hand side of (3.10) satisfies where we used the summation by parts formula given by As for the second term on the left-hand side of (3.10) using Taylor's series expansion on ) where , and Qk+1 Adding and subtracting G(x j , c) and G(x j−1 , c) in (3.12) we get Now adding and subtracting sign 0 (Q since Here the index set Jump1 is defined by We also use the fact that g(x, Q) > 0 and φ has compact support to get the inequality in the entropy solution.
Next we consider the third term on the left hand side of (3.10) Here, we used the fact where (3.18) Consider the right hand term of (3.10), we have from the second equation of (2.8) Here we used the fact that Hence,

FINITE DIFFERENCE APPROXIMATIONS FOR MEASURE-VALUED SOLUTIONS 247
Thus, where where To this end, we define γk+1 (3.25) Then by (3.11), (3.15), (3.16) and (3.25), we have the following inequality Using Lemma 3.6 one can show that as ∆t, ∆x → 0 + .Using similar techniques as in the proof of Lemma 16.10 page 279 in [18] we can show that the subsequence Q ∆t,∆x satisfies for t ∈ [0, T ] and every test function φ ∈ C 1 ([0, T ] × [0, 1]).Here 1 → 0 as ∆t and ∆x → 0. By the bounded convergence theorem as ∆t, ∆x → 0 + the limit Q for any convergent subsequence defined in Theorem 3.3, satisfies the inequality in Definition 3.5.
Similarly, from (2.8) where | → 0 as ∆x, ∆t → 0 + .From this it follows that the limit Q satisfies Definition 3.5 and hence is an entropy solution.
Combining Theorem 3.7 and the uniqueness of the entropy solution for (2.9) established in [4], we obtain the following theorem: Theorem 3.8.The sequence Q ∆t,∆x (t, x) constructed via the difference scheme (2.1) converges to the unique entropy solution Q(t, x) of (2.9).
3.1.Weak star convergence.Our goal for this section is to show that the approximations u ∆t,∆x converges in weak* topology to a unique limit u and establish the relation (1.2) between the limits Q(t, x) and u(t, x).To this end, define a family of functions {u ∆t,∆x } by Theorem 3.9.The sequence u ∆t,∆x (t, x) constructed via our difference scheme converges in weak* topology to a unique limit u(t, x) which satisfies Proof.First we will establish the uniqueness of the measure valued limit u.Since Q(t, x) is unique, suppose that u(t, x) and û(t, x) are two finite radon measures satisfying (3.27), i.e., and Subtracting (3.28) and (3.29), we have The equation (3.30) is true if and only if u = û, since any two such measures agree on all compact sets of the interval [0, 1].So we conclude that u(t, x) is unique. Clearly, Thus, using Lemma 3.1, we obtain by the Banach-Alaoglu-Bourbaki Theorem that there exist a subsequence {u ∆tr,∆xr } of {u ∆t,∆x } and a limit u such that as r → ∞ (i.e., as ∆t r , ∆x r → 0) the following hold: {u ∆tr,∆xr } converges to u in C[0, 1] * uniformly in t, i.e., and for all φ ∈ C[0, 1].Next we will prove the relation between Q(t, x) and u(t, x).To this end, From the definition of the family of functions {u ∆tr,∆xr }, we have (3.33) taking the limit as r → ∞ in (3.33) we obtain From the theory of distributions we have

4.
A second order explicit finite-difference scheme.To approximate the measure valued solutions, in this section we will develop a second order finitedifference approximation similar to that presented in [16].The following notation will be used throughout this section: ∆x = 1 N and ∆t = T M denote the spatial and time mesh size respectively.For x j = j∆x, j = 0, 1, 2, • • • , N and t k = k∆t, k = 0, 1, 2, • • • , M , we will denote by u k j and Q k j the difference approximations of u(t k , x j ) and Q(t k , x j ) and we define We define the finite difference operators We discretize the first equation of the model (1.1) using the following secondorder explicit high-resolution finite difference approximation where the numerical flux is defined by Notice that This finite difference scheme is second-order accurate in space except at the boundary (first and last two grid intervals), where it is first-order accurate.But this provides second-order accuracy in the global L 1 norm.The integral in the boundary condition of the model (1.1) (second equation of (1.1)) is implemented by a second order composite trapezoid rule, except at the first interval which is computed by a first order method using the right-hand point.That is, A similar technique can be used to compute Q k j for 0 ≤ j ≤ N as We now rewrite the scheme (4.1) together with the initial and boundary conditions as the following system of equations: We set It can be easily shown that We now impose the following CFL type condition concerning ∆x and ∆t : where ζ is defined in (3.2).
In this case g ).In Figure 1 the function u(1.5, x) is plotted for several values of N and M using the first-order semi-implicit (FOSI) and second-order explicit (SOE) schemes.Figure 1 indicates that Dirac delta measures are forming at x ≈ 0.61 and x ≈ 0.75.In Figure 2 the three dimensional dynamics of the solution approximated using SOE scheme is presented (where N = 640, M = 6400).In Figure 3 (left) the two dimensional distributions of Q(0, x) and Q(1.5, x) using the FOSI and SOE scheme are presented (where N = 640, M = 6400).In Figure 3 (right) we zoom on the values x ∈ [0.55, 0.8] in the distribution Q(1.5, x) to clearly demonstrate the emerging discontinuities at x ≈ 0.61 and x ≈ 0.75.Since, u = −(1 − α)Q x , then each discontinuity in Q(t, x) corresponds to a point measures in u, as was observed from Figures 1 and 3. Figures 1 and 3 clearly demonstrate that SOE scheme is superior to FOSI in capturing the singular solutions with high accuracy.
To further compare the performance of the SOE scheme and FOSI scheme, we consider a test problem with initial condition u 0 (x) = exp(−x) and vital rates g(x, Q) = 2 − 2 exp(x − 1), β(x, Q) = 2, and m(x, Q) = 1.One can easily verify that for this choice of parameters the equation (1.1) along with boundary condition has an exact smooth solution u(t, x) = exp(t − x).In Table 1, we present errors and order of convergence (in L 1 norm) resulting from the FOSI scheme and the SOE scheme.We compute the L 1 norm of the error by We also calculated the CPU time required for each simulation.From Table 1, it is clear that the SOE scheme requires significantly less computational time than the FOSI scheme to achieve similar accuracy.In particular, we see that to compute an approximate solution with the L 1 norm of the error equal to 0.000087 the FOSI scheme requires 2331.11seconds, but using the SOE scheme we get better accuracy in just 9.61 seconds; thus saving 99.4% of CPU time.In Figure 4 we plotted the L 1 norm of the errors in the FOSI and SOE schemes against the computation time using logarithmic scales.6. Concluding remarks.In this article we have developed a first-order semiimplicit scheme for approximating the measure-valued solution of a hierarchically size-structured population model.This approach provided not only existence of solutions but also a numerical method for computing these measure-valued solutions.Convergence of the first order finite difference approximations to a measure valued solution is established.We also developed a second order high resolution  explicit scheme for approximating the measure-valued solution.In our numerical simulations, we showed that the scheme is non-oscillatory in the presence of discontinuous (or singular) solutions.We compared the first order semi-implicit scheme with the second order explicit scheme.The comparison of these two schemes suggests that although the first order scheme captures the point measures, the second order scheme yields similar accuracy with much less time and space grid points and hence significantly less computational time.
Our future efforts will focus on studying the convergence of the high resolution second order explicit scheme developed in this paper.We think that the general approach used here to study the convergence of first order scheme could be adapted for the second order scheme.However, novel techniques may be needed to obtain estimates analogous to those obtained in Section 3.
Finally, we point out that the analysis may be extended to a more general environment than (1.2) given by Q(t, x) = α

1 0β
(y, Q(t, y))u(t, y)dy, t ∈ [0, T ], u(0, x) = u 0 (x),x∈ [0, 1].(1.1)Here, x ∈ [0, 1], t ∈ [0, T ],and u(t, x) is the density of individuals of size x at time t.That is, b a u(t, x)dx represents the number of individuals having sizes between a and b at time t.The function g represents the growth rate of an individual of size x, while m represents the mortality rate of an individual of size x.The function β is the reproduction rate of an individual of size x and C(t) is the inflow of the smallest size individuals by an external source (for example, in a population of trees C(t) represents the rate at which seeds carried by the wind enter the population).The function Q is given by Q(t, x) = α x 0 u(t, y)dy + 1 x u(t, y)dy, 0 ≤ α < 1, (1.2)

and the 1
norm by ) and the minmod function mm is defined by mm(a, b) = sign(a) + sign(b) 2 min(|a|, |b|).

Figure 1 .
Figure 1.Numerical solutions of u(1.5, x) obtained by the FOSI and SOE schemes for various sets of values of N and M .

Figure 2 .Figure 3 .
Figure 2. The three dimensional distribution of u(t, x) approximated by the SOE scheme using N = 640 and M = 6400.

Figure 4 .
Figure 4.The loglog plot of the L 1 norm of the errors for the FOSI and SOE schemes against the computational time.

Table 1 .
L 1 errors, order of accuracy and computational time (in seconds) for the FOSI and SOE schemes.