Finite Range Method of Approximation for Balance Laws in Measure Spaces

In the following paper we reconsider a recently introduced numerical scheme. The method was designed for a wide class of size structured population models as a variation of the Escalator Boxcar Train (EBT) method, which is commonly used in computational biology. The scheme under consideration bases on the kinetic approach and the split-up technique - it approximates a solution by a sum of Dirac measures at each discrete time moment. In the current paper we propose a modification of this algorithm, which prevents (possible) exponential growth of the number of Dirac Deltas approximating the solution. Our approach bases on the finite range approximation of a coefficient which describes birth processes in a population. We provide convergence results, including the convergence speed. Moreover, some results of numerical simulations for several test cases are shown.


Introduction
The main purpose of this paper is to present a numerical scheme designed for a class of one-sex structured population models, which corresponds with a current trend basing on a kinetic approach to the population dynamics problems [2,3,14,17,18,19,22].This approach is a convenient starting point for both analytical and numerical studies and certainly is coherent with empirical data, since a result of an experiment or a measurement in the population dynamics is usually a number of individuals which state is within a specified range.A broad group of mesh-free methods originated from the kinetic theory are particle methods, which are frequently used in models describing behaviour of large groups of interacting particles or individuals.They were originally designed for conservative problems and successfully applied for solving numerically kinetic models originated from physics, see [20,21,23,24] and references therein.In particular, the particle methods were adopted to the Euler equations in fluid mechanics [27], isentropic Euler equations [5,15], Vlasov equation in plasma physics, Boltzmann equation, or Fokker-Planck equation.Recently, they are also used in problems related to the crowd dynamics, pedestrians flow [14,22], collective motion of large groups of agents [7] and population dynamics [8].
A challenge associated with an application of the kinetic theory in the population dynamics models is a nonconservative character of these problems.Furthermore, new particles may appear at arbitrary points of the ambient space.This is the main reason why natural distances for probability measures, like Wasserstein distances, cannot be exploited.Indeed, let µ, ν be finite Radon measures on Ê + such that µ(Ê + ) = ν(Ê + ).Then, according to [25,Definition 6.1] and thus W 1 (µ, ν) = +∞.The problem has been recently overcome by using the flat metric, which is obtained by restriction of the class of test functions in (1) to the space of Lipschitz and bounded functions.Within the latter setting, results concerning well posedness of the following nonlinear size-structured equation were obtained in [8] (see also [18,19]).Variables t ∈ [0, T ] and x ≥ 0 denote, respectively, time and the size of an individual, b, c, η are vital functions, and µ is a finite, nonnegative Radon measure describing a distribution of the individuals with respect to the variable x.Here, b(t, µ) describes a dynamics of a transformation of the individual's state.More precisely, the individual changes its state according to the following ODE ẋ = b(t, µ)(x).
By c(t, µ)(x) we denote a rate of evolution (growth or death rate).The integral on the right hand side of (2) accounts for an influx of the new individuals into the system.In case all new born individuals have the same physiological state x b , then and the integral in (2) transforms into a boundary condition.Thus, the equation ( 2) can be reduced to the following classical nonlinear renewal equation with the nonlocal boundary condition where D λ µ(x b ) is the Radon-Nikodym derivative of µ with respect to the Lebesgue measure at x b .Studying the population dynamics problems in the measure setting is a relatively new approach and thus, formal convergence of particle-based schemes for (2) was difficult to establish for a long period of time.One of the first steps in this direction was made in [18,19], where existence, uniqueness and stability of solutions to (5) in the space of finite, nonnegative Radon measures equipped with the flat metric were proved.Within the latter framework a proof of convergence of the Escalator Boxcar Train (EBT) method was conducted in [4].The EBT scheme applies to the classical nonlinear renewal equation (5).A concept of this method bases on approximating the solution by a sum of Dirac measures, each one of which represents the average state x i and the number of individuals m i within the i-th group (cohort).m i changes its value due to the process of evolution (growth or death), while x i evolves according to the ODE (3).Although the EBT method was described for the first time in 80's in [11], the proof of convergence [4,16] and the convergence rate [16] is very recent.
For the general size-structured equation (2), a numerical scheme based on the splitting technique and the particle methods was developed in [9].Similarly as in the EBT case, the output of the splitting-particle method can be understood as a description of collective behaviour of the individuals divided into cohorts.The main idea of the splitting technique is to separate the dynamics induced by the transport operator from the dynamics induced by the nonlocal term.Formally, the algorithm bases on representing a semigroup generated by the solution as a product of two semigroups related to simpler equations.One of the essential assumptions in [9] is the particular form of the η function, namely which means that an individual at the state y gives rise to offsprings being at the states {f p (y)}, p = 1, . . ., r.A use of the numerical scheme presented in [9] implies that, in most cases, the number of Dirac Deltas approximating the solution increases, in the worst case, exponentially.For instance, setting r = 1 and f 1 (y) = 1 2 y in (6), which is the case of the equal cell division model, implies that the number of Dirac Deltas is doubled at each time step.In order to prevent this phenomenon authors proposed a reconstruction procedure, which is an appropriate approximation of a sum of Dirac Deltas by a sum consisting of a smaller number of Dirac measures.Nevertheless, the reconstruction has to be performed once per several time steps, which may influence the accuracy of the numerical solution.
The essential difference between the scheme presented in [9] and the one developed in the present paper lies in a different method of approximation of the η function.More precisely, we approximate f p functions in ( 6) by piecewise constant functions.As a consequence, new particles appear only at these points of the ambient space, which belong to the set of values of the approximating functions.This procedure prevents the exponential growth of particles, so the reconstruction procedure introduced in [9] is no longer needed.The main problem with such approximation is that the class of piecewise constant functions is not a subclass of the Lipschitz functions, therefore assumptions on the regularity of coefficients from [9] are not fulfilled.In this paper we develop the estimate which does not require Lipschitz continuity of f p functions (see Theorem 2).Our approximation is valid for piecewise constant functions satisfying the condition f (x) ≤ x.It is worth to mention that this limitation is not very restrictive, since the sublinearity of f is a typical assumption in variety of fragmentation models.
This paper is organised as follows.In Section 2 we describe the numerical scheme and the method of finite range approximation.Section 3 consists of some basic facts about the space of finite, nonnegative Radon measures equipped with the flat metric.We also justify considering (2) in the latter space instead of the Banach space (W 1,∞ ) * .In Section 4 we prove convergence of the algorithm and provide the convergence rate.Section 5 is devoted to the results of numerical simulations for several test cases.

Splitting-Particle Method with Finite Range Approximation
The numerical scheme described in this section is a modification of the method developed in [9].It bases on the splitting algoritm, that is, solving subsequently sub-problems associated to the transport and integral operators.As it has already been stated above, the essential difference between both methods lies behind the approximation of the η function.In the setting discussed here, the approximation ensures at most linear growth of the number of Dirac Deltas approximating the solution.In order to keep the paper self-contained, we briefly recall the method and describe in details the new element, that is, the finite range approximation.
The next lemma shows how to construct a function f ε p with the finite image, which approximates accurately a Lipschitz continuous function f p appearing in (6).
Lemma 1.Let ε > 0, M > 0, p = 1, . . ., r, and f p : Ê + → Ê + be a Lipschitz function satisfying the condition f p (x) ≤ x.Then, for each p = 1, . . ., r there exists a function f ε p with a finite range, f ε p : Ê + → Ê + , such that Proof of Lemma 1.Since , for j = 1, . . ., J and p = 1, . . ., r.It follows directly from the construction that ∪ J j=1 A jp = [0, M).Since f ε p is supposed to be defined for all x ∈ Ê + , we redefine the set A Jp := A Jp ∪ [M, +∞).The approximation f ε p is given by the following formula where x j * = (j−1)ε and χ A jp is the characteristic function of the set A jp .Conditions postulated in (7) are thus straightforward consequences of the definition of f ε p .
Remark 1.For our purposes, it is sufficient to consider a finite interval [0, M) in Lemma 1.Note that the equation admits a finite propagation speed property.Assume that the support of the initial data µ(0) is contained in the interval [0, M o ], for some M o > 0.Then, the support of a solution to (9) on a finite time interval [0, T ] is a subset of the interval [0, M), where M depends on M o , an appropriate norm of b, and the length of the time interval [0, T ].Substituting the zero right hand side of (9) by η, where η is given by (6) and f p (x) ≤ x, for p = 1, . . ., r, does not affect the support of the solution.

Description of the algorithm
Fix T > 0 and N ∈ AE.Define the length of a time step ∆t = T /N and a set of discrete time points t k = k∆t, where k = 0, . . ., N. Assume that a numerical solution µ k at time t k is a sum of N k Dirac Deltas, that is, where m i k is a mass of the i-th Dirac Delta and x i k denotes its location.In particular, we assume that the initial data is a sum of Dirac measures.If it is not the case, the initial data can be approximated by such a sum (see Lemma 3 for details).Let M > 0 be a constant, for which the support of a solution to (2) is contained in the interval [0, M), for all t ∈ [0, T ].Such a constant exists according to Remark 1.A procedure of obtaining the numerical solution µ k+1 at time t k+1 is divided into two steps described below.
Step 1.The first step is to calculate new locations of Dirac Deltas.It is obtained by solving the following ODEs system Step 2. The second step is to determine locations of Dirac Deltas which correspond to the new particles and to recalculate masses of all Dirac measures.In order to determine the locations of the new particles we use the finite range approximation of f p , for p = 1, . . ., r, provided by Lemma 1.According to (8), new particles appear at {x j * } J j=1 and these locations are fixed in time.As a consequence, we obtain a set of N k+1 = N k + J Dirac measures.In odder to unify the notation we set x N k +j k = x j * , for j = 1, . . ., J. The resulting system of ODEs, which has to be solved on the time interval [t k , t k+1 ] is the following where c k = c(t k , μk ), β pk = β p (t k , μk ) and After a reassignment of indexes (e.g., arrange all Dirac measures in ascending order with respect to location) we obtain a measure µ k+1 = N k+1 i=1 m i k+1 δ x i k+1 , which is the output of the algorithm at time t k+1 .The essential feature of this algorithm is that after k time steps the number of Dirac Deltas approximating the solution is equal to N o + Jk, where N o is the initial number of Dirac measures.
As it has already been mentioned in the Introduction, well posedness of ( 2) is established for functions f p which are Lipschitz.Clearly, the functions f ε p constructed in Lemma 1 do not satisfy this assumption.Nevertheless, we can construct functions f ε p , which depend only on the initial data approximation µ o and coefficients b, c, β, such that the result of solving ( 10)-( 11) is the same for f ε p and f ε p .This is the subject of the next remark.
Without loss of generality we assume that D ⊂ [0, M), where M is the same constant as in Remark 1.Let f ε p be defined by (8), for p = 1, . . ., r.For a fixed p, f ε p is a piecewise constant function, which can be rewritten in the following form where i=1 is a finite, strictly increasing sequence with d J D +1 = +∞.Define and a piecewise linear function Remark 3. Note that the solutions to the systems of ODEs (10) and (11) are in fact the solutions (according to the definition [8, Definition 2.2]) to the following PDEs and 3 The Space of Measures (M + (Ê + ), ρ F ) Henceforth, M + (Ê + ) denotes the space of nonnegative Radon measures with bounded total variation on Ê + := {x ∈ Ê : x ≥ 0}.We equip M + (Ê + ) with the flat metric 13) can be replaced by W 1,∞ (Ê + ) through a standard mollifying sequence argument applied to the test function ϕ, as its derivative is not involved in the value of the integral, which implies that ρ F is the metric dual to the • (W 1,∞ ) * distance.Note that in this paper, the space M + (Ê + ) is equipped with the metric ρ F and this shall remain until said differently.The space (M + (Ê + ), ρ F ) is complete and separable.Using the standard mollification procedure, ρ F can be equivalently rewritten as and therefore, for all This equality gives a rise to a question about the possibility of setting the model in the Banach space W 1,∞ , • (W 1,∞ ) * .It is a natural question, considering the previous works in which the problem of continuity of solutions with respect to time were addressed.In [26] the author established the continuity in the L 1 topology of solutions to structured population models set in the space of densities.Next it was proved by Diekmann (see [13]) that solutions to balance laws in the space of measures are continuous in weak-* topology of the Radon measures space.Unfortunately, it turns out that the answer to this question is negative, because shift operators are not continuous in W 1,∞ , • (W 1,∞ ) * contrary to the L 1 topology.This is the obstacle one cannot overcome, since the continuity of the shift operators is essential for obtaining the continuity of solutions to the transport equation.
Lemma 2. One-parameter semigroup {T } t≥0 of the shift operators is not a strongly continuous semigroup on the space (W 1,∞ (Ê + )) * .Moreover, for all t > 0 it holds that where In what follows we construct Ψ λ and φ, such that Let δ 1 be a Dirac measure in x = 1, i.e. the functional on the space of bounded continuous functions (BC . This functional can be extended to the linear space ) and f has one-sided limits in x = 1} ⊂ L ∞ (Ê + ), where B(Ê + ) denotes the space of measurable and bounded functions on Ê + .This extension is obtained by putting where 0 ≤ λ ≤ 1.It holds that The Hahn-Banach extension theorem guaratees existence of an extension of δ λ 1 to L ∞ (Ê + ), which will be denoted as δλ 1 .Therefore, setting yields a family of bounded linear functionals Clearly, f ∈ L ∞ (Ê + ), φ ∈ W 1,∞ (Ê + ) and φ W 1,∞ = 1.Substituting Ψ λ defined by (15) and φ defined as above into ( 14) yields Taking supremum over λ ∈ [0, 1] finishes the proof.
Here, BC α,1 ([0, T ] ×M + (Ê + ); W 1,∞ (Ê + )) is the space of W 1,∞ (Ê + ) valued functions which are bounded in the • W 1,∞ norm, Hölder continuous with exponent 0 < α ≤ 1 with respect to time and Lipschitz continuous in ρ F with respect to the measure variable.This space is equipped with the • BC α,1 norm defined by where To simplify the notation, we define x i X where x = (x 1 , . . ., x n ) ∈ X n .
Since in our case η has the specific form (6), we rewrite the original estimates from [8, Theorem 2.11] in terms of β p and f p instead of η.Regularity of β p and f p imposed in (16) guarantees that η defined by ( 6) fulfills the assumptions of [8, Theorem 2.11] and thus, ( 2) is well posed.to (2).Moreover, the following properties are satisfied: 1.For all 0 ≤ t 1 ≤ t 2 ≤ T there exist constants K 1 and K 2 , such that Here, C 2 and C 3 depend on the coefficients (b i , c i , β i , f i ) for i = 1, 2. In particular, the constants depend on the Lipschitz constants of f 1 and f 2 , which makes this estimate useless for us.Indeed, we need to apply Theorem 1 to where f ε p is the approximation of f p defined by (12), and the Lipschitz constant of f ε p can be arbitrarily large.Nevertheless, it is possible to improve (18) and, as a consequence, obtain this estimate with the constants depending on the coefficients of one of the equations.Moreover, the norm • BC can be replaced by the weaker one.This is the subject of the next theorem.
Theorem 2. Let the assumptions of Theorem 1, claim 2, be satisfied.Then, there exists a constant C, which depends only on the coefficients (b 1 , c 1 , β 1 , f 1 ) and µ o (Ê + ), such that where The proof of Theorem 2 bases on formula (20) (see [6, Theorem 2.9]), which allows to consider equations locally in time.Before we proceed, we introduce some preliminary notions.
Proof of Theorem 2. Since C 1 in (18) depends only on (b 1 , c 1 , β 1 , f 1 ), we may assume that µ 1 (t) and µ 2 (t) are solutions to (2) with the same initial data µ o .As the first step, we define the time dependent functions for j = 1, 2, p = 1, . . ., r. Fix n ∈ AE, define ∆t = T /2 n , t i n = i∆t for i = 0, 1, . . ., 2 n , and approximate b j , c j and β j,p as follows: Note, that on each interval [t i n , t i+1 n ) functions defined above do not depend on t.Therefore, according to [8,Theorem 2.8], solving (2) with coefficients b j n , c j n , and β j,p n on the time interval [t i n , t i+1 n ) yields a Lipschitz semigroup.Call S j,i,n the corresponding semigroup and define the map It follows from the construction that F j,n is a Lipschitz semiflow.Without loss of generality we assume that t = t k n .Substituting F 1,n and F 2,n to (20) yields According to estimates [8, proof of Theorem 2.8] for the linear autonomous problem, it holds that Summing over i = 0, . . ., k − 1 yields According to the estimate from [8, proof of Theorem 2.10, claim iii)] it holds that , where C depends only on (b 1 , c 1 , β 1 , f 1 ).Therefore, there exists a constant C, which depends on the latter set of coefficients and µ o , such that According to [8] (see the proof of Theorem 2.10), the map F j,n t µ o converges uniformly with respect to time to µ j (t), as n → +∞.Therefore, passing to the limit in (22) ends the proof.

Error Estimates in ρ F
The following theorem provides the estimate on the speed of convergence of the numerical method described in Subsection 2.1.
Theorem 3. Let µ be a solution to (2) with the coefficients (b, c, β, f ) and initial data µ o .Let µ k be obtained by solving (10) - (11) with the coefficients (b, c, β, f ε ) and initial data µ δ o , where f ε is the approximation of f defined by (12) and µ δ o is a sum of Dirac Deltas.Then, there exists a constant C, which depends on (b, c, β, f ), µ o , and T , such that where ∆t is the length of a time step, ε is the error of the finite range approximation, i.e. f − f ε L ∞ ≤ ε, and I(µ o ) is the error of the initial data approximation, i.e.Proof of Theorem 3. The proof is divided into several steps.For simplicity, in all estimates below, we will use a generic constant C, without specifying its exact form that may change from line to line.
Step 1: The auxiliary scheme.Let us define the auxiliary semi-continuous scheme, which consists in solving subsequently the problems and where µ k ∈ M + (Ê + ), μk is the solution to (24) at time t k+1 and bk , ck , and ηk are defined as bk A solution to the second equation at time t k+1 is denoted by µ k+1 .Denote by ν k+1 a solution to (25) with ηk defined as r p=1 β p (t k , μk )(y) δ x=fp(y) .
Step 2: Error of the finite range approximation.According to (19), it holds that where C is such that C e Ch ≤ C for all h ∈ [0, T ].
Step 3: Error of splitting.
Let ν(t) be a solution to (2) on a time interval [t k , t k+1 ] with initial datum µ k and parameter functions bk , ck , ηk , where bk is defined by (26), According to [10,Proposition 2.7] and [8,Proposition 2.7], the distance between ν k+1 and ν(t k+1 ), that is, the error coming from the application of the splitting algorithm can be estimated as where C depends on the • W 1,∞ norm of b, c, β and Lip(f p ), p = 1, . . ., r. To Using the Lipschitz continuity of the solution µ(t) (Theorem 2) yields Substituting the latter expression into (34) yields Using this inequality in (33) yields Combining the inequality above with (32) and redefining C leads to Finally, putting together (28), (31), and (35) we conclude that Step 4: Adding the errors.Application of the discrete Gronwall's inequality to (36) yields There exists a constant C * depending only on T such that e Ck∆t − 1 < C * k∆t, for each k∆t ∈ [0, T ].Therefore, we deduce and thus, Since k∆t = t k , the assertion is proved.
In the following lemma we show that any measure ν ∈ M + (Ê + ) can be approximated in ρ F with an arbitrarily small error by a sum of Dirac Deltas.
Lemma 3. Let ν ∈ M + (Ê + ) be such that M ν = Ê + dν = 0.Then, for each δ > 0 there exists K ∈ AE and a measure ν = K i=1 m i δ x i , such that reference solution µ ref was calculated with finite range method with parameters ε = 3.90625 • 10 −05 and ∆t = ∆x = 0.00078125.We recall, that ε denotes the accuracy of the η function approximation, see Lemma 1.According to our leading assumption ( 6), η has the following form η(t, µ)(y) = r p=1 β p (t, µ)(y)δ x=fp(y) .In our case, p = 1 and f p (y) = 1 2 y. ∆t denotes the length of a time step, and (x max − 1 2 x o )/∆x is the number of Dirac Deltas approximating the initial data.It is inversly proportional to the error of the approximation of the initial data, see Remak 4. The error is defined by the following formula Err(T, ∆t, ε) where k is such that k∆t = T and T is the final time.Functions ρ is defined as follows where M µ = Ê + dµ and W 1 is the 1-Wasserstien metric on the space of probability measures, which in the one dimensional case can be calculated with the formula Here, F µ i denotes the distribution function of the measure µ i and ρ is a function equivalent to the flat metric ρ F in the sense that there exists a constant C such that Cρ(µ 1 , µ 2 ) ≤ ρ F (µ 1 , µ 2 ) ≤ ρ(µ 1 , µ 2 ), see [9, Lemma 2.1] for details.The order of the method is given by q := lim ∆t→0 log Err(T,2∆t,ε) Err(T,∆t,ε) log 2 .
In Tables 1 -3 we present results of numerical simulations with different values of ∆t and ε.Parameter ∆x is always equal to the corresponding value of ∆t.
Conclusions of our numerical tests are the following.For a fixed value of ε, the order of convergence tends to zero as ∆t = ∆x << ε, see Table 1.It is intuitively clear, since ε denotes the accuracy of the finite range approximation of η coefficient.Hence, decreasing the time step ∆t and parameter ∆x does not decrease the error, since one of the model coefficient is not sufficiently accurately approximated.When ε is relatively small (10 −3 or 10 −4 ), the order of convergence tends to 1, see Table 2 and Table 3. Summing up, the optional approach is to use ε proportional to ∆t and ∆x, which is consistent with the theoretical error estimate from Theorem 3.

Remark 4 .
The error estimate(23) accounts for different error sources.More specifically, the error of the order O(∆t) is a consequence of the splitting algorithm.The term of order O((∆t) α ) follows from the fact that we solve ODEs with parameter functions independent of time, while b, c and η are in fact Hölder continuous with exponent α with respect to time.The error of the initial data approximation I(µ o ) is inversly proportional to the number of Dirac Deltas approximating µ o .According to Remark 3, I(µ o ) can be arbitrarily small.