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

1.   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:

where

Here, given a Borel subset $A\subset \mathbb{R}^+$, $\mu(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 $\beta$ represents the reproduction rate of these individuals. More precisely, at time $t$ and distribution $\mu(t)$ of the structure, an individual with structure $x$ produces offspring at rate $\beta(t, \mu(t))(x)$. Thus, $D_{dx}\mu (0)$ denotes the Radon-Nikodym derivative of $\mu(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$, $\mu(t)(x)$ informally, change in time $t$ due to the combination of two effects: the transport term $\partial_x(g(t, \mu)\mu)$ which moves the distribution $\mu$ 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, \beta,$ are assumed to be influenced by both time, $t$, and solution to the population model, $\mu(t)$. In applications (e.g., see ^[6,10,11]), it is common to choose $\beta, \; g$ and $\mu$ to depend on a weighted mean of the population in the following form:

for given maps $B, D, G: [0, T] \times \mathbb{R}^+ \times \mathbb{R}^+ \to \mathbb{R}^+$ and $K_B, K_G, K_D : \mathbb{R}^+ \to \mathbb{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

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, \beta, 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.

2.   Notation and background

We employ the following notation for commonly used spaces: $\mathbb{R}^+ = [0, \infty)$, $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$, $\mathcal{M^+}(X)$ denotes the space of all finite non-negative Radon measures on $X$, and $W^{1, \infty}(X; Y)$ denotes the usual Sobolev space on $X$ with codomain $Y$. When talking about spaces of functions, if the codomain $Y = \mathbb{R}$, we will omit $Y$ from the notation above.

Unless otherwise specified, we equipped $\mathcal{M^+}(\mathbb{R}^+)$ with the Bounded-Lipschitz Norm (BL-norm) defined as

where $\|\phi\|_{ W^{1, \infty}} = \max\{ \| \phi\|_{L^\infty}, \| \partial_x \phi\|_{L^\infty} \}$. 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 $\mathcal{M^+}(\mathbb{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 $\nu, \mu \in \mathcal{M^+}(\mathbb{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, $(\mu_n)$, converges weakly if there is a $\mu \in \mathcal{M^+}(\mathbb{R}^+)$ such that for every $f \in C_b(\mathbb{R}^+)$,

as $n \longrightarrow \infty$. We say $(\mu_n)$ is tight if

In $\mathcal{M^+}(\mathbb{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 well-known Riesz representation and Alaoglu Theorems:

Theorem 2.1. For any compact set $K\subset \mathbb{R}^+$ and any $C > 0$, the set $\left \lbrace \mu \in \mathcal{M^+}(K) : \|\mu \|_{TV} \leq C \right \rbrace$ is compact for weak convergence.

3.   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 $\| \mu_i\|_{TV}\leq R$ and $t_i\in [0, T]$ ($i = 1, 2$) the following hold

(A2) There exists $\zeta > 0$ such that

(A3) For all $(t, \mu)\in [0, T]\times \mathcal{M^+}(\mathbb{R}^+)$,

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, $\omega_R(|t_1 - t_2|).$

Remark 3.2. The functions given in (1.3) satisfy assumptions (A1) and (A2) provided that $B, G, D \in W^{1, \infty} ([0, T] \times \mathbb{R}^+ \times \mathbb{R}^+)$ and $K_B, K_G, K_D \in W^{1, \infty} (\mathbb{R}^+)$. Furthermore, Eq (1.4) provides a biologically relevant example of a growth function which satisfies (A1) - (A3).

We define a solution to (1.1) as follows:

Definition 3.1. Given $T\geq 0$, we say a function $\mu \in C([0, T], \mathcal{M^+}(\mathbb{R}^+))$ is a weak solution to (1.1) if for all $\phi \in (C^1 \cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$, the following holds:

4.   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. For fixed $\Delta x$ and $\Delta t$, we discretize $\mathbb{R}^+$ and $[0, T]$ with points $x_j = j\Delta x$ and $t_k = k\Delta t$ where $j = 0, 1, \dots$ and $k = 0, 1, \dots K-1$ such that $K\Delta t = T$.

We assume that the initial measure $\mu(0) \in \mathcal{M^+}(\mathbb{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+\zeta 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 $\mathcal{M^+}(\mathbb{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].

Remark 4.2. The schemes and proofs below hold for the applicable case where the variable $x$ is contained in a finite interval, $[0, X]$, provided we additionally have $g(t, \mu)(X) = 0$ for all $(t, \mu) \in [0, T]\times \mathcal{M}^+(\mathbb{R}^+)$.

4.1.   Upwind scheme

We approximate the initial measure by a linear combination of Dirac Delta masses. More precisely,

where

Letting $j_L = \lfloor\frac{L}{\Delta x}\rfloor$, we assume additionally that $L$ is chosen large enough such that $m_{j}^0 = 0$ for all $j\geq j_L$ (it is possible since $\mu(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 $\mu(k\Delta t)$ with

We let $g_j^k = g(t_k, \mu^k_{\Delta x})(x_j)$, $d_j^k = d(t_k, \mu^k_{\Delta x})(x_j)$, and $\beta^k_{j} = \beta (t_k, \mu^k_{\Delta x})(x_j).$ We remark that due to the model functions dependency of $\mu$, their discretization at time $k$ depend on every $m_j^k$, $j = 1, 2, 3, \dots$. We present two different explicit schemes. From here on, we additionally assume for all $(t, \mu)\in [0, T]\times \mathcal{M^+}(\mathbb{R}^+)$ and $x\in (0, X)$,

Remark 4.3. 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, \dots$ and $k = 1, 2, \dots, K-1$. Notice that $m_0^k$, which approximate the boundary value $D_{dx}\mu(k\Delta t)$, depends on all the $m_j^k$, $j\ge 1$, and that $m_j^{k+1}$, the approximation of $\mu(k\Delta t+\Delta t)((x_{j-1}, x_{j}])$, is computed explicitely and directly from the knowledge of the result of the previous time step $m_0^k, m_1^k, ...$.

Throughout this section, we assume the following Courant-Friedrichs-Lewy (CFL) condition on our mesh: for $\zeta$ 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$, $\mu^k_{\Delta x} \in \mathcal{M^+}(\mathbb{R}^+)$.

Proof. Since $\mu^k_{\Delta x}$ is a linear combination of Dirac measures with summable coefficients $m_j^k$, it suffices to show that for every $k$, $m_j^k \geq 0$ for all $j$. Notice that since $\mu(0) \in \mathcal{M^+}(\mathbb{R}^+)$, we have $m_j^0 \geq 0$ for all $j = 1, 2, \dots$ and by the positivity assumptions on $g_0^k$ and $\beta$, we have $m_0^0 \geq 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:

Corollary 4.1. The total variation $\Vert \mu_{\Delta x}^k \Vert_{TV}$ of the approximated solution is

In the next lemma, we show that the approximations $\mu^k_{\Delta 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.

Lemma 4.2. For every $k = 0, 1, \dots, K$, let $\mathcal{C}_k = [0, j_{L}+k]$. Then,

More precisely $m_j^k = 0$ for all $j \geq j_L+k$.

Proof. We will show that for every $k$ there is an index $j_k$ such that $\sum_{j = j_k}^\infty m_j^k \leq 0.$ Since the $m_j^k$ are non-negative, $m_j^k = 0$ for all $j \geq 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_{j_L +k} = 0$ we have

In view of this result, we restrict the approximation to a finite sum. That is

for some positive integer $J > j_L+K$.

Lemma 4.3. Let $R > 0$ be such that $\Vert \mu(0) \Vert_{TV} \leq R$. Then there is a constant $C^*$ independent of $\Delta x$ and $\Delta t$ such that for every $k = 0, 1, \dots, K$,

Proof. From scheme (4.2) we have

Notice the "telescoping" term in (4.5), the second equation in scheme (4.2), and the CFL condition (4.4) imply

Therefore,

We next show the regularity of the scheme.

Lemma 4.4. There exists an $\mathcal{L} > 0$ independent of $\Delta x$ and $\Delta t$ such that for $l > p$,

Proof. Let $\phi \in W^{1, \infty}(\mathbb{R}^+)$ with $\Vert \phi \Vert_{W^{1, \infty}} \leq 1$. Letting $\phi_j = \phi(x_j)$ then we have

Using "summation-by-parts" and the second equation in (4.3),

Using that $|\phi_j|\le 1$ and $|\phi_{j+1}-\phi_j|\le |x_{j+1}-x_j| = \Delta x$, we obtain

where we used Lemma 4.3.

Our main goal is to use Ascoli-Arzela's Theorem to guarantee the existence of a convergent sub-sequence 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 $\mu^{\Delta t}_{\Delta x}:[0, T]\to \mathcal{M}^+(\mathbb{R}_+)$ linearly interpolating the $\mu^{k}_{\Delta x}$:

where $\chi_E(t)$ denotes the characteristic function over the set $E$ and where $\Delta t$ and $\Delta x$ satisfy (4.4). One can check that $\mu^{\Delta t}_{\Delta x} \in C([0, T], \mathcal{M^+}([0, X]))$ and, with the assistance of Lemmas 4.3 and 4.4, the following two lemmas hold:

Lemma 4.5. Assume $\Vert \mu (0) \Vert_{TV} \leq R$. Then

for all $t \in [0, T]$. Here the constant $C^{*}$ is the same from Lemma 4.3.

Lemma 4.6. There exist an $\mathcal{L} > 0$ independent of $\Delta x$ and $\Delta t$ such that for all $t_1, t_2 \in [0, T]$,

With Theorem 2.1, we have for each $t \in [0, T]$, the set $M(t) = \{\mu^{\Delta t}_{\Delta x}(t) \} \subset \mathcal{M^+}([0, X])$ is precompact for weak convergence. Since $[0, X]$ is compact, this implies $M(t)$ is precompact for convergence in the BL-norm as well. Owing to Ascoli-Arzela, we obtain the following Theorem.

Theorem 4.1. For any $\Delta t, \Delta x \longrightarrow 0$, there exists a subsequence $\left(\mu^{\Delta t_i}_{\Delta x_i} \right)_{i \in \mathbb{N}}$ of $\left(\mu^{\Delta t}_{\Delta x} \right)$, the family of measures defined above, which converges to some $\mu$ in $C([0, T], \mathcal{M^+}([0, X]))$, i.e., $\sup_{t\in[0, T]} \| \mu_{\Delta x_i}^{\Delta t_i} (t) - \mu(t) \|_{BL} \to 0$.

In fact we can prove that

Theorem 4.2. As $\Delta x, \Delta t\to 0$ the sequence $\mu_{\Delta x}^{\Delta t}$ converges in $C([0, T], \mathcal{M^+}([0, X]))$ to the solution of (1.1), i.e., $\sup_{t\in[0, T]} \| \mu_{\Delta x}^{\Delta t} (t) - \mu(t) \|_{BL} \to 0$.

Proof. 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 $\phi \in (C^2 \cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$. We may assume $\phi$ is compactly supported as Lemma 4.2 shows the measures $\mu^{\Delta t}_{\Delta x}$ are also of compact support. We multiply (4.2) by $\phi_j^k \Delta t = \phi(t_k, x_j)\Delta t$ and sum over $j = 1, 2, \dots, J$ and $k = 0, 1, \dots, 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_J^k = 0$ for all $k$ we have

and

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

Moreover

for some $\theta(x)\in (0, 1)$. Assuming that $\phi$ has compact support in $[0, T]\times \mathbb{R}^+$, we deduce

where $o(1)\to 0$ as $\Delta t\to 0$ uniformly in $k$ and $x$. In the same way

where $o(1)\to 0$ as $\Delta x\to 0$ uniformly in $k$ and $x$. Moreover using Lemma 4.3 and recalling that $K\Delta t = T$, we have $\Delta t \sum_{k = 0}^{K-1} \int_{ \mathbb{R}^+}o(1)d\mu^k_{\Delta x}(x) = o(1)$. By the above results and with adding and subtracting some terms, we obtain

where $o(1) \longrightarrow 0$ as $\Delta t, \Delta x \longrightarrow 0$.

We now verify that for any $k$ and any $t\in [t_k, t_{k+1}]$ there holds

where $o(1)\to 0$ as $\Delta t, \Delta x\to 0$ uniformly in $k, t, x$. Indeed since $\partial_x\phi(t, x) = \partial_x\phi(t_k, x)+o(1)$ with $o(1)\to 0$ as $\Delta t\to 0$ uniformly in $t, x$, we have

Moreover in view of assumption (A1) and Lemma 4.4, we see

and so

Finally, using assumption (A2) and Lemma 4.4, we have

Since $\phi \in C^2([0, T]\times \mathbb{R}^+)$, we have that $\| \partial_x\phi(t_k, \cdot) \|_{W^{1, \infty}}\le C$, and so we obtain

From the above calculations, we deduce (4.11).

Integrating (4.11) for $t\in [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 $\Delta t, \Delta x \longrightarrow 0$ along a converging subsequence, we then obtain that Eq (3.1) holds for any $\phi \in (C^2\cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$ with compact support. A standard density argument shows that Eq (3.1) holds for any $\phi \in (C^1\cap W^{1, \infty})([0, T]\times \mathbb{R}^+)$.

Borrowing uniqueness of the solution of Eq (3.1) from ^[15], we have in fact that the whole sequence $\mu_{\Delta x}^{\Delta t}$ converges to the solution of (1.1).

4.2.   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_j^0$ 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

and the sum on the boundary is given by

Here we denote by ${\rm mm}(a, b)$ the minmod function given by

and we denote by $\Delta_+ m_j^k$ and $\Delta_- m_j^k$ the forwards and backwards difference in the variable $x$ (index $j$), respectively. We impose the following CFL condition on the scheme: for $\zeta$ defined in (A2),

In a similar fashion to the work cited before, it is useful to define the following coefficients:

and

Notice, $|A_j^k| \leq \frac{3\Delta t}{2\Delta x}\zeta$ and $A_j^k -B_j^k \geq 0$ as

Scheme (4.12) can then be rewritten as

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.

Lemma 4.7. For every $k = 0, 1, \dots, K$, $\mu_{\Delta x}^k \in \mathcal{M^+(\mathbb{R}^+)}$.

Lemma 4.8. For every $k = 0, 1, \dots, K$, let $\mathcal{C}_k = [0, j_{L}+k]$. Then

Following the ideas from the Section 4.1, we next prove that our approximations are uniformly bounded and Lipschitz in time.

Lemma 4.9. Assume $\Vert \mu(0) \Vert_{TV} \leq R$. Then there is a constant $C^{**}$ independent of $\Delta x$ and $\Delta t$ such that for every $k = 0, 1, \dots, K$,

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

Therefore,

Lemma 4.10. There exists an $\mathcal{L^*} > 0$ independent of $\Delta x$ and $\Delta t$ such that for $l > p$,

Proof. Let $\phi \in W^{1, \infty}(\mathbb{R}^+)$ be with $\| \phi \|_{W^{1, \infty}} \leq 1$. It follows from the same arguments presented in Lemma 4.4 that

Notice from definition (4.13),

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

Theorem 4.3. As $\Delta x, \Delta t\to 0$ the sequence $\mu_{\Delta x}^{\Delta t}$ converges in $C([0, T], \mathcal{M^+}([0, X]))$ to the unique solution of (1.1), i.e., $\sup_{t\in[0, T]} \| \mu_{\Delta x}^{\Delta t} (t) - \mu(t) \|_{BL} \to 0$.

Proof. We follow the same argument presented in the upwind scheme. Multiplying (4.12) by $\phi_j^k \Delta t = \phi(t_k, x_j)\Delta t$ where $\phi \in (C^2_c \cap W^{1, \infty})(\mathbb{R}^+)$, 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

Notice that since $g$ and the $m_j^k$ are bounded, the terms from (4.13), $\frac{1}{2}(g_{j+1}^k - g_{j}^k)m_j^k + \frac{1}{2} g_j^k {\rm mm}(\Delta_+ m_j^k, \Delta_- m_j^k) \longrightarrow 0$ as $\Delta x \longrightarrow 0$. Thus, (4.18) reduces to Eq (4.9) in the upwind scheme with the addition of an $O(\Delta t)$ term which arises from the boundary. The proof follows as in Theorem 4.2.

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

with

and compute

with

5.   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 $\mu$ represents the exact solution of the examples considered. As proved in ^[34], when the support of $\mu$ and $\nu$ 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+\frac{1}{2}}$, suggests it is more accurate to adjust the mesh size for the interval $(x_1, x_2]$ from $\Delta x$ to $\frac{3}{2}\Delta x$ and the mesh size for the interval $(x_{J-2}, x_{J-1}]$ from $\Delta x$ to $\frac{1}{2} \Delta 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.

5.1.   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 \in [0, 1]$. One can check that the solution to this problem is given by $\mu(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.

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 \times 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 \times 10^{-4}$ in around 0.11 seconds.

Following ^[6], we test the influence of nonlocal terms in the model function with the example

where $\mu(0) = (1+0.5\cos(x))dx$ and $x \in [0, 1]$. The exact solution is given by $\mu(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.

5.2.   $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 $\frac{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

whose solution is $\mu(t) = \chi_{[0.5+t, 1+t]}(x)dx$, and test the accuracy of the methods in Table 3. Since $d = 0$ and $\beta = 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 $\rho$ 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 $\frac{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 $\frac{4}{3}$.

5.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:

which has exact solution $\mu(t) = \delta_{x = t+0.5}+\delta_{x = t+1.5}$. We present the approximate solutions below by measuring intervals of length 0.04 (the agammaegation 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 $\frac{1}{2}$ for non-smooth initial conditions in $L^1$. Analogously, the numerical results in Table 4 also suggest order $\frac{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 $\frac{2}{3}$.

6.   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.

6.1.   Population grouping

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) = \| \mu (t) \|_{TV}$, we take

with $\mu (0) = e^{-x}dx$ for $x \in [0, 1]$. In Figure 4, we mesh the measure over the time interval [0, 10] agammaegated at each step over intervals of length 0.01. Notice that $\beta (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:

Therefore, through a standard phase plane argument, one can show that any solution with $P(0) = \|\mu(0)\|_{TV} > 0$ converges to the equilibrium $P^* = 5.4$.

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^* = \frac{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, $\mu_{\infty}$, is absolutely continuous with respect to the Lebesgue measure regardless of the regularity of $\mu(0)$. As $t\to +\infty$ we expect this stationary state to have the equilibrium mass of Eq (6.1), that is $\int_0^1 d\, \mu_\infty = P^* = 27/5$. Then $g(P^*)(x) = 0$ for $x\ge x^*$ where $x^*$ is such that $x^*f(P^*) = 1$ i.e., $x^* = 5/17\simeq 0.29$. Thus the mass in $[x^*, 1]$ as $t\to +\infty$ is flowing at a diminishing rate and at the same time is disappearing due to death. This implies $\mu_\infty((x^*, 1]) = 0$. We can then calculate the stationary state by looking for a solution in the form $\mu_\infty = u_\infty(x)dx$ for $x \in [0, x^*]$ and $0$ otherwise.

Using Eq (6.1), we have $\mu_{\infty}$ satisfies

Since $g(P^*)(0) = 1$ the boundary condition reveals the initial condition

Taking a smooth $\phi$ with compact support in $(0, x^*)$ we have

where we integrated by parts and used that $g(P^*)'(x) = -f(P^*)$. Since this holds for any $\phi$ we see that $u_\infty$ solves the differential equation

whose solution is

with $C = u_\infty(0) = \beta(P^*)P^*$. For convenience, here $P^* = 5.4$, $x^* = 5/17$, $d(P^*) = \beta(P^*) = 9/5$, $f(P^*) = 17/5$, $\beta(P^*)P^* = 9.72$, $\frac{d(P^*)-f(P^*)}{f(P^*)} = -8/17$. For comparison, in Figure 5 we have plotted the formula above against the density from the high-resolution scheme approximation.

6.2.   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 $\mu(0) = \delta_{0.1} + \delta_{0.5} + \delta_{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_\infty (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.

7.   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 ^[15,6].

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].

Acknowledgments

The work of N. Saintier was supported by grant UBACYT 20020170200256BA from the University of Buenos Aires. The authors would like to thank two anonymous referees for their thorough review and useful comments that led to an improved version.

Conflict of interest

The authors declare there is no conflict of interest.